diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index c953718bcaf85061b36629e20abf01cf585e08c1..d819cf4cf2e145a7cf6e6c5f324ce05b448d4bfc 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -15,47 +15,36 @@ */ #include "CbmL1.h" +#include "CbmKF.h" #include "CbmKFVertex.h" #include "CbmL1PFFitter.h" - #include "CbmMCDataManager.h" -#include "L1Algo/L1Algo.h" -#include "L1Algo/L1Branch.h" -#include "L1Algo/L1Field.h" -#include "L1Algo/L1StsHit.h" -#include "L1AlgoInputData.h" - - -#include "CbmKF.h" +#include "CbmMuchGeoScheme.h" +#include "CbmMuchModuleGem.h" +#include "CbmMuchPad.h" +#include "CbmMuchStation.h" +#include "CbmMvdDetector.h" +#include "CbmMvdStationPar.h" #include "CbmStsFindTracks.h" #include "CbmStsParSetModule.h" #include "CbmStsParSetSensor.h" #include "CbmStsParSetSensorCond.h" #include "CbmStsSetup.h" #include "CbmStsStation.h" -#include "FairEventHeader.h" -#include "FairRunAna.h" - +#include "CbmTofCell.h" +#include "CbmTofDigiBdfPar.h" +#include "CbmTofDigiPar.h" // in tof/TofParam #include "CbmTrdParModDigi.h" // for CbmTrdModule #include "CbmTrdParSetDigi.h" // for CbmTrdParSetDigi -#include "CbmTofCell.h" -#include "CbmTofDigiPar.h" // in tof/TofParam +#include "FairEventHeader.h" +#include "FairRunAna.h" #include "TGeoArb8.h" #include "TGeoBoolNode.h" #include "TGeoCompositeShape.h" #include "TGeoManager.h" - -#include "CbmMuchGeoScheme.h" -#include "CbmMuchModuleGem.h" -#include "CbmMuchPad.h" -#include "CbmMuchStation.h" -#include "TGeoManager.h" #include "TGeoNode.h" -#include <list> - -#include "L1Event.h" #include "TMatrixD.h" #include "TROOT.h" #include "TRandom3.h" @@ -63,14 +52,19 @@ #include "TVectorD.h" #include <TFile.h> -#include "CbmMvdDetector.h" -#include "CbmMvdStationPar.h" - #include <fstream> #include <iomanip> #include <iostream> +#include <list> #include <vector> +#include "L1Algo/L1Algo.h" +#include "L1Algo/L1Branch.h" +#include "L1Algo/L1Field.h" +#include "L1Algo/L1StsHit.h" +#include "L1AlgoInputData.h" +#include "L1Event.h" + using std::cout; using std::endl; using std::ios; @@ -167,11 +161,14 @@ CbmL1::CbmL1() , fTofHitDigiMatches(nullptr) , fTofHits(nullptr) , fDigiPar(nullptr) + , fTofDigiBdfPar(nullptr) , fPerfFile(nullptr) , fHistoDir(nullptr) , vStsHits() + , SortedIndex(0) + , StsIndex(0) , vMCTracks() , vHitMCRef() , dFEI2vMCPoints() @@ -186,15 +183,12 @@ CbmL1::CbmL1() , fExtrapolateToTheEndOfSTS(false) , fTimesliceMode(0) , fTopoPerformance(nullptr) - , fEventEfficiency() { + , fEventEfficiency() +{ if (!fInstance) fInstance = this; } -CbmL1::CbmL1(const char* name, - Int_t iVerbose, - Int_t _fPerformance, - int fSTAPDataMode_, - TString fSTAPDataDir_, +CbmL1::CbmL1(const char* name, Int_t iVerbose, Int_t _fPerformance, int fSTAPDataMode_, TString fSTAPDataDir_, int findParticleMode_) : FairTask(name, iVerbose) , algo(0) @@ -220,9 +214,7 @@ CbmL1::CbmL1(const char* name, fPerformance(_fPerformance) , fSTAPDataMode(fSTAPDataMode_) , fSTAPDataDir(fSTAPDataDir_) - , - - fTrackingLevel(2) + , fTrackingLevel(2) , // really doesn't used fMomentumCutOff(0.1) , // really doesn't used @@ -276,11 +268,14 @@ CbmL1::CbmL1(const char* name, , fTofHitDigiMatches(nullptr) , fTofHits(nullptr) , fDigiPar(nullptr) + , fTofDigiBdfPar(nullptr) , fPerfFile(nullptr) , fHistoDir(nullptr) , vStsHits() + , SortedIndex(0) + , StsIndex(0) , vMCTracks() , vHitMCRef() , dFEI2vMCPoints() @@ -295,16 +290,19 @@ CbmL1::CbmL1(const char* name, , fExtrapolateToTheEndOfSTS(false) , fTimesliceMode(0) , fTopoPerformance(nullptr) - , fEventEfficiency() { + , fEventEfficiency() +{ if (!fInstance) fInstance = this; } -CbmL1::~CbmL1() { +CbmL1::~CbmL1() +{ if (fInstance == this) fInstance = 0; } -void CbmL1::SetParContainers() { +void CbmL1::SetParContainers() +{ FairRunAna* ana = FairRunAna::Instance(); FairRuntimeDb* rtdb = ana->GetRuntimeDb(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); @@ -314,21 +312,23 @@ void CbmL1::SetParContainers() { // Trigger loading of STS parameters // If L1 runs standalone the parameters are not required by any STS task // but they are needed by CbmStsSetup - fStsParSetModule = - dynamic_cast<CbmStsParSetModule*>(rtdb->getContainer("CbmStsParSetModule")); - fStsParSetSensor = - dynamic_cast<CbmStsParSetSensor*>(rtdb->getContainer("CbmStsParSetSensor")); - fStsParSetSensorCond = dynamic_cast<CbmStsParSetSensorCond*>( - rtdb->getContainer("CbmStsParSetSensorCond")); + fStsParSetModule = dynamic_cast<CbmStsParSetModule*>(rtdb->getContainer("CbmStsParSetModule")); + fStsParSetSensor = dynamic_cast<CbmStsParSetSensor*>(rtdb->getContainer("CbmStsParSetSensor")); + fStsParSetSensorCond = dynamic_cast<CbmStsParSetSensorCond*>(rtdb->getContainer("CbmStsParSetSensorCond")); + + fTofDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); + rtdb->initContainers(ana->GetRunId()); // needed to get tracking stations for ToF hits } -InitStatus CbmL1::ReInit() { +InitStatus CbmL1::ReInit() +{ SetParContainers(); return Init(); } -InitStatus CbmL1::Init() { +InitStatus CbmL1::Init() +{ fData = new L1AlgoInputData(); if (fVerbose > 1) { @@ -339,39 +339,19 @@ InitStatus CbmL1::Init() { Y[0] = y[0] = W[0] = o[0] = 0x1B; // escape character cout << endl << endl; - cout << " " << W - << " " - << o; - cout << " " << W - << " ===////====================================================== " - << o; - cout << " " << W - << " = = " - << o; - cout << " " << W << " = " << Y << "L1 on-line finder" - << W << " = " << o; - cout << " " << W - << " = = " - << o; - cout << " " << W << " = " << W << "Cellular Automaton 3.1 Vector" << y - << " with " << W << "KF Quadro" << y << " technology" << W << " = " - << o; - cout << " " << W - << " = = " - << o; - cout << " " << W << " = " << y << "Designed for CBM collaboration" << W - << " = " << o; - cout << " " << W << " = " << y << "All rights reserved" << W - << " = " << o; - cout << " " << W - << " = = " - << o; - cout << " " << W - << " ========================================================////= " - << o; - cout << " " << W - << " " - << o; + cout << " " << W << " " << o; + cout << " " << W << " ===////====================================================== " << o; + cout << " " << W << " = = " << o; + cout << " " << W << " = " << Y << "L1 on-line finder" << W << " = " << o; + cout << " " << W << " = = " << o; + cout << " " << W << " = " << W << "Cellular Automaton 3.1 Vector" << y << " with " << W << "KF Quadro" << y + << " technology" << W << " = " << o; + cout << " " << W << " = = " << o; + cout << " " << W << " = " << y << "Designed for CBM collaboration" << W << " = " << o; + cout << " " << W << " = " << y << "All rights reserved" << W << " = " << o; + cout << " " << W << " = = " << o; + cout << " " << W << " ========================================================////= " << o; + cout << " " << W << " " << o; cout << endl << endl; } @@ -379,9 +359,8 @@ InitStatus CbmL1::Init() { FairRunAna* Run = FairRunAna::Instance(); { - fUseMVD = 1; - CbmStsFindTracks* FindTask = - L1_DYNAMIC_CAST<CbmStsFindTracks*>(Run->GetTask("STSFindTracks")); + fUseMVD = 1; + CbmStsFindTracks* FindTask = L1_DYNAMIC_CAST<CbmStsFindTracks*>(Run->GetTask("STSFindTracks")); if (FindTask) fUseMVD = FindTask->MvdUsage(); } @@ -425,8 +404,7 @@ InitStatus CbmL1::Init() { LOG(info) << GetName() << ": running in time-slice mode."; fTimeSlice = NULL; fTimeSlice = (CbmTimeSlice*) fManger->GetObject("TimeSlice."); - if (fTimeSlice == NULL) - LOG(fatal) << GetName() << ": No time slice branch in tree!"; + if (fTimeSlice == NULL) LOG(fatal) << GetName() << ": No time slice branch in tree!"; } //? time-slice mode @@ -434,12 +412,9 @@ InitStatus CbmL1::Init() { LOG(info) << GetName() << ": running in event mode."; - listStsClusters = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsCluster")); - listStsHitMatch = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHitMatch")); - listStsClusterMatch = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsClusterMatch")); + listStsClusters = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsCluster")); + listStsHitMatch = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHitMatch")); + listStsClusterMatch = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsClusterMatch")); listStsHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("StsHit")); @@ -452,8 +427,8 @@ InitStatus CbmL1::Init() { fMuchPoints = 0; listMuchHitMatches = 0; - - } else { + } + else { fMuchPixelHits = (TClonesArray*) fManger->GetObject("MuchPixelHit"); @@ -464,8 +439,8 @@ InitStatus CbmL1::Init() { fTrdHitMatches = 0; fTrdPoints = 0; listTrdHits = 0; - - } else { + } + else { listTrdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("TrdHit")); } @@ -474,16 +449,15 @@ InitStatus CbmL1::Init() { fTofPoints = 0; fTofHitDigiMatches = 0; fTofHits = 0; - - } else { + } + else { fTofHits = (TClonesArray*) fManger->GetObject("TofHit"); } if (fPerformance) { - CbmMCDataManager* mcManager = - (CbmMCDataManager*) fManger->GetObject("MCDataManager"); + CbmMCDataManager* mcManager = (CbmMCDataManager*) fManger->GetObject("MCDataManager"); if (NULL == mcManager) LOG(fatal) << GetName() << ": No CbmMCDataManager!"; fStsPoints = mcManager->InitBranch("StsPoint"); @@ -496,32 +470,28 @@ InitStatus CbmL1::Init() { if (fTimesliceMode) { fEventList = (CbmMCEventList*) fManger->GetObject("MCEventList."); - if (NULL == fEventList) - LOG(fatal) << GetName() << ": No MCEventList data!"; + if (NULL == fEventList) LOG(fatal) << GetName() << ": No MCEventList data!"; } if (!fUseMVD) { listMvdPts = 0; listMvdHitMatches = 0; - } else { - listMvdPts = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdPoint")); - listMvdDigiMatches = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdDigiMatch")); - listMvdHitMatches = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHitMatch")); + } + else { + listMvdPts = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdPoint")); + listMvdDigiMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdDigiMatch")); + listMvdHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHitMatch")); if (!listMvdHitMatches && listMvdPts) - LOG(error) - << "No listMvdHitMatches provided, performance is not done correctly"; + LOG(error) << "No listMvdHitMatches provided, performance is not done correctly"; } if (!fUseTRD) { fTrdPoints = 0; fTrdHitMatches = 0; fTrdPoints = 0; - - } else { + } + else { fTrdHitMatches = (TClonesArray*) fManger->GetObject("TrdHitMatch"); fTrdPoints = mcManager->InitBranch("TrdPoint"); } @@ -529,28 +499,27 @@ InitStatus CbmL1::Init() { if (!fUseMUCH) { fMuchPoints = 0; listMuchHitMatches = 0; + } + else { - } else { - - fDigisMuch = (TClonesArray*) fManger->GetObject("MuchDigi"); - fDigiMatchesMuch = (TClonesArray*) fManger->GetObject("MuchDigiMatch"); - fClustersMuch = (TClonesArray*) fManger->GetObject("MuchCluster"); - fMuchPoints = mcManager->InitBranch("MuchPoint"); - listMuchHitMatches = - L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MuchPixelHitMatch")); + fDigisMuch = (TClonesArray*) fManger->GetObject("MuchDigi"); + fDigiMatchesMuch = (TClonesArray*) fManger->GetObject("MuchDigiMatch"); + fClustersMuch = (TClonesArray*) fManger->GetObject("MuchCluster"); + fMuchPoints = mcManager->InitBranch("MuchPoint"); + listMuchHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MuchPixelHitMatch")); } if (!fUseTOF) { fTofPoints = 0; fTofHitDigiMatches = 0; - } else { - - fTofPoints = mcManager->InitBranch("TofPoint"); - fTofHitDigiMatches = - static_cast<TClonesArray*>(fManger->GetObject("TofCalDigiMatch")); } + else { - } else { + fTofPoints = mcManager->InitBranch("TofPoint"); + fTofHitDigiMatches = static_cast<TClonesArray*>(fManger->GetObject("TofHitMatch")); + } + } + else { listMvdPts = 0; listMvdHitMatches = 0; fTrdPoints = 0; @@ -561,9 +530,8 @@ InitStatus CbmL1::Init() { fTofPoints = 0; fTofHitDigiMatches = 0; } - if (!fUseMVD) { - listMvdHits = 0; - } else { + if (!fUseMVD) { listMvdHits = 0; } + else { listMvdHits = L1_DYNAMIC_CAST<TClonesArray*>(fManger->GetObject("MvdHit")); } @@ -574,7 +542,6 @@ InitStatus CbmL1::Init() { NTOFStation = 0; NStation = 0; - if (fUseTOF) NTOFStation = 1; algo = &algo_static; @@ -584,8 +551,7 @@ InitStatus CbmL1::Init() { for (int i = 0; i < 3; i++) { Double_t point[3] = {0, 0, 2.5 * i}; Double_t B[3] = {0, 0, 0}; - if (CbmKF::Instance()->GetMagneticField()) - CbmKF::Instance()->GetMagneticField()->GetFieldValue(point, B); + if (CbmKF::Instance()->GetMagneticField()) CbmKF::Instance()->GetMagneticField()->GetFieldValue(point, B); geo.push_back(2.5 * i); geo.push_back(B[0]); geo.push_back(B[1]); @@ -610,8 +576,7 @@ InitStatus CbmL1::Init() { if (fUseTRD) { Int_t layerCounter = 0; TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes(); - for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); - iTopNode++) { + for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) { TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode)); if (TString(topNode->GetName()).Contains("trd")) { TObjArray* layers = topNode->GetNodes(); @@ -626,60 +591,48 @@ InitStatus CbmL1::Init() { // count ToF parameters - float z_average = 0; - - float z_min = 100000; - float z_max = 0; + int maxTofStation = 0; + int St = 0; - float r_max = 0; - float r_min = 100000; - - if (fUseTOF) { + vector<int> NHits; + vector<float> Z_pos; - Int_t nrOfCells = fDigiPar->GetNrOfModules(); + NHits.resize(10, 0); + Z_pos.resize(10, 0); - for (Int_t icell = 0; icell < nrOfCells; ++icell) { - Int_t cellId = fDigiPar->GetCellId(icell); - CbmTofCell* fCellInfo = fDigiPar->GetCell(cellId); + if (fUseTOF) { - Double_t x = fCellInfo->GetX(); - Double_t y = fCellInfo->GetY(); - Double_t z = fCellInfo->GetZ(); - // Double_t dx = fCellInfo->GetSizex(); - // Double_t dy = fCellInfo->GetSizey(); + if (fTofHits) { + for (int j = 0; j < fTofHits->GetEntries(); j++) { - if (z < z_min) z_min = z; - if (z > z_max) z_max = z; + CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j)); + St = fTofDigiBdfPar->GetNbTrackingStations(); - if (x < r_min) r_min = x; - if (x > r_max) r_max = x; + if (maxTofStation < St) maxTofStation = St; - if (y < r_min) r_min = y; - if (y > r_max) r_max = y; + TVector3 pos; + mh->Position(pos); + Z_pos[St] += pos.Z(); + NHits[St]++; + } + } + if (fUseTOF) NTOFStation = fTofDigiBdfPar->GetNbTrackingStations(); - z_average += z; - } - z_average = z_average / nrOfCells; + for (int i = 0; i < (maxTofStation + 1); i++) + Z_pos[i] = Z_pos[i] / NHits[i]; } CbmStsSetup* stsSetup = CbmStsSetup::Instance(); if (!stsSetup->IsInit()) { stsSetup->Init(nullptr); } - if (!stsSetup->IsModuleParsInit()) { - stsSetup->SetModuleParameters(fStsParSetModule); - } - if (!stsSetup->IsSensorParsInit()) { - stsSetup->SetSensorParameters(fStsParSetSensor); - } - if (!stsSetup->IsSensorCondInit()) { - stsSetup->SetSensorConditions(fStsParSetSensorCond); - } + if (!stsSetup->IsModuleParsInit()) { stsSetup->SetModuleParameters(fStsParSetModule); } + if (!stsSetup->IsSensorParsInit()) { stsSetup->SetSensorParameters(fStsParSetSensor); } + if (!stsSetup->IsSensorCondInit()) { stsSetup->SetSensorConditions(fStsParSetSensorCond); } NMvdStations = (fUseMVD) ? CbmKF::Instance()->vMvdMaterial.size() : 0; NStsStations = CbmStsSetup::Instance()->GetNofStations(); - NStation = - NMvdStations + NStsStations + NMuchStations + NTrdStations + NTOFStation; + NStation = NMvdStations + NStsStations + NMuchStations + NTrdStations + NTOFStation; geo.push_back(NStation); geo.push_back(NMvdStations); geo.push_back(NStsStations); @@ -760,8 +713,7 @@ InitStatus CbmL1::Init() { geo.push_back(t.R); geo.push_back(t.RadLength); - fscal f_phi = 0, f_sigma = mvdStationPar->GetXRes(ist) / 10000, - b_phi = 3.14159265358 / 2., + fscal f_phi = 0, f_sigma = mvdStationPar->GetXRes(ist) / 10000, b_phi = 3.14159265358 / 2., b_sigma = mvdStationPar->GetYRes(ist) / 10000; geo.push_back(f_phi); geo.push_back(f_sigma); @@ -773,15 +725,12 @@ InitStatus CbmL1::Init() { if (ist >= NMvdStations && ist < (NMvdStations + NStsStations)) { - CbmStsStation* station = - CbmStsSetup::Instance()->GetStation(ist - NMvdStations); + CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ist - NMvdStations); geo.push_back(0); geo.push_back(station->GetZ()); geo.push_back(station->GetSensorD()); geo.push_back(0); - geo.push_back(station->GetYmax() < station->GetXmax() - ? station->GetXmax() - : station->GetYmax()); + geo.push_back(station->GetYmax() < station->GetXmax() ? station->GetXmax() : station->GetYmax()); geo.push_back(station->GetRadLength()); double Pi = 3.14159265358; @@ -806,17 +755,14 @@ InitStatus CbmL1::Init() { Ymax = station->GetYmax(); } - if ((ist < (NMvdStations + NStsStations + NMuchStations)) - && (ist >= (NMvdStations + NStsStations))) { + if ((ist < (NMvdStations + NStsStations + NMuchStations)) && (ist >= (NMvdStations + NStsStations))) { int iStation = (ist - NMvdStations - NStsStations) / 3; - CbmMuchStation* station = - (CbmMuchStation*) fGeoScheme->GetStation(iStation); + CbmMuchStation* station = (CbmMuchStation*) fGeoScheme->GetStation(iStation); - CbmMuchLayer* layer = - station->GetLayer((ist - NMvdStations - NStsStations) % 3); + CbmMuchLayer* layer = station->GetLayer((ist - NMvdStations - NStsStations) % 3); // CbmMuchModuleGem* module = (CbmMuchModuleGem*) CbmMuchGeoScheme::Instance()->GetModule(0,0,0,0); @@ -856,8 +802,7 @@ InitStatus CbmL1::Init() { int ModuleId = fTrdDigiPar->GetModuleId(num); - CbmTrdParModDigi* module = - (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(ModuleId); + CbmTrdParModDigi* module = (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(ModuleId); // if (!true_station[ist]) continue; @@ -870,8 +815,7 @@ InitStatus CbmL1::Init() { geo.push_back(2 * module->GetSizeX()); geo.push_back(10); - fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., - b_sigma = 1 / 10000; + fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., b_sigma = 1 / 10000; geo.push_back(f_phi); geo.push_back(f_sigma); geo.push_back(b_phi); @@ -879,20 +823,24 @@ InitStatus CbmL1::Init() { Xmax = Ymax = 20; } - if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations - + NTOFStation)) - && (ist - >= (NMvdStations + NStsStations + NMuchStations + NTrdStations))) { + if ((ist < (NMvdStations + NStsStations + NTrdStations + NMuchStations + NTOFStation)) + && (ist >= (NMvdStations + NStsStations + NMuchStations + NTrdStations))) { geo.push_back(4); - geo.push_back(z_average); - geo.push_back(z_max - z_min); + + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 0)) geo.push_back(245); + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 1)) geo.push_back(249); + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 2)) geo.push_back(261); + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 3)) geo.push_back(266); + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 4)) geo.push_back(278); + if (ist == (NMvdStations + NStsStations + NTrdStations + NMuchStations + 5)) geo.push_back(283); + + geo.push_back(10); /// TODO: add Tof width dz geo.push_back(0); - geo.push_back(r_max); + geo.push_back(150); /// TODO: add Tof max radius geo.push_back(10); - fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., - b_sigma = 1 / 10000; + fscal f_phi = 0, f_sigma = 1 / 10000, b_phi = 3.14159265358 / 2., b_sigma = 1 / 10000; geo.push_back(f_phi); geo.push_back(f_sigma); geo.push_back(b_phi); @@ -981,8 +929,7 @@ InitStatus CbmL1::Init() { int ind2, ind = geo.size(); ReadSTAPGeoData(geo, ind2); if (ind2 != ind) - LOG(error) << "-E- CbmL1: Read geometry from file " - << fSTAPDataDir + "geo_algo.txt was NOT successful."; + LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful."; } algo->Init(geo, fUseHitErrors, fmCBMmode); @@ -1003,31 +950,25 @@ InitStatus CbmL1::Init() { stationNameMvd += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameMvd); if (!hStaRadLen) { - cout << "L1: incorrect " << fMvdMatBudgetFileName << " file. No " - << stationNameMvd << "\n"; + cout << "L1: incorrect " << fMvdMatBudgetFileName << " file. No " << stationNameMvd << "\n"; exit(1); } - const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y - const float RMax = - hStaRadLen->GetXaxis()->GetXmax(); // should be same as min + const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y + const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min algo->fRadThick[iSta].SetBins(NBins, RMax); algo->fRadThick[iSta].table.resize(NBins); for (int iB = 0; iB < NBins; iB++) { algo->fRadThick[iSta].table[iB].resize(NBins); for (int iB2 = 0; iB2 < NBins; iB2++) { - algo->fRadThick[iSta].table[iB][iB2] = - 0.01 * hStaRadLen->GetBinContent(iB, iB2); + algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2); // Correction for holes in material map - if (algo->fRadThick[iSta].table[iB][iB2] - < algo->vStations[iSta].materialInfo.RadThick[0]) + if (algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0]) if (iB2 > 0 && iB2 < NBins - 1) - algo->fRadThick[iSta].table[iB][iB2] = - TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), - 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); + algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), + 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); // Correction for the incorrect harcoded value of RadThick of MVD stations - if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) - algo->fRadThick[iSta].table[iB][iB2] = 0.0015; + if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = 0.0015; // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1035,7 +976,8 @@ InitStatus CbmL1::Init() { rlFile->Close(); rlFile->Delete(); gFile = oldfile; - } else { + } + else { LOG(warn) << "No MVD material budget file is found. Homogenious budget " "will be used"; for (int iSta = 0; iSta < algo->NMvdStations; iSta++) { @@ -1044,8 +986,7 @@ InitStatus CbmL1::Init() { 100); // mvd should be in +-100 cm square algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); - algo->fRadThick[iSta].table[0][0] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0]; } } } @@ -1053,32 +994,25 @@ InitStatus CbmL1::Init() { TFile* oldfile = gFile; TFile* rlFile = new TFile(fStsMatBudgetFileName); cout << "STS Material budget file is " << fStsMatBudgetFileName << ".\n"; - for (int j = 1, iSta = algo->NMvdStations; - iSta < (algo->NMvdStations + NStsStations); - iSta++, j++) { + for (int j = 1, iSta = algo->NMvdStations; iSta < (algo->NMvdStations + NStsStations); iSta++, j++) { TString stationNameSts = stationName; stationNameSts += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts); if (!hStaRadLen) { - cout << "L1: incorrect " << fStsMatBudgetFileName << " file. No " - << stationNameSts << "\n"; + cout << "L1: incorrect " << fStsMatBudgetFileName << " file. No " << stationNameSts << "\n"; exit(1); } - const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y - const float RMax = - hStaRadLen->GetXaxis()->GetXmax(); // should be same as min + const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y + const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min algo->fRadThick[iSta].SetBins(NBins, RMax); algo->fRadThick[iSta].table.resize(NBins); for (int iB = 0; iB < NBins; iB++) { algo->fRadThick[iSta].table[iB].resize(NBins); for (int iB2 = 0; iB2 < NBins; iB2++) { - algo->fRadThick[iSta].table[iB][iB2] = - 0.01 * hStaRadLen->GetBinContent(iB, iB2); - if (algo->fRadThick[iSta].table[iB][iB2] - < algo->vStations[iSta].materialInfo.RadThick[0]) - algo->fRadThick[iSta].table[iB][iB2] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2); + if (algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0]) + algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0]; // cout << " iSta " << iSta << " iB " << iB << "iB2 " << iB2 << " RadThick " << algo->fRadThick[iSta].table[iB][iB2] << endl; } } @@ -1086,17 +1020,15 @@ InitStatus CbmL1::Init() { rlFile->Close(); rlFile->Delete(); gFile = oldfile; - } else { + } + else { LOG(warn) << "No STS material budget file is found. Homogenious budget " "will be used"; - for (int iSta = algo->NMvdStations; - iSta < (algo->NMvdStations + NStsStations); - iSta++) { + for (int iSta = algo->NMvdStations; iSta < (algo->NMvdStations + NStsStations); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); - algo->fRadThick[iSta].table[0][0] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1104,24 +1036,20 @@ InitStatus CbmL1::Init() { if (fMuchMatBudgetFileName != "") { TFile* oldfile = gFile; TFile* rlFile = new TFile(fMuchMatBudgetFileName); - cout << "Much Material budget file is " << fMuchMatBudgetFileName - << ".\n"; - for (int j = 1, iSta = (NStsStations + NMvdStations); - iSta < (NStsStations + NMvdStations + NMuchStations); + cout << "Much Material budget file is " << fMuchMatBudgetFileName << ".\n"; + for (int j = 1, iSta = (NStsStations + NMvdStations); iSta < (NStsStations + NMvdStations + NMuchStations); iSta++, j++) { TString stationNameSts = stationName; stationNameSts += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts); if (!hStaRadLen) { - cout << "L1: incorrect " << fMuchMatBudgetFileName << " file. No " - << stationNameSts << "\n"; + cout << "L1: incorrect " << fMuchMatBudgetFileName << " file. No " << stationNameSts << "\n"; exit(1); } - const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y - const float RMax = - hStaRadLen->GetXaxis()->GetXmax(); // should be same as min + const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y + const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min algo->fRadThick[iSta].SetBins(NBins, RMax); algo->fRadThick[iSta].table.resize(NBins); @@ -1129,21 +1057,17 @@ InitStatus CbmL1::Init() { algo->fRadThick[iSta].table[iB].resize(NBins); float hole = 0.15; for (int iB2 = 0; iB2 < NBins; iB2++) { - algo->fRadThick[iSta].table[iB][iB2] = - 0.01 * hStaRadLen->GetBinContent(iB, iB2); + algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2); // Correction for holes in material map //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0]) if (iB2 > 0 && iB2 < NBins - 1) - algo->fRadThick[iSta].table[iB][iB2] = - TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), - 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); + algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), + 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); // Correction for the incorrect harcoded value of RadThick of MVD stations - if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) - hole = algo->fRadThick[iSta].table[iB][iB2]; - if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) - algo->fRadThick[iSta].table[iB][iB2] = hole; + if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2]; + if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole; // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1151,17 +1075,15 @@ InitStatus CbmL1::Init() { rlFile->Close(); rlFile->Delete(); gFile = oldfile; - } else { + } + else { LOG(warn) << "No Much material budget file is found. Homogenious budget " "will be used"; - for (int iSta = (NStsStations + NMvdStations); - iSta < (NStsStations + NMvdStations + NMuchStations); - iSta++) { + for (int iSta = (NStsStations + NMvdStations); iSta < (NStsStations + NMvdStations + NMuchStations); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); - algo->fRadThick[iSta].table[0][0] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1171,21 +1093,18 @@ InitStatus CbmL1::Init() { TFile* rlFile = new TFile(fTrdMatBudgetFileName); cout << "TRD Material budget file is " << fTrdMatBudgetFileName << ".\n"; for (int j = 1, iSta = (NStsStations + NMvdStations + NMuchStations); - iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations); - iSta++, j++) { + iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations); iSta++, j++) { TString stationNameSts = stationName; stationNameSts += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts); if (!hStaRadLen) { - cout << "L1: incorrect " << fTrdMatBudgetFileName << " file. No " - << stationNameSts << "\n"; + cout << "L1: incorrect " << fTrdMatBudgetFileName << " file. No " << stationNameSts << "\n"; exit(1); } - const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y - const float RMax = - hStaRadLen->GetXaxis()->GetXmax(); // should be same as min + const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y + const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min algo->fRadThick[iSta].SetBins(NBins, RMax); algo->fRadThick[iSta].table.resize(NBins); @@ -1193,21 +1112,17 @@ InitStatus CbmL1::Init() { algo->fRadThick[iSta].table[iB].resize(NBins); float hole = 0.15; for (int iB2 = 0; iB2 < NBins; iB2++) { - algo->fRadThick[iSta].table[iB][iB2] = - 0.01 * hStaRadLen->GetBinContent(iB, iB2); + algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2); // Correction for holes in material map //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0]) if (iB2 > 0 && iB2 < NBins - 1) - algo->fRadThick[iSta].table[iB][iB2] = - TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), - 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); + algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), + 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); // Correction for the incorrect harcoded value of RadThick of MVD stations - if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) - hole = algo->fRadThick[iSta].table[iB][iB2]; - if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) - algo->fRadThick[iSta].table[iB][iB2] = hole; + if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2]; + if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole; // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1215,17 +1130,16 @@ InitStatus CbmL1::Init() { rlFile->Close(); rlFile->Delete(); gFile = oldfile; - } else { + } + else { LOG(warn) << "No TRD material budget file is found. Homogenious budget " "will be used"; for (int iSta = (NStsStations + NMvdStations + NMuchStations); - iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations); - iSta++) { + iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); - algo->fRadThick[iSta].table[0][0] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1234,25 +1148,19 @@ InitStatus CbmL1::Init() { TFile* oldfile = gFile; TFile* rlFile = new TFile(fTofMatBudgetFileName); cout << "TOF Material budget file is " << fTofMatBudgetFileName << ".\n"; - for (int j = 1, - iSta = - (NStsStations + NMvdStations + NMuchStations + NTrdStations); - iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations - + NTOFStation); - iSta++, j++) { + for (int j = 1, iSta = (NStsStations + NMvdStations + NMuchStations + NTrdStations); + iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations + NTOFStation); iSta++, j++) { TString stationNameSts = stationName; stationNameSts += j; TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(stationNameSts); if (!hStaRadLen) { - cout << "L1: incorrect " << fTofMatBudgetFileName << " file. No " - << stationNameSts << "\n"; + cout << "L1: incorrect " << fTofMatBudgetFileName << " file. No " << stationNameSts << "\n"; exit(1); } - const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y - const float RMax = - hStaRadLen->GetXaxis()->GetXmax(); // should be same as min + const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y + const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min algo->fRadThick[iSta].SetBins(NBins, RMax); algo->fRadThick[iSta].table.resize(NBins); @@ -1260,21 +1168,17 @@ InitStatus CbmL1::Init() { algo->fRadThick[iSta].table[iB].resize(NBins); float hole = 0.0015; for (int iB2 = 0; iB2 < NBins; iB2++) { - algo->fRadThick[iSta].table[iB][iB2] = - 0.01 * hStaRadLen->GetBinContent(iB, iB2); + algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2); // Correction for holes in material map //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0]) if (iB2 > 0 && iB2 < NBins - 1) - algo->fRadThick[iSta].table[iB][iB2] = - TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), - 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); + algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1), + 0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1)); // Correction for the incorrect harcoded value of RadThick of MVD stations - if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) - hole = algo->fRadThick[iSta].table[iB][iB2]; - if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) - algo->fRadThick[iSta].table[iB][iB2] = hole; + if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2]; + if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole; // algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0]; } } @@ -1282,19 +1186,16 @@ InitStatus CbmL1::Init() { rlFile->Close(); rlFile->Delete(); gFile = oldfile; - } else { + } + else { LOG(warn) << "No TOF material budget file is found. Homogenious budget " "will be used"; - for (int iSta = - (NStsStations + NMvdStations + NMuchStations + NTrdStations); - iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations - + NTOFStation); - iSta++) { + for (int iSta = (NStsStations + NMvdStations + NMuchStations + NTrdStations); + iSta < (NStsStations + NMvdStations + NMuchStations + NTrdStations + NTOFStation); iSta++) { algo->fRadThick[iSta].SetBins(1, 100); algo->fRadThick[iSta].table.resize(1); algo->fRadThick[iSta].table[0].resize(1); - algo->fRadThick[iSta].table[0][0] = - algo->vStations[iSta].materialInfo.RadThick[0]; + algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0]; } } return kSUCCESS; @@ -1303,10 +1204,41 @@ InitStatus CbmL1::Init() { void CbmL1::Exec(Option_t* /*option*/) {} -void CbmL1::Reconstruct(CbmEvent* event) { +void CbmL1::Reconstruct(CbmEvent* event) +{ static int nevent = 0; vFileEvent.clear(); + std::vector<std::pair<double, int>> SortStsHits; + + float start_t = 10000000000; + + bool newTS = 1; // whole TS processed? + float TsStart = 0; // starting time of sub-TS + float TsLength = 10000; // length of sub-TS + float TsOverlap = 15; // min length of overlap region + int FstHitinTs = 0; // 1st hit index in TS + + /// sort input hits by time + for (Int_t j = 0; j < listStsHits->GetEntries(); j++) { + CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(j)); + double t = sh->GetTime(); + if (t < start_t) start_t = t; + SortStsHits.push_back(std::pair<double, int>(t, j)); + } + + TsStart = start_t; ///reco TS start time is set to smallest hit time + + + std::sort(SortStsHits.begin(), SortStsHits.end()); + StsIndex.resize(0); + + for (int i = 0; i < SortStsHits.size(); i++) { + int j = SortStsHits[i].second; + StsIndex.push_back(j); + }; + + if (fTimesliceMode && fPerformance) { int nofEvents = fEventList->GetNofEvents(); @@ -1316,122 +1248,333 @@ void CbmL1::Reconstruct(CbmEvent* event) { int eventId = fEventList->GetEventIdByIndex(iE); vFileEvent.insert(DFSET::value_type(fileId, eventId)); } - - } else { + } + else { Int_t iFile = FairRunAna::Instance()->GetEventHeader()->GetInputFileId(); Int_t iEvent = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber(); vFileEvent.insert(DFSET::value_type(iFile, iEvent)); } - if (fVerbose > 1) { - cout << endl << "CbmL1::Exec event " << ++nevent << " ..." << endl << endl; - } + if (fVerbose > 1) { cout << endl << "CbmL1::Exec event " << ++nevent << " ..." << endl << endl; } #ifdef _OPENMP omp_set_num_threads(1); #endif // repack data - fData->Clear(); + vector<CbmL1Track> vRTracksCur; // reconstructed tracks - if (fSTAPDataMode >= 2) { // 2,3 - fData->ReadHitsFromFile(fSTAPDataDir.Data(), 1, fVerbose); - - algo->SetData(fData->GetStsHits(), - fData->GetNStsStrips(), - fData->GetStsZPos(), - fData->GetSFlag(), - fData->GetStsHitsStartIndex(), - fData->GetStsHitsStopIndex()); - } else - ReadEvent(fData, event); - - if (0) { // correct hits on MC // dbg - TRandom3 random; - vector<int> strips, zP; - strips.clear(); - zP.clear(); - for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) { - L1StsHit& h = const_cast<L1StsHit&>((*algo->vStsHits)[iH]); + + while (newTS) { + + fData->Clear(); + + if (!fTimesliceMode) { + + fData->vSFlag.clear(); + + newTS = 0; + TsStart = 0; + TsLength = 2000000000; + TsOverlap = 0; + FstHitinTs = 0; + } + + if (fSTAPDataMode >= 2) { // 2,3 + fData->ReadHitsFromFile(fSTAPDataDir.Data(), 1, fVerbose); + + algo->SetData(fData->GetStsHits(), fData->GetNStsStrips(), fData->GetStsZPos(), fData->GetSFlag(), + fData->GetStsHitsStartIndex(), fData->GetStsHitsStopIndex(), listStsHits->GetEntries()); + } + else + ReadEvent(fData, TsStart, TsLength, TsOverlap, FstHitinTs, newTS, event); + + if (0) { // correct hits on MC // dbg + TRandom3 random; + vector<int> strips, zP; + strips.clear(); + zP.clear(); + for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) { + L1StsHit& h = const_cast<L1StsHit&>((*algo->vStsHits)[iH]); #ifdef USE_EVENT_NUMBER - h.n = -1; + h.n = -1; #endif - if (vStsHits[iH].mcPointIds.size() == 0) continue; + if (vStsHits[iH].mcPointIds.size() == 0) continue; - const CbmL1MCPoint& mcp = vMCPoints[vStsHits[iH].mcPointIds[0]]; + const CbmL1MCPoint& mcp = vMCPoints[vStsHits[iH].mcPointIds[0]]; #ifdef USE_EVENT_NUMBER - h.n = mcp.event; + h.n = mcp.event; #endif - const int ista = (*algo->vSFlag)[h.f] / 4; - const L1Station& sta = algo->vStations[ista]; - if (std::find(strips.begin(), strips.end(), h.f) - != strips.end()) { // separate strips + const int ista = (*algo->vSFlag)[h.f] / 4; + const L1Station& sta = algo->vStations[ista]; + if (std::find(strips.begin(), strips.end(), h.f) != strips.end()) { // separate strips - const_cast<vector<unsigned char>*>(algo->vSFlag) - ->push_back((*algo->vSFlag)[h.f]); + const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back((*algo->vSFlag)[h.f]); - h.f = algo->NStsStrips; - algo->NStsStrips++; - } - strips.push_back(h.f); - if (std::find(strips.begin(), strips.end(), h.b) != strips.end()) { - const_cast<vector<unsigned char>*>(algo->vSFlag) - ->push_back((*algo->vSFlag)[h.b]); - h.b = algo->NStsStrips; - algo->NStsStrips++; - } - strips.push_back(h.b); - if (std::find(zP.begin(), zP.end(), h.iz) - != zP.end()) { // TODO why do we need it??gives prob=0 - h.iz = (*algo->vStsZPos).size(); - const_cast<std::vector<float>*>(algo->vStsZPos)->push_back(0); - } - zP.push_back(h.iz); + h.f = algo->NStsStrips; + algo->NStsStrips++; + } + strips.push_back(h.f); + if (std::find(strips.begin(), strips.end(), h.b) != strips.end()) { + const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back((*algo->vSFlag)[h.b]); + h.b = algo->NStsStrips; + algo->NStsStrips++; + } + strips.push_back(h.b); + if (std::find(zP.begin(), zP.end(), h.iz) != zP.end()) { // TODO why do we need it??gives prob=0 + h.iz = (*algo->vStsZPos).size(); + const_cast<std::vector<float>*>(algo->vStsZPos)->push_back(0); + } + zP.push_back(h.iz); - double u = - mcp.x * sta.frontInfo.cos_phi[0] + mcp.y * sta.frontInfo.sin_phi[0]; - double v = - mcp.x * sta.backInfo.cos_phi[0] + mcp.y * sta.backInfo.sin_phi[0]; + double u = mcp.x * sta.frontInfo.cos_phi[0] + mcp.y * sta.frontInfo.sin_phi[0]; + double v = mcp.x * sta.backInfo.cos_phi[0] + mcp.y * sta.backInfo.sin_phi[0]; #if 1 // GAUSS - u += random.Gaus(0, sqrt(sta.frontInfo.sigma2)[0]); - v += random.Gaus(0, sqrt(sta.backInfo.sigma2)[0]); + u += random.Gaus(0, sqrt(sta.frontInfo.sigma2)[0]); + v += random.Gaus(0, sqrt(sta.backInfo.sigma2)[0]); #else // UNIFORM - u += 3.5 * sqrt(sta.frontInfo.sigma2)[0] * random.Uniform(-1, 1); - v += 3.5 * sqrt(sta.backInfo.sigma2)[0] * random.Uniform(-1, 1); + u += 3.5 * sqrt(sta.frontInfo.sigma2)[0] * random.Uniform(-1, 1); + v += 3.5 * sqrt(sta.backInfo.sigma2)[0] * random.Uniform(-1, 1); #endif - h.u = u; - h.v = v; - const_cast<float&>((*algo->vStsZPos)[h.iz]) = mcp.z; + h.u = u; + h.v = v; + const_cast<float&>((*algo->vStsZPos)[h.iz]) = mcp.z; + } } + + + if (fPerformance) { + HitMatch(); + // calculate the max number of Hits\mcPoints on continuous(consecutive) stations + + for (vector<CbmL1MCTrack>::iterator it = vMCTracks.begin(); it != vMCTracks.end(); ++it) + it->Init(); + } + + if (fSTAPDataMode % 2 == 1) { // 1,3 + WriteSTAPAlgoData(); + WriteSTAPPerfData(); + }; + if (fSTAPDataMode >= 2) { // 2,3 + //ReadSTAPAlgoData(); + + ReadSTAPPerfData(); + }; + + if ((fPerformance) && (fSTAPDataMode < 2)) { InputPerformance(); } + + // FieldApproxCheck(); + // FieldIntegralCheck(); + + for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) { +#ifdef USE_EVENT_NUMBER + L1StsHit& h = const_cast<L1StsHit&>((*algo->vStsHits)[iH]); + h.n = -1; +#endif + if (vStsHits[iH].mcPointIds.size() == 0) continue; +#ifdef USE_EVENT_NUMBER + const CbmL1MCPoint& mcp = vMCPoints[vStsHits[iH].mcPointIds[0]]; + h.n = mcp.event; +#endif + } + + for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) { + CbmL1MCTrack& MC = *i; + + if (!MC.IsReconstructable()) continue; + if (!(MC.ID >= 0)) continue; + + if (MC.StsHits.size() < 4) continue; + vector<int> hitIndices(algo->NStations, -1); + + for (unsigned int iH = 0; iH < MC.StsHits.size(); iH++) { + const int hitI = MC.StsHits[iH]; + CbmL1StsHit& hit = const_cast<CbmL1StsHit&>(vStsHits[hitI]); + + hit.event = MC.iEvent; + + // const int iStation = vMCPoints[hit.mcPointIds[0]].iStation; + // hitIndices[iStation] = hitI; + } + } + + + if (fVerbose > 1) { cout << "L1 Track finder..." << endl; } + algo->CATrackFinder(); + // IdealTrackFinder(); + + + if (fVerbose > 1) { cout << "L1 Track finder ok" << endl; } + // algo->L1KFTrackFitter( fExtrapolateToTheEndOfSTS ); + + + if (fmCBMmode || fGlobalMode) { + + L1FieldValue fB0 = algo->GetVtxFieldValue(); + + if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001) && (fabs(fB0.z[0]) < 0.0000001)) + algo->KFTrackFitter_simple(); + + else + algo->L1KFTrackFitterMuch(); + } + else { + + L1FieldValue fB0 = algo->GetVtxFieldValue(); + + if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001) && (fabs(fB0.z[0]) < 0.0000001)) + algo->KFTrackFitter_simple(); + + else + algo->L1KFTrackFitter(); + } + + + if (fVerbose > 1) { cout << "L1 Track fitter ok" << endl; } + + // save recontstructed tracks + if (!fTimesliceMode) vRTracksCur.clear(); + int start_hit = 0; + + float TsStart_new = TsStart + TsLength - TsOverlap; + + for (vector<L1Track>::iterator it = algo->vTracks.begin(); it != (algo->vTracks.begin() + algo->NTracksIsecAll); + it++) { + + CbmL1Track t; + for (int i = 0; i < 7; i++) + t.T[i] = it->TFirst[i]; + for (int i = 0; i < 21; i++) + t.C[i] = it->CFirst[i]; + for (int i = 0; i < 7; i++) + t.TLast[i] = it->TLast[i]; + for (int i = 0; i < 21; i++) + t.CLast[i] = it->CLast[i]; + for (int k = 0; k < 7; k++) + t.Tpv[k] = it->Tpv[k]; + for (int k = 0; k < 21; k++) + t.Cpv[k] = it->Cpv[k]; + t.chi2 = it->chi2; + t.NDF = it->NDF; + //t.T[4] = it->Momentum; + t.StsHits.clear(); + t.fTrackTime = it->fTrackTime; + + vector<int> StsHitsLocal; + + for (int i = 0; i < it->NHits; i++) { + int start_hit1 = start_hit; + if (algo->vRecoHits[start_hit1] > vStsHits.size() - 1) start_hit++; + else if (fTimesliceMode) + t.StsHits.push_back(((*algo->vStsHits)[algo->vRecoHits[start_hit]]).ID); + else + t.StsHits.push_back(algo->vRecoHits[start_hit]); + + StsHitsLocal.push_back(algo->vRecoHits[start_hit++]); + } + + t.mass = 0.1395679; // pion mass + t.is_electron = 0; + + t.SetId(vRTracksCur.size()); + // CbmL1Track* prtra = &t; + + int indd = 0; + bool isInOverlap = 0; + + + for (int i = 0; i < StsHitsLocal.size(); i++) { + // if ((*ih) > int(vStsHits.size() - 1)) { + // indd = 1; + // break; + // } + + if (vStsHits[StsHitsLocal[i]].t >= (TsStart + TsLength - TsOverlap)) { + isInOverlap = 1; + if (TsStart_new > vStsHits[StsHitsLocal[i]].t) TsStart_new = vStsHits[StsHitsLocal[i]].t; + } + + int nMCPoints = vStsHits[StsHitsLocal[i]].mcPointIds.size(); + for (int iP = 0; iP < nMCPoints; iP++) { + int iMC = vStsHits[StsHitsLocal[i]].mcPointIds[iP]; + if (iMC > int(vMCPoints.size() - 1)) { + // cout << " iMC " << iMC << " vMCPoints.size() " << vMCPoints.size() << endl; + indd = 1; + } + } + } + + if (indd == 1) continue; + + if ((fTimesliceMode) && (isInOverlap == 1)) { + + continue; ///Discard tracks from overlap region + + /// set strips as unused + for (int i = 0; i < StsHitsLocal.size(); i++) { + algo->SetFUnUsed(const_cast<unsigned char&>((*algo->vSFlag)[vStsHits[StsHitsLocal[i]].f])); + algo->SetFUnUsed(const_cast<unsigned char&>((*algo->vSFlag)[vStsHits[StsHitsLocal[i]].b])); + } + } + vRTracksCur.push_back(t); + } + + + for (int i = 0; i < listStsHits->GetEntries(); i++) { + + CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(StsIndex[i])); + float time = sh->GetTime(); + + if (TsStart_new <= time) { + FstHitinTs = i; + break; + } + }; + + if (fTimesliceMode) TsStart = TsStart_new; ///Set new TS strat to earliest discarted track + + if (fTimesliceMode) cout << "CA Track Finder: " << algo->CATime << " s/sub-ts" << endl << endl; } if (fPerformance) { + + float start = 0; + float end = 10000000000.f; + int fHit = 0; + bool stop = 0; + + ReadEvent(fData, start, end, start, fHit, stop, event); HitMatch(); // calculate the max number of Hits\mcPoints on continuous(consecutive) stations - for (vector<CbmL1MCTrack>::iterator it = vMCTracks.begin(); - it != vMCTracks.end(); - ++it) + for (vector<CbmL1MCTrack>::iterator it = vMCTracks.begin(); it != vMCTracks.end(); ++it) it->Init(); } + // + // if (fSTAPDataMode % 2 == 1) { // 1,3 + // WriteSTAPAlgoData(); + // WriteSTAPPerfData(); + // }; + // if (fSTAPDataMode >= 2) { // 2,3 + // //ReadSTAPAlgoData(); + // + // ReadSTAPPerfData(); + // }; - if (fSTAPDataMode % 2 == 1) { // 1,3 - WriteSTAPAlgoData(); - WriteSTAPPerfData(); - }; - if (fSTAPDataMode >= 2) { // 2,3 - //ReadSTAPAlgoData(); + for (int iTrack = 0; iTrack < vRTracksCur.size(); iTrack++) { + + for (int iHit = 0; iHit < vRTracksCur[iTrack].StsHits.size(); iHit++) + if (fTimesliceMode) vRTracksCur[iTrack].StsHits[iHit] = SortedIndex[vRTracksCur[iTrack].StsHits[iHit]]; + + vRTracks.push_back(vRTracksCur[iTrack]); + } - ReadSTAPPerfData(); - }; if ((fPerformance) && (fSTAPDataMode < 2)) { InputPerformance(); } - // FieldApproxCheck(); - // FieldIntegralCheck(); for (unsigned int iH = 0; iH < (*algo->vStsHits).size(); ++iH) { #ifdef USE_EVENT_NUMBER @@ -1445,9 +1588,7 @@ void CbmL1::Reconstruct(CbmEvent* event) { #endif } - for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); - i != vMCTracks.end(); - ++i) { + for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) { CbmL1MCTrack& MC = *i; if (!MC.IsReconstructable()) continue; @@ -1461,111 +1602,10 @@ void CbmL1::Reconstruct(CbmEvent* event) { CbmL1StsHit& hit = const_cast<CbmL1StsHit&>(vStsHits[hitI]); hit.event = MC.iEvent; - - // const int iStation = vMCPoints[hit.mcPointIds[0]].iStation; - // hitIndices[iStation] = hitI; } } - if (fVerbose > 1) { cout << "L1 Track finder..." << endl; } - algo->CATrackFinder(); - // IdealTrackFinder(); - - - if (fVerbose > 1) { cout << "L1 Track finder ok" << endl; } - // algo->L1KFTrackFitter( fExtrapolateToTheEndOfSTS ); - - - if (fmCBMmode || fGlobalMode) { - - L1FieldValue fB0 = algo->GetVtxFieldValue(); - - if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001) - && (fabs(fB0.z[0]) < 0.0000001)) - algo->KFTrackFitter_simple(); - - else - algo->L1KFTrackFitterMuch(); - } else { - - L1FieldValue fB0 = algo->GetVtxFieldValue(); - - if ((fabs(fB0.x[0]) < 0.0000001) && (fabs(fB0.y[0]) < 0.0000001) - && (fabs(fB0.z[0]) < 0.0000001)) - algo->KFTrackFitter_simple(); - - else - algo->L1KFTrackFitter(); - } - - - if (fVerbose > 1) { cout << "L1 Track fitter ok" << endl; } - - // save recontstructed tracks - vRTracks.clear(); - int start_hit = 0; - - for (vector<L1Track>::iterator it = algo->vTracks.begin(); - it != (algo->vTracks.begin() + algo->NTracksIsecAll); - it++) { - - CbmL1Track t; - for (int i = 0; i < 7; i++) - t.T[i] = it->TFirst[i]; - for (int i = 0; i < 21; i++) - t.C[i] = it->CFirst[i]; - for (int i = 0; i < 7; i++) - t.TLast[i] = it->TLast[i]; - for (int i = 0; i < 21; i++) - t.CLast[i] = it->CLast[i]; - for (int k = 0; k < 7; k++) - t.Tpv[k] = it->Tpv[k]; - for (int k = 0; k < 21; k++) - t.Cpv[k] = it->Cpv[k]; - t.chi2 = it->chi2; - t.NDF = it->NDF; - //t.T[4] = it->Momentum; - t.StsHits.clear(); - t.fTrackTime = it->fTrackTime; - - for (int i = 0; i < it->NHits; i++) { - int start_hit1 = start_hit; - - if (algo->vRecoHits[start_hit1] > vStsHits.size() - 1) - start_hit++; - else - t.StsHits.push_back(algo->vRecoHits[start_hit++]); - } - t.mass = 0.1395679; // pion mass - t.is_electron = 0; - - t.SetId(vRTracks.size()); - CbmL1Track* prtra = &t; - - int indd = 0; - - for (vector<int>::iterator ih = (prtra->StsHits).begin(); - ih != (prtra->StsHits).end(); - ++ih) { - if ((*ih) > int(vStsHits.size() - 1)) { - indd = 1; - break; - } - int nMCPoints = vStsHits[*ih].mcPointIds.size(); - for (int iP = 0; iP < nMCPoints; iP++) { - int iMC = vStsHits[*ih].mcPointIds[iP]; - if (iMC > int(vMCPoints.size() - 1)) { - // cout << " iMC " << iMC << " vMCPoints.size() " << vMCPoints.size() << endl; - indd = 1; - } - } - } - - if (indd == 1) continue; - - vRTracks.push_back(t); - } // output performance if (fPerformance) { if (fVerbose > 1) { cout << "Performance..." << endl; } @@ -1593,10 +1633,13 @@ void CbmL1::Reconstruct(CbmEvent* event) { if ((symbol == 'e') || (symbol == 'q')) exit(0); } while (symbol != '\n'); } + + // } } // ----- Finish CbmStsFitPerformanceTask task ----------------------------- -void CbmL1::Finish() { +void CbmL1::Finish() +{ TDirectory* curr = gDirectory; TFile* currentFile = gFile; @@ -1612,9 +1655,9 @@ void CbmL1::Finish() { } -void CbmL1::writedir2current(TObject* obj) { - if (!obj->IsFolder()) - obj->Write(); +void CbmL1::writedir2current(TObject* obj) +{ + if (!obj->IsFolder()) obj->Write(); else { TDirectory* cur = gDirectory; TDirectory* sub = cur->mkdir(obj->GetName()); @@ -1629,13 +1672,12 @@ void CbmL1::writedir2current(TObject* obj) { /// ----- Ideal Tracking ----------------------------- -void CbmL1::IdealTrackFinder() { +void CbmL1::IdealTrackFinder() +{ algo->vTracks.clear(); algo->vRecoHits.clear(); - for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); - i != vMCTracks.end(); - ++i) { + for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) { CbmL1MCTrack& MC = *i; if (!MC.IsReconstructable()) continue; @@ -1651,7 +1693,8 @@ void CbmL1::IdealTrackFinder() { for (unsigned int iH = 0; iH < MC.StsHits.size(); iH++) { const int hitI = MC.StsHits[iH]; const CbmL1StsHit& hit = vStsHits[hitI]; - const int iStation = vMCPoints[hit.mcPointIds[0]].iStation; + + const int iStation = vMCPoints[hit.mcPointIds[0]].iStation; if (iStation >= 0) hitIndices[iStation] = hitI; } @@ -1691,7 +1734,8 @@ void CbmL1::IdealTrackFinder() { /// ----- STandAlone Package service-functions ----------------------------- -void CbmL1::WriteSTAPGeoData(const vector<float>& geo_) { +void CbmL1::WriteSTAPGeoData(const vector<float>& geo_) +{ // write geo in file TString fgeo_name = fSTAPDataDir + "geo_algo.txt"; ofstream fgeo(fgeo_name); @@ -1716,8 +1760,7 @@ void CbmL1::WriteSTAPAlgoData() // must be called after ReadEvent // if ( vNEvent <= maxNEvent ) { if (1) { - if (vNEvent == 1) - fadata.open(fadata_name, fstream::out); // begin new file + if (vNEvent == 1) fadata.open(fadata_name, fstream::out); // begin new file else fadata.open(fadata_name, fstream::out | fstream::app); @@ -1780,14 +1823,12 @@ void CbmL1::WriteSTAPAlgoData() // must be called after ReadEvent // write StsHitsStartIndex and StsHitsStopIndex n = 20; for (int i = 0; i < n; i++) { - if (int(algo->MaxNStations) + 1 > i) - fadata << algo->StsHitsStartIndex[i] << endl; + if (int(algo->MaxNStations) + 1 > i) fadata << algo->StsHitsStartIndex[i] << endl; else fadata << 0 << endl; }; for (int i = 0; i < n; i++) { - if (int(algo->MaxNStations) + 1 > i) - fadata << algo->StsHitsStopIndex[i] << endl; + if (int(algo->MaxNStations) + 1 > i) fadata << algo->StsHitsStopIndex[i] << endl; else fadata << 0 << endl; }; @@ -1795,8 +1836,8 @@ void CbmL1::WriteSTAPAlgoData() // must be called after ReadEvent fadata.close(); } - cout << "-I- CbmL1: CATrackFinder data for event number " << vNEvent - << " have been written in file " << fadata_name << endl; + cout << "-I- CbmL1: CATrackFinder data for event number " << vNEvent << " have been written in file " << fadata_name + << endl; vNEvent++; } // void CbmL1::WriteSTAPAlgoData() @@ -1813,8 +1854,7 @@ void CbmL1::WriteSTAPPerfData() // must be called after ReadEvent // if ( vNEvent <= maxNEvent ) { if (1) { - if (vNEvent == 1) - fpdata.open(fpdata_name, fstream::out); // begin new file + if (vNEvent == 1) fpdata.open(fpdata_name, fstream::out); // begin new file else fpdata.open(fpdata_name, fstream::out | fstream::app); @@ -1946,8 +1986,8 @@ void CbmL1::WriteSTAPPerfData() // must be called after ReadEvent } fpdata.close(); } - cout << "-I- CbmL1: Data for performance of event number " << vNEvent - << " have been written in file " << fpdata_name << endl; + cout << "-I- CbmL1: Data for performance of event number " << vNEvent << " have been written in file " << fpdata_name + << endl; vNEvent++; } // void CbmL1::WriteSTAPPerfData() @@ -1965,7 +2005,8 @@ istream& CbmL1::eatwhite(istream& is) // skip spaces //void CbmL1::ReadSTAPGeoData(vector<float> geo_, int &size) //void CbmL1::ReadSTAPGeoData(vector<fscal> geo_, int &size) -void CbmL1::ReadSTAPGeoData(vector<fscal>& geo_, int& size) { +void CbmL1::ReadSTAPGeoData(vector<fscal>& geo_, int& size) +{ TString fgeo_name = fSTAPDataDir + "geo_algo.txt"; ifstream fgeo(fgeo_name); @@ -1981,7 +2022,8 @@ void CbmL1::ReadSTAPGeoData(vector<fscal>& geo_, int& size) { fgeo.close(); } // void CbmL1::ReadSTAPGeoData(void* geo_, int &size) -void CbmL1::ReadSTAPAlgoData() { +void CbmL1::ReadSTAPAlgoData() +{ static int nEvent = 1; static fstream fadata; TString fadata_name = fSTAPDataDir + "data_algo.txt"; @@ -1989,11 +2031,9 @@ void CbmL1::ReadSTAPAlgoData() { if (1) { if (nEvent == 1) fadata.open(fadata_name, fstream::in); - if (algo->vStsHits) - const_cast<std::vector<L1StsHit>*>(algo->vStsHits)->clear(); + if (algo->vStsHits) const_cast<std::vector<L1StsHit>*>(algo->vStsHits)->clear(); algo->NStsStrips = 0; - if (algo->vStsZPos) - const_cast<std::vector<float>*>(algo->vStsZPos)->clear(); + if (algo->vStsZPos) const_cast<std::vector<float>*>(algo->vStsZPos)->clear(); if (algo->vSFlag) const_cast<vector<unsigned char>*>(algo->vSFlag)->clear(); // check correct position in file @@ -2001,9 +2041,7 @@ void CbmL1::ReadSTAPAlgoData() { int nEv; fadata >> s; fadata >> nEv; - if (nEv != nEvent) - cout << "-E- CbmL1: Can't read event number " << nEvent << " from file " - << fadata_name << endl; + if (nEv != nEvent) cout << "-E- CbmL1: Can't read event number " << nEvent << " from file " << fadata_name << endl; int n; // number of elements // read algo->vStsStrips @@ -2030,8 +2068,7 @@ void CbmL1::ReadSTAPAlgoData() { for (int i = 0; i < n; i++) { int element; fadata >> element; - const_cast<vector<unsigned char>*>(algo->vSFlag) - ->push_back(static_cast<unsigned char>(element)); + const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back(static_cast<unsigned char>(element)); } if (fVerbose >= 4) { cout << "vSFlag[" << n << "]" @@ -2045,8 +2082,7 @@ void CbmL1::ReadSTAPAlgoData() { int element_iz; for (int i = 0; i < n; i++) { L1StsHit element; - fadata >> element_f >> element_b >> element_n >> element_iz >> element.u - >> element.v >> element.t_reco; + fadata >> element_f >> element_b >> element_n >> element_iz >> element.u >> element.v >> element.t_reco; element.f = static_cast<THitI>(element_f); element.b = static_cast<THitI>(element_b); element.iz = static_cast<TZPosI>(element_iz); @@ -2061,24 +2097,22 @@ void CbmL1::ReadSTAPAlgoData() { for (int i = 0; i < n; i++) { int tmp; fadata >> tmp; - if (int(algo->MaxNStations) + 1 > i) - (const_cast<unsigned int&>(algo->StsHitsStartIndex[i]) = tmp); + if (int(algo->MaxNStations) + 1 > i) (const_cast<unsigned int&>(algo->StsHitsStartIndex[i]) = tmp); } for (int i = 0; i < n; i++) { int tmp; fadata >> tmp; - if (int(algo->MaxNStations) + 1 > i) - (const_cast<unsigned int&>(algo->StsHitsStopIndex[i]) = tmp); + if (int(algo->MaxNStations) + 1 > i) (const_cast<unsigned int&>(algo->StsHitsStopIndex[i]) = tmp); } - cout << "-I- CbmL1: CATrackFinder data for event " << nEvent - << " has been read from file " << fadata_name << " successfully." - << endl; + cout << "-I- CbmL1: CATrackFinder data for event " << nEvent << " has been read from file " << fadata_name + << " successfully." << endl; } nEvent++; } // void CbmL1::ReadSTAPAlgoData() -void CbmL1::ReadSTAPPerfData() { +void CbmL1::ReadSTAPPerfData() +{ static int nEvent = 1; static fstream fpdata; TString fpdata_name = fSTAPDataDir + "data_perfo.txt"; @@ -2100,8 +2134,7 @@ void CbmL1::ReadSTAPPerfData() { fpdata >> nEv; if (nEv != nEvent) - cout << "-E- CbmL1: Performance: can't read event number " << nEvent - << " from file " + cout << "-E- CbmL1: Performance: can't read event number " << nEvent << " from file " << "data_perfo.txt" << endl; // vMCPoints int n; // number of elements @@ -2246,15 +2279,15 @@ void CbmL1::ReadSTAPPerfData() { // fpdata.close(); // cout << " -I- Performance: data read from file " << "data_perfo.txt" << " successfully"<< endl; // } - cout << "-I- CbmL1: L1Performance data for event " << nEvent - << " has been read from file " << fpdata_name << " successfully." - << endl; + cout << "-I- CbmL1: L1Performance data for event " << nEvent << " has been read from file " << fpdata_name + << " successfully." << endl; } // if (nEvent <= maxNEvent) nEvent++; } // void CbmL1::ReadSTAPPerfData() -void CbmL1::WriteSIMDKFData() { +void CbmL1::WriteSIMDKFData() +{ static bool first = 1; /// Write geometry info @@ -2273,23 +2306,19 @@ void CbmL1::WriteSIMDKFData() { rfg[1] = 0.; rfg[2] = 0.; dMF->GetFieldValue(rfg, bfg); - FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " - << endl; + FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " << endl; rfg[0] = 0.; rfg[1] = 0.; rfg[2] = 2.5; dMF->GetFieldValue(rfg, bfg); - FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " - << endl; + FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " << endl; rfg[0] = 0.; rfg[1] = 0.; rfg[2] = 5.0; dMF->GetFieldValue(rfg, bfg); - FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " - << endl - << endl; + FileGeo << rfg[2] << " " << bfg[0] << " " << bfg[1] << " " << bfg[2] << " " << endl << endl; FileGeo << NStation << endl; const int M = 5; // polinom order @@ -2309,12 +2338,12 @@ void CbmL1::WriteSIMDKFData() { b_sigma = 5.e-4; z = t.z; Xmax = Ymax = t.R; - } else { - CbmStsStation* station = - CbmStsSetup::Instance()->GetStation(ist - NMvdStations); - f_phi = station->GetSensorRotation(); - b_phi = f_phi; - double Pi = 3.14159265358; + } + else { + CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ist - NMvdStations); + f_phi = station->GetSensorRotation(); + b_phi = f_phi; + double Pi = 3.14159265358; f_phi += station->GetSensorStereoAngle(0) * Pi / 180.; b_phi += station->GetSensorStereoAngle(1) * Pi / 180.; f_sigma = station->GetSensorPitch(0) / TMath::Sqrt(12); @@ -2398,9 +2427,9 @@ void CbmL1::WriteSIMDKFData() { FileGeo << t.z << " "; FileGeo << t.dz << " "; FileGeo << t.RadLength << " "; - } else if (ist < (NStsStations + NMvdStations)) { - CbmStsStation* station = - CbmStsSetup::Instance()->GetStation(ist - NMvdStations); + } + else if (ist < (NStsStations + NMvdStations)) { + CbmStsStation* station = CbmStsSetup::Instance()->GetStation(ist - NMvdStations); FileGeo << station->GetZ() << " "; FileGeo << station->GetSensorD() << " "; FileGeo << station->GetRadLength() << " "; @@ -2436,16 +2465,15 @@ void CbmL1::WriteSIMDKFData() { McTracksIn.open("mctracksin.dat", fstream::out); McTracksOut.open("mctracksout.dat", fstream::out); first = 0; - } else { + } + else { Tracks.open("tracks.dat", fstream::out | fstream::app); McTracksCentr.open("mctrackscentr.dat", fstream::out | fstream::app); McTracksIn.open("mctracksin.dat", fstream::out | fstream::app); McTracksOut.open("mctracksout.dat", fstream::out | fstream::app); } - for (vector<CbmL1Track>::iterator RecTrack = vRTracks.begin(); - RecTrack != vRTracks.end(); - ++RecTrack) { + for (vector<CbmL1Track>::iterator RecTrack = vRTracks.begin(); RecTrack != vRTracks.end(); ++RecTrack) { if (RecTrack->IsGhost()) continue; CbmL1MCTrack* MCTrack = RecTrack->GetMCTrack(); @@ -2464,7 +2492,8 @@ void CbmL1::WriteSIMDKFData() { y[jHit] = MvdH->GetY(); z[jHit] = MvdH->GetZ(); jHit++; - } else { + } + else { CbmStsHit* StsH = (CbmStsHit*) listStsHits->At(h.ExtIndex); x[jHit] = StsH->GetX(); y[jHit] = StsH->GetY(); @@ -2473,40 +2502,32 @@ void CbmL1::WriteSIMDKFData() { } } - Tracks << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " - << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " " - << MCTrack->pz << " " << MCTrack->q << " " << NHits << endl; + Tracks << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " " + << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NHits << endl; float AngleXAxis = 0, AngleYAxis = 0; for (int i = 0; i < NHits; i++) - Tracks << " " << st[i] << " " << x[i] << " " << y[i] << " " << z[i] - << " " << AngleXAxis << " " << AngleYAxis << endl; + Tracks << " " << st[i] << " " << x[i] << " " << y[i] << " " << z[i] << " " << AngleXAxis << " " << AngleYAxis + << endl; Tracks << endl; int NMCPoints = (MCTrack->Points).size(); - McTracksIn << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " - << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " " - << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl; - McTracksOut << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " - << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py << " " - << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl; - McTracksCentr << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " - << MCTrack->z << " " << MCTrack->px << " " << MCTrack->py - << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints - << endl; + McTracksIn << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " " + << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl; + McTracksOut << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px << " " + << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl; + McTracksCentr << TrNumber << " " << MCTrack->x << " " << MCTrack->y << " " << MCTrack->z << " " << MCTrack->px + << " " << MCTrack->py << " " << MCTrack->pz << " " << MCTrack->q << " " << NMCPoints << endl; for (int iPoint = 0; iPoint < NMCPoints; iPoint++) { CbmL1MCPoint& MCPoint = vMCPoints[MCTrack->Points[iPoint]]; - McTracksIn << " " << MCPoint.iStation << " " << MCPoint.xIn << " " - << MCPoint.yIn << " " << MCPoint.zIn << " " << MCPoint.pxIn - << " " << MCPoint.pyIn << " " << MCPoint.pzIn << endl; - McTracksOut << " " << MCPoint.iStation << " " << MCPoint.xOut << " " - << MCPoint.yOut << " " << MCPoint.zOut << " " << MCPoint.pxOut - << " " << MCPoint.pyOut << " " << MCPoint.pzOut << endl; - McTracksCentr << " " << MCPoint.iStation << " " << MCPoint.x << " " - << MCPoint.y << " " << MCPoint.z << " " << MCPoint.px << " " - << MCPoint.py << " " << MCPoint.pz << endl; + McTracksIn << " " << MCPoint.iStation << " " << MCPoint.xIn << " " << MCPoint.yIn << " " << MCPoint.zIn << " " + << MCPoint.pxIn << " " << MCPoint.pyIn << " " << MCPoint.pzIn << endl; + McTracksOut << " " << MCPoint.iStation << " " << MCPoint.xOut << " " << MCPoint.yOut << " " << MCPoint.zOut + << " " << MCPoint.pxOut << " " << MCPoint.pyOut << " " << MCPoint.pzOut << endl; + McTracksCentr << " " << MCPoint.iStation << " " << MCPoint.x << " " << MCPoint.y << " " << MCPoint.z << " " + << MCPoint.px << " " << MCPoint.py << " " << MCPoint.pz << endl; } McTracksIn << endl; McTracksOut << endl; diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index bfd59ff358684b79da178be750991ee99babbf2a..c94341998e95c6c6e95c93f1c6005d9cb463ebb8 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -26,24 +26,21 @@ #include "CbmL1Vtx.h" //#include "L1Algo/L1Algo.h" +#include "CbmEvent.h" #include "CbmL1MCPoint.h" #include "CbmL1MCTrack.h" #include "CbmL1StsHit.h" - - -#include "FairDetector.h" -#include "FairRootManager.h" -#include "FairTask.h" - -#include "CbmEvent.h" #include "CbmMCTrack.h" +#include "CbmMvdHit.h" +#include "CbmMvdPoint.h" #include "CbmStsCluster.h" #include "CbmStsDigi.h" #include "CbmStsHit.h" #include "CbmStsPoint.h" -#include "CbmMvdHit.h" -#include "CbmMvdPoint.h" +#include "FairDetector.h" +#include "FairRootManager.h" +#include "FairTask.h" //#include "CbmMCEventHeader.h" //#include "L1AlgoInputData.h" @@ -90,6 +87,7 @@ class KFTopoPerformance; class CbmStsParSetSensor; class CbmStsParSetSensorCond; class CbmStsParSetModule; +class CbmTofDigiBdfPar; class CbmL1HitStore { public: @@ -140,30 +138,16 @@ public: * @param fSTAPDataMode_ - way to work with files for the standalone package: 0 - off , 1 - write, 2 - read data and work only with it, 3 - write and read (debug) * @param findParticleMode_ : 0 - don't run FindParticles; 1 - run, all MC particle is reco-able; 2 - run, MC particle is reco-able if created from reco-able tracks; 3 - run, MC particle is reco-able if created from reconstructed tracks */ - CbmL1(const char* name, - Int_t iVerbose = 1, - Int_t _fPerformance = 0, - int fSTAPDataMode_ = 0, - TString fSTAPDataDir_ = "./", - int findParticleMode_ = 0); + CbmL1(const char* name, Int_t iVerbose = 1, Int_t _fPerformance = 0, int fSTAPDataMode_ = 0, + TString fSTAPDataDir_ = "./", int findParticleMode_ = 0); ~CbmL1(/*if (targetFieldSlice) delete;*/); - void SetStsMaterialBudgetFileName(TString fileName) { - fStsMatBudgetFileName = fileName; - } - void SetMvdMaterialBudgetFileName(TString fileName) { - fMvdMatBudgetFileName = fileName; - } - void SetMuchMaterialBudgetFileName(TString fileName) { - fMuchMatBudgetFileName = fileName; - } - void SetTrdMaterialBudgetFileName(TString fileName) { - fTrdMatBudgetFileName = fileName; - } - void SetTofMaterialBudgetFileName(TString fileName) { - fTofMatBudgetFileName = fileName; - } + void SetStsMaterialBudgetFileName(TString fileName) { fStsMatBudgetFileName = fileName; } + void SetMvdMaterialBudgetFileName(TString fileName) { fMvdMatBudgetFileName = fileName; } + void SetMuchMaterialBudgetFileName(TString fileName) { fMuchMatBudgetFileName = fileName; } + void SetTrdMaterialBudgetFileName(TString fileName) { fTrdMatBudgetFileName = fileName; } + void SetTofMaterialBudgetFileName(TString fileName) { fTofMatBudgetFileName = fileName; } void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; } void SetDataMode(int TimesliceMode) { fTimesliceMode = TimesliceMode; } void SetMuchPar(TString fileName) { fDigiFile = fileName; } @@ -197,33 +181,30 @@ private: vector<CbmL1MCPoint> vMCPoints; int nMvdPoints; vector<int> vMCPoints_in_Time_Slice; - void - IdealTrackFinder(); // just copy all reconstructable MCTracks into RecoTracks. + void IdealTrackFinder(); // just copy all reconstructable MCTracks into RecoTracks. /// Read information about hits, mcPoints and mcTracks into L1 classes - void ReadEvent(L1AlgoInputData*, CbmEvent* event = NULL); + + void ReadEvent(L1AlgoInputData* fData_, float& fTsStart, float& fTsLength, float& fTsOverlap, int& fFstHitinTs, + bool& fnewTS, CbmEvent* event = NULL); + bool ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int MVD); // help procedure bool ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int MVD); // static bool compareZ(const int &a, const int &b ); // bool compareZ(const int &a, const int &b ); void Fill_vMCTracks(); /// Input Performance - void HitMatch(); // Procedure for match hits and MCPoints. - void - FieldApproxCheck(); // Build histos with difference between Field map and approximated field - void - FieldIntegralCheck(); // Build 2D histo: dependence of the field integral on phi and theta - void - InputPerformance(); // Build histos about input data, like hit pulls, etc. + void HitMatch(); // Procedure for match hits and MCPoints. + void FieldApproxCheck(); // Build histos with difference between Field map and approximated field + void FieldIntegralCheck(); // Build 2D histo: dependence of the field integral on phi and theta + void InputPerformance(); // Build histos about input data, like hit pulls, etc. void TimeHist(); /// Reconstruction Performance - void - TrackMatch(); // Procedure for match Reconstructed and MC Tracks. Should be called before Performances + void TrackMatch(); // Procedure for match Reconstructed and MC Tracks. Should be called before Performances void EfficienciesPerformance(); // calculate efficiencies - void - TrackFitPerformance(); // pulls & residuals. Can be called only after Performance() - void HistoPerformance(); // fill some histograms and calculate efficiencies + void TrackFitPerformance(); // pulls & residuals. Can be called only after Performance() + void HistoPerformance(); // fill some histograms and calculate efficiencies /// STandAlone Package service-functions void WriteSTAPGeoData(const vector<float>& geo); // create geo_algo.dat @@ -242,9 +223,8 @@ private: static void writedir2current(TObject* obj); // help procedure int NStation, NMvdStations, NStsStations, NMuchStations, NTrdStations, - NTOFStation; // number of detector stations (all\sts\mvd) - Int_t - fPerformance; // 0 - w\o perf. 1 - L1-Efficiency definition. 2 - QA-Eff.definition + NTOFStation; // number of detector stations (all\sts\mvd) + Int_t fPerformance; // 0 - w\o perf. 1 - L1-Efficiency definition. 2 - QA-Eff.definition int fSTAPDataMode; // way to work with file for standalone package. 0 (off) , 1 (write), 2 (read data and work only with it), 3 (debug - write and read) TString fSTAPDataDir; @@ -313,6 +293,7 @@ private: TClonesArray* fTofHitDigiMatches; // CbmMatches array TClonesArray* fTofHits; // CbmMatches array CbmTofDigiPar* fDigiPar; + CbmTofDigiBdfPar* fTofDigiBdfPar; struct TH1FParameters { @@ -333,7 +314,9 @@ private: //CbmMCEventHeader* fEvent; /// Used data = Repacked input data vector<CbmL1StsHit> vStsHits; // hits with hit-mcpoint match information - // vector<CbmL1MCPoint> vMCPoints; + // vector<CbmL1MCPoint> vMCPoints; + vector<int> SortedIndex; + vector<int> StsIndex; vector<CbmL1MCTrack> vMCTracks; vector<int> vHitMCRef; // indices of MCPoints in vMCPoints, indexed by index of hit in algo->vStsHits array. According to StsMatch. Used for IdealResponce @@ -343,9 +326,7 @@ private: typedef std::map<Double_t, Int_t> DFEI2I; DFEI2I dFEI2vMCPoints; DFEI2I dFEI2vMCTracks; - inline Double_t dFEI(int file, int event, int idx) { - return (1000 * idx) + file + (0.0001 * event); - } + inline Double_t dFEI(int file, int event, int idx) { return (1000 * idx) + file + (0.0001 * event); } // DFEI2I::iterator map_it; L1AlgoInputData* fData; diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 462520bff617fe85d2dbf0cf35dfbc4f645c82ed..f2411dac8b18f2da562eca7cea3e4cc6c419c6c8 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -17,21 +17,21 @@ #include "CbmKF.h" #include "CbmL1.h" -#include "CbmMatch.h" -#include "CbmStsAddress.h" -#include "CbmStsSetup.h" -#include "L1Algo/L1Algo.h" -#include "L1AlgoInputData.h" - #include "CbmMatch.h" #include "CbmMuchGeoScheme.h" #include "CbmMuchPixelHit.h" #include "CbmMuchPoint.h" +#include "CbmStsAddress.h" +#include "CbmStsSetup.h" +#include "CbmTofDigiBdfPar.h" #include "CbmTofHit.h" #include "CbmTofPoint.h" #include "CbmTrdHit.h" #include "CbmTrdPoint.h" +#include "L1Algo/L1Algo.h" +#include "L1AlgoInputData.h" + //#include "CbmMvdHitMatch.h" @@ -50,7 +50,8 @@ using std::vector; //#define STSIDEALHITS -static bool compareZ(const int& a, const int& b) { +static bool compareZ(const int& a, const int& b) +{ // return (CbmL1::Instance()->vMCPoints[a].z < CbmL1::Instance()->vMCPoints[b].z); const CbmL1* l1 = CbmL1::Instance(); return l1->Get_Z_vMCPoint(a) < l1->Get_Z_vMCPoint(b); @@ -60,7 +61,7 @@ struct TmpHit { // used for sort Hits before writing in the normal arrays int iStripF, iStripB; // indices of strips int iStation; - int ExtIndex; // index of hit in the TClonesArray array ( negative for MVD ) + int ExtIndex; // index of hit in the TClonesArray array ( negative for MVD ) double u_front, u_back; // positions of strips double x, y; // position of hit double dx, dy, dxy; @@ -70,14 +71,16 @@ struct TmpHit { // used for sort Hits before writing in the normal arrays int Det; int id; int track; - static bool Compare(const TmpHit& a, const TmpHit& b) { - return (a.iStation < b.iStation) - || ((a.iStation == b.iStation) && (a.y < b.y)); + static bool Compare(const TmpHit& a, const TmpHit& b) + { + return (a.iStation < b.iStation) || ((a.iStation == b.iStation) && (a.y < b.y)); } }; /// Repack data from Clones Arrays to L1 arrays -void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { +void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, float& TsOverlap, int& FstHitinTs, + bool& newTS, CbmEvent* event) +{ if (fVerbose >= 10) cout << "ReadEvent: start." << endl; // clear arrays for next event @@ -122,9 +125,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fPerformance) { Fill_vMCTracks(); - for (DFSET::iterator set_it = vFileEvent.begin(); - set_it != vFileEvent.end(); - ++set_it) { + for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) { Int_t iFile = set_it->first; Int_t iEvent = set_it->second; @@ -155,8 +156,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { vMCPoints.push_back(MC); vMCPoints_in_Time_Slice.push_back(0); - dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC), - vMCPoints.size() - 1)); + dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC), vMCPoints.size() - 1)); nMvdPoints++; } } @@ -173,8 +173,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { MC.iStation = -1; L1Station* sta = algo->vStations + NMvdStations; for (Int_t iSt = 0; iSt < NStsStations; iSt++) - MC.iStation = - (MC.z > sta[iSt].z[0] - 2.5) ? (NMvdStations + iSt) : MC.iStation; + MC.iStation = (MC.z > sta[iSt].z[0] - 2.5) ? (NMvdStations + iSt) : MC.iStation; Double_t dtrck = dFEI(iFile, iEvent, MC.ID); DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck); @@ -186,8 +185,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { vMCPoints.push_back(MC); vMCPoints_in_Time_Slice.push_back(0); - dFEI2vMCPoints.insert(DFEI2I::value_type( - dFEI(iFile, iEvent, iMC + nMvdPoints), vMCPoints.size() - 1)); + dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints), vMCPoints.size() - 1)); nStsPoints++; } } @@ -203,9 +201,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { MC.iStation = -1; L1Station* sta = algo->vStations + NMvdStations + NStsStations; for (Int_t iSt = 0; iSt < NMuchStations; iSt++) - MC.iStation = (MC.z > sta[iSt].z[0] - 2.5) - ? (NMvdStations + NStsStations + iSt) - : MC.iStation; + MC.iStation = (MC.z > sta[iSt].z[0] - 2.5) ? (NMvdStations + NStsStations + iSt) : MC.iStation; Double_t dtrck = dFEI(iFile, iEvent, MC.ID); DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck); @@ -218,9 +214,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { vMCPoints.push_back(MC); vMCPoints_in_Time_Slice.push_back(0); - dFEI2vMCPoints.insert(DFEI2I::value_type( - dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints), - vMCPoints.size() - 1)); + dFEI2vMCPoints.insert( + DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints), vMCPoints.size() - 1)); nMuchPoints++; } } @@ -236,14 +231,11 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) { - MC.iStation = -1; - L1Station* sta = - algo->vStations + NMvdStations + NStsStations + NMuchStations; + MC.iStation = -1; + L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations; for (Int_t iSt = 0; iSt < NTrdStations; iSt++) MC.iStation = - (MC.z > sta[iSt].z[0] - 4.0) - ? (NMvdStations + NStsStations + NMuchStations + iSt) - : MC.iStation; + (MC.z > sta[iSt].z[0] - 4.0) ? (NMvdStations + NStsStations + NMuchStations + iSt) : MC.iStation; Double_t dtrck = dFEI(iFile, iEvent, MC.ID); @@ -257,9 +249,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { vMCPoints.push_back(MC); vMCPoints_in_Time_Slice.push_back(0); - dFEI2vMCPoints.insert(DFEI2I::value_type( - dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints + nMuchPoints), - vMCPoints.size() - 1)); + dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints + nMuchPoints), + vMCPoints.size() - 1)); nTrdPoints++; } } @@ -271,21 +262,19 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j)); - CbmMatch* matchHitMatch = - L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(j)); + CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(j)); if (matchHitMatch->GetNofLinks() > 0) { CbmLink* link = (CbmLink*) &matchHitMatch->GetLink(0); - CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get( - link->GetFile(), link->GetEntry(), link->GetIndex()); + CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(link->GetFile(), link->GetEntry(), link->GetIndex()); for (int iLink = 1; iLink < matchHitMatch->GetNofLinks(); iLink++) { CbmLink* link1 = (CbmLink*) &(matchHitMatch->GetLink(iLink)); - CbmTofPoint* pt_cur = (CbmTofPoint*) fTofPoints->Get( - link1->GetFile(), link1->GetEntry(), link1->GetIndex()); + CbmTofPoint* pt_cur = + (CbmTofPoint*) fTofPoints->Get(link1->GetFile(), link1->GetEntry(), link1->GetIndex()); TVector3 pos_cur, pos_old, pos_hit; @@ -294,8 +283,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { pt->Position(pos_old); mh->Position(pos_hit); - if (fabs(pos_cur.Z() - pos_hit.Z()) - < fabs(pos_old.Z() - pos_hit.Z())) { + if (fabs(pos_cur.Z() - pos_hit.Z()) < fabs(pos_old.Z() - pos_hit.Z())) { pt = pt_cur; link = link1; } @@ -315,14 +303,15 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (eventNr != iEvent) continue; - if (!ReadMCPoint(&MC, - ToFPointsMatch[iMC]->GetIndex(), - ToFPointsMatch[iMC]->GetFile(), - ToFPointsMatch[iMC]->GetEntry(), - 4)) { + if (!ReadMCPoint(&MC, ToFPointsMatch[iMC]->GetIndex(), ToFPointsMatch[iMC]->GetFile(), + ToFPointsMatch[iMC]->GetEntry(), 4)) { - MC.iStation = - NMvdStations + NStsStations + NMuchStations + NTrdStations; + MC.iStation = -1; + L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations + NTrdStations; + for (Int_t iSt = 0; iSt < NTOFStation; iSt++) + MC.iStation = (MC.z > sta[iSt].z[0] - 1.5) + ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt) + : MC.iStation; Double_t dtrck = dFEI(iFile, eventNr, MC.ID); DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck); @@ -338,10 +327,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { vMCPoints_in_Time_Slice.push_back(0); dFEI2vMCPoints.insert(DFEI2I::value_type( - dFEI(iFile, - eventNr, - ToFPointsMatch[iMC]->GetIndex() + nMvdPoints + nStsPoints - + nMuchPoints + nTrdPoints), + dFEI(iFile, eventNr, + ToFPointsMatch[iMC]->GetIndex() + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints), vMCPoints.size() - 1)); nTofPoints++; } @@ -352,16 +339,12 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { for (unsigned int iTr = 0; iTr < vMCTracks.size(); iTr++) { - sort( - vMCTracks[iTr].Points.begin(), vMCTracks[iTr].Points.end(), compareZ); + sort(vMCTracks[iTr].Points.begin(), vMCTracks[iTr].Points.end(), compareZ); if (vMCTracks[iTr].mother_ID >= 0) { - Double_t dtrck = dFEI(vMCTracks[iTr].iFile, - vMCTracks[iTr].iEvent, - vMCTracks[iTr].mother_ID); + Double_t dtrck = dFEI(vMCTracks[iTr].iFile, vMCTracks[iTr].iEvent, vMCTracks[iTr].mother_ID); DFEI2I::iterator map_it = dFEI2vMCTracks.find(dtrck); - if (map_it == dFEI2vMCTracks.end()) - vMCTracks[iTr].mother_ID = -2; + if (map_it == dFEI2vMCTracks.end()) vMCTracks[iTr].mother_ID = -2; else vMCTracks[iTr].mother_ID = map_it->second; } @@ -390,6 +373,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { { CbmMvdHit* mh = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(j)); th.ExtIndex = -(1 + j); + th.id = j; th.iStation = mh->GetStationNr(); // th.iSector = 0; int iStripF = j; @@ -412,10 +396,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.y = pos.Y(); L1Station& st = algo->vStations[th.iStation]; - th.u_front = - th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - th.u_back = - th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } th.Det = 0; th.iMC = -1; @@ -430,8 +412,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { // } if (fPerformance) { if (listMvdHitMatches) { - CbmMatch* mvdHitMatch = - L1_DYNAMIC_CAST<CbmMatch*>(listMvdHitMatches->At(j)); + CbmMatch* mvdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(listMvdHitMatches->At(j)); if (mvdHitMatch->GetNofLinks() > 0) @@ -482,33 +463,35 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fVerbose >= 10) cout << "ReadEvent: mvd hits are gotten." << endl; + Int_t nEntSts = 0; if (listStsHits) { - Int_t nEnt = 0; - if (fTimesliceMode) - nEnt = listStsHits->GetEntries(); + + if (fTimesliceMode) nEntSts = listStsHits->GetEntries(); else - nEnt = (event ? event->GetNofData(ECbmDataType::kStsHit) - : listStsHits->GetEntries()); + nEntSts = (event ? event->GetNofData(ECbmDataType::kStsHit) : listStsHits->GetEntries()); int firstDetStrip = NStrips; - for (Int_t j = 0; j < nEnt; j++) { + for (Int_t j = FstHitinTs; j < nEntSts; j++) { Int_t hitIndex = 0; - if (fTimesliceMode) - hitIndex = j; + if (fTimesliceMode) hitIndex = j; else hitIndex = (event ? event->GetIndex(ECbmDataType::kStsHit, j) : j); - CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex)); + int hitIndexSort = 0; + if (fTimesliceMode) hitIndexSort = StsIndex[hitIndex]; + else + hitIndexSort = j; + + CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndexSort)); TmpHit th; { - CbmStsHit* mh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex)); - th.ExtIndex = hitIndex; + CbmStsHit* mh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndexSort)); + th.ExtIndex = hitIndexSort; th.Det = 1; - th.iStation = NMvdStations - + CbmStsSetup::Instance()->GetStationNumber( - mh->GetAddress()); //mh->GetStationNr() - 1; + th.iStation = + NMvdStations + CbmStsSetup::Instance()->GetStationNumber(mh->GetAddress()); //mh->GetStationNr() - 1; th.iStripF = firstDetStrip + mh->GetFrontClusterId(); th.iStripB = firstDetStrip + mh->GetBackClusterId(); @@ -519,6 +502,15 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.time = mh->GetTime(); th.t_er = mh->GetTimeError(); + if (fTimesliceMode) th.id = hitIndex; + else + th.id = j + nMvdHits; + + + if ((th.time > (TsStart + TsLength)) && ((nEntSts - hitIndex) > 300)) + break; /// stop if reco TS ends or few hits left + if (hitIndex == (nEntSts - 1)) newTS = 0; ///stop while if all hits are processed + th.iStripF += nMvdHits; th.iStripB += nMvdHits; @@ -538,10 +530,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.dv = mh->GetDv(); L1Station& st = algo->vStations[th.iStation]; - th.u_front = - th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - th.u_back = - th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } th.iMC = -1; @@ -549,20 +539,16 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fPerformance) { if (listStsClusterMatch) { - const CbmMatch* frontClusterMatch = static_cast<const CbmMatch*>( - listStsClusterMatch->At(sh->GetFrontClusterId())); - const CbmMatch* backClusterMatch = static_cast<const CbmMatch*>( - listStsClusterMatch->At(sh->GetBackClusterId())); + const CbmMatch* frontClusterMatch = + static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetFrontClusterId())); + const CbmMatch* backClusterMatch = + static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetBackClusterId())); CbmMatch stsHitMatch; - for (Int_t iFrontLink = 0; - iFrontLink < frontClusterMatch->GetNofLinks(); - iFrontLink++) { + for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) { const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink); - for (Int_t iBackLink = 0; - iBackLink < backClusterMatch->GetNofLinks(); - iBackLink++) { + for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); iBackLink++) { const CbmLink& backLink = backClusterMatch->GetLink(iBackLink); if (frontLink == backLink) { stsHitMatch.AddLink(frontLink); @@ -580,8 +566,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry(); // if(!fTimesliceMode) //TODO Fix the event number in links // iEvent+=1; - Int_t iIndex = - stsHitMatch.GetLink(iLink).GetIndex() + nMvdPoints; + Int_t iIndex = stsHitMatch.GetLink(iLink).GetIndex() + nMvdPoints; Double_t dtrck = dFEI(iFile, iEvent, iIndex); DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck); if (trk_it == dFEI2vMCPoints.end()) continue; @@ -589,7 +574,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { } } } - } else + } + else iMC = sh->GetRefId(); // TODO1: don't need this with FairLinks } //fPerformance @@ -614,23 +600,21 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { TmpHit th; { - CbmMuchPixelHit* mh = - static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(j)); + CbmMuchPixelHit* mh = static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(j)); th.ExtIndex = j; th.Det = 2; - Int_t stationNumber = - CbmMuchGeoScheme::GetStationIndex(mh->GetAddress()); - Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(mh->GetAddress()); + Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(mh->GetAddress()); + Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(mh->GetAddress()); int DetId = stationNumber * 3 + layerNumber; th.iStation = DetId + NMvdStations + NStsStations; //Get time - th.time = mh->GetTime() - 17; + th.time = mh->GetTime() - 14.5; th.t_er = mh->GetTimeError(); @@ -640,7 +624,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; } th.x = mh->GetX(); - th.y = mh->GetY(); + th.y = mh->GetY() - 0.5; th.z = mh->GetZ(); th.dx = mh->GetDx(); @@ -652,10 +636,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { L1Station& st = algo->vStations[th.iStation]; - th.u_front = - th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - th.u_back = - th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } th.iMC = -1; int iMC = -1; @@ -663,8 +645,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fPerformance) { if (listMuchHitMatches) { - CbmMatch* matchHitMatch = - L1_DYNAMIC_CAST<CbmMatch*>(listMuchHitMatches->At(j)); + CbmMatch* matchHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(listMuchHitMatches->At(j)); for (Int_t iLink = 0; iLink < matchHitMatch->GetNofLinks(); iLink++) { @@ -679,32 +660,32 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck); if (trk_it == dFEI2vMCPoints.end()) continue; th.iMC = trk_it->second; - // th.track = vMCPoints[th.iMC].ID; - // th.qp = vMCPoints[iMC].p; - // if(matchHitMatch -> GetNofLinks() == 0) continue; - // Float_t bestWeight = 0.f; - // Float_t totalWeight = 0.f; - // int iMCPoint = -1; - // CbmLink link; - - // CbmMuchPoint* pt = (CbmMuchPoint*) fMuchPoints->Get( - // matchHitMatch->GetLink(0).GetFile(), - // matchHitMatch->GetLink(0).GetEntry(), - // matchHitMatch->GetLink(0).GetIndex()); - // double mom = sqrt(pt->GetPxOut()*pt->GetPxOut()+pt->GetPyOut()*pt->GetPyOut()+pt->GetPzOut()*pt->GetPzOut()); - // th.p = mom; - // th.q = pt->GetTrackID();//(L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(link.GetFile(),link.GetEntry(), pt->GetTrackID()) ))->GetCharge(); - - // th.x = pt->GetX( th.z ) + gRandom->Gaus(0,th.dx); + // th.track = vMCPoints[th.iMC].ID; + // th.qp = vMCPoints[iMC].p; + // if(matchHitMatch -> GetNofLinks() == 0) continue; + // Float_t bestWeight = 0.f; + // Float_t totalWeight = 0.f; + // int iMCPoint = -1; + // CbmLink link; + // + // CbmMuchPoint* pt = (CbmMuchPoint*) fMuchPoints->Get( + // matchHitMatch->GetLink(0).GetFile(), + // matchHitMatch->GetLink(0).GetEntry(), + // matchHitMatch->GetLink(0).GetIndex()); + // // double mom = sqrt(pt->GetPxOut()*pt->GetPxOut()+pt->GetPyOut()*pt->GetPyOut()+pt->GetPzOut()*pt->GetPzOut()); + // // th.p = mom; + // // th.q = pt->GetTrackID();//(L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(link.GetFile(),link.GetEntry(), pt->GetTrackID()) ))->GetCharge(); + // + // th.x = pt->GetX( th.z );// + gRandom->Gaus(0,th.dx); + // + // th.y = pt->GetY(th.z);// + gRandom->Gaus(0,th.dy); + // th.time = pt->GetTime(); //+ gRandom->Gaus(0,th.t_er); // - // th.y = pt->GetY(th.z) + gRandom->Gaus(0,th.dy); - // th.time = pt->GetTime(); //+ gRandom->Gaus(0,th.t_er); - - // L1Station& st = algo->vStations[th.iStation]; - // th.u_front = - // th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - // th.u_back = - // th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + // L1Station& st = algo->vStations[th.iStation]; + // th.u_front = + // th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + // th.u_back = + // th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } } } @@ -726,22 +707,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(j)); th.ExtIndex = j; th.Det = 3; - /* num variable not used - int num =mh->GetPlaneId(); - - if ((mh->GetPlaneId())==0) num = 0; - if ((mh->GetPlaneId())==1) num = 1; - if ((mh->GetPlaneId())==2) num = 2; - if ((mh->GetPlaneId())==3) num = 3; - if ((mh->GetPlaneId())==4) num = 4; -// if (num == 1) continue; -// if (num == 3) continue; - */ - // for (int k=0; k <TrdHitsOnStation[num+1].size(); k++ ){ - - th.iStation = - NMvdStations + mh->GetPlaneId() + NStsStations + NMuchStations; + th.iStation = NMvdStations + mh->GetPlaneId() + NStsStations + NMuchStations; th.time = mh->GetTime(); @@ -758,6 +725,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.x = pos.X(); th.y = pos.Y(); + th.z = pos.Z(); th.dx = fabs(mh->GetDx()); @@ -767,37 +735,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.du = fabs(mh->GetDx()); th.dv = fabs(mh->GetDy()); - // CbmTrdHit mh1 = TrdHitsOnStation[num+1][k]; - // - // TVector3 pos1, err1; - // mh1.Position(pos1); - // mh1.PositionError(err1); - // - // if (th.dx>th.dy){ - // - // th.x = pos1.X(); - // - // - // th.dx = mh1.GetDx(); - // - // - // th.du = mh1.GetDx(); - // - // } - // - // if (th.dy>th.dx){ - // - // th.y = pos1.Y(); - // - // th.dy = mh1.GetDy(); - // - // th.dv = mh1.GetDy(); - // } - L1Station& st = algo->vStations[th.iStation]; - th.u_front = - th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; th.iMC = -1; int iMC = -1; @@ -806,53 +746,38 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fTrdHitMatches) { - CbmMatch* trdHitMatch = - L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j)); - - // CbmMatch *trdHitMatch1 = L1_DYNAMIC_CAST<CbmMatch*>( fTrdHitMatches->At(TrdHitsOnStationIndex[num+1][k]) ); + CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j)); if (trdHitMatch->GetNofLinks() > 0) if (trdHitMatch->GetLink(0).GetIndex() < nTrdPoints) { iMC = trdHitMatch->GetLink(0).GetIndex(); th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints; - // th.track = vMCPoints[th.iMC].ID; - - // CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get( - // trdHitMatch->GetLink(0).GetFile(), - // trdHitMatch->GetLink(0).GetEntry(), - // trdHitMatch->GetLink(0).GetIndex()); - - // float min = 0.1; - - - // if( trdHitMatch1->GetNofLinks() > 0 ) - // if( trdHitMatch1->GetLink(0).GetIndex() < nTrdPoints ) - // { - // - // CbmTrdPoint* pt1 = (CbmTrdPoint*) fTrdPoints->Get(trdHitMatch1->GetLink(0).GetFile(),trdHitMatch1->GetLink(0).GetEntry(),trdHitMatch1->GetLink(0).GetIndex()); + // th.track = vMCPoints[th.iMC].ID; // - // if (mh->GetDx()>mh->GetDy()){ + // CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get( + // trdHitMatch->GetLink(0).GetFile(), + // trdHitMatch->GetLink(0).GetEntry(), + // trdHitMatch->GetLink(0).GetIndex()); // - // th.dx = mh1.GetDx(); - // th.du = mh1.GetDx(); - // th.x = pt1->GetXOut()+ gRandom->Gaus(0,th.dx); - // } + // float min = 0.1; // - // if (mh->GetDy()>mh->GetDx()){ // - // th.dy = mh1.GetDy(); - // th.dv = mh1.GetDy(); - // th.y = pt1->GetYOut()+ gRandom->Gaus(0,th.dy); - // } + // th.dx = min; + // th.du = min; + // th.x = 0.5*(pt->GetXOut()+pt->GetXIn())+ gRandom->Gaus(0,th.dx); // + // th.dy = min; + // th.dv = min; + // th.y = 0.5*(pt->GetYOut()+pt->GetYIn())+ gRandom->Gaus(0,th.dy); // - // L1Station& st = algo->vStations[th.iStation]; - // th.u_front = - // th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - // th.u_back = - // th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + // th.time = pt->GetTime(); // + // L1Station& st = algo->vStations[th.iStation]; + // th.u_front = + // th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + // th.u_back = + // th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; } } } @@ -873,12 +798,12 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j)); th.ExtIndex = j; th.Det = 4; - th.iStation = NMvdStations + mh->GetPlaneId() + NStsStations - + NMuchStations + NTrdStations; + + th.iStation = fTofDigiBdfPar->GetTrackingStation(mh) + NMvdStations + NStsStations + NMuchStations + NTrdStations; th.time = mh->GetTime(); - th.t_er = mh->GetTimeError(); + th.t_er = 5; //mh->GetTimeError(); th.dx = mh->GetDx(); th.dy = mh->GetDy(); @@ -898,11 +823,11 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { th.x = pos.X(); th.y = pos.Y(); + th.z = pos.Z(); L1Station& st = algo->vStations[th.iStation]; - th.u_front = - th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; - th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; + th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0]; + th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0]; th.iMC = -1; @@ -918,8 +843,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { Int_t iFile = ToFPointsMatch[j]->GetFile(); Int_t iEvent = ToFPointsMatch[j]->GetEntry(); - Int_t iIndex = ToFPointsMatch[j]->GetIndex() + nMvdPoints + nStsPoints - + nMuchPoints + nTrdPoints; + Int_t iIndex = ToFPointsMatch[j]->GetIndex() + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints; Double_t dtrck = dFEI(iFile, iEvent, iIndex); DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck); @@ -959,7 +883,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { fData_->NStsStrips = NStrips; fData_->vSFlag.resize(NStrips, 0); for (int ih = 0; ih < nHits; ih++) { - TmpHit& th = tmpHits[ih]; + TmpHit& th = tmpHits[ih]; char flag = th.iStation * 4; fData_->vSFlag[th.iStripF] = flag; fData_->vSFlag[th.iStripB] = flag; @@ -969,8 +893,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { // -- save hits -- int nEffHits = 0; - vector<float> - vStsZPos_temp; // temp array for unsorted z positions of detectors segments + SortedIndex.resize(max(nEntSts, nHits)); + + vector<float> vStsZPos_temp; // temp array for unsorted z positions of detectors segments for (int i = 0; i < nHits; i++) { TmpHit& th = tmpHits[i]; @@ -985,6 +910,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { s.dxy = th.dxy; s.time = th.time; + SortedIndex[th.id] = i; + assert(th.iStripF >= 0 || th.iStripF < NStrips); assert(th.iStripB >= 0 || th.iStripB < NStrips); @@ -992,6 +919,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { h.f = th.iStripF; h.b = th.iStripB; + h.ID = th.id; + h.t_reco = th.time; h.t_er = th.t_er; // h.track = th.track; @@ -1011,43 +940,36 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { int ist = th.iStation; if (ist < NMvdStations) { #ifdef MVDIDEALHITS - CbmMvdPoint* point = - L1_DYNAMIC_CAST<CbmMvdPoint*>(listMvdPts->At(-s.ExtIndex - 1)); - z_tmp = 0.5 * (point->GetZOut() + point->GetZ()); + CbmMvdPoint* point = L1_DYNAMIC_CAST<CbmMvdPoint*>(listMvdPts->At(-s.ExtIndex - 1)); + z_tmp = 0.5 * (point->GetZOut() + point->GetZ()); #else - CbmMvdHit* mh_m = - L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(-s.ExtIndex - 1)); - z_tmp = mh_m->GetZ(); + CbmMvdHit* mh_m = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(-s.ExtIndex - 1)); + z_tmp = mh_m->GetZ(); #endif } if ((ist >= NMvdStations) && (ist < (NStsStations + NMvdStations))) { #ifdef STSIDEALHITS - CbmStsPoint* point = - L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex)); - z_tmp = 0.5 * (point->GetZOut() + point->GetZIn()); + CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex)); + z_tmp = 0.5 * (point->GetZOut() + point->GetZIn()); #else - CbmStsHit* mh_m = - L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(s.ExtIndex)); - z_tmp = mh_m->GetZ(); + CbmStsHit* mh_m = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(s.ExtIndex)); + z_tmp = mh_m->GetZ(); #endif } - if ((ist >= NStsStations + NMvdStations) - && (ist < (NStsStations + NMvdStations + NMuchStations))) { + if ((ist >= NStsStations + NMvdStations) && (ist < (NStsStations + NMvdStations + NMuchStations))) { //#ifdef STSIDEALHITS // CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex)); // z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() ); //#else - CbmMuchPixelHit* mh = - static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(s.ExtIndex)); - z_tmp = mh->GetZ(); + CbmMuchPixelHit* mh = static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(s.ExtIndex)); + z_tmp = mh->GetZ(); //#endif } if ((ist >= NStsStations + NMvdStations + NMuchStations) - && (ist - < (NStsStations + NMvdStations + NMuchStations + NTrdStations))) { + && (ist < (NStsStations + NMvdStations + NMuchStations + NTrdStations))) { //#ifdef STSIDEALHITS // CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex)); // z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() ); @@ -1057,10 +979,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { //#endif } - if ((ist >= NStsStations + NMvdStations + NMuchStations + NTrdStations - + NTOFStation) - && (ist - < (NStsStations + NMvdStations + NMuchStations + NTrdStations))) { + if ((ist >= NStsStations + NMvdStations + NMuchStations + NTrdStations) + && (ist < (NStsStations + NMvdStations + NMuchStations + NTrdStations + NTOFStation))) { //#ifdef STSIDEALHITS // CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex)); // z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() ); @@ -1083,20 +1003,23 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { } // save hit - vStsHits.push_back( - CbmL1StsHit(fData->vStsHits.size(), th.ExtIndex, th.Det)); + vStsHits.push_back(CbmL1StsHit(fData->vStsHits.size(), th.ExtIndex, th.Det)); vStsHits[vStsHits.size() - 1].x = th.x; vStsHits[vStsHits.size() - 1].y = th.y; vStsHits[vStsHits.size() - 1].t = th.time; + vStsHits[vStsHits.size() - 1].ID = th.id; + + vStsHits[vStsHits.size() - 1].f = th.iStripF; + vStsHits[vStsHits.size() - 1].b = th.iStripB; + fData_->vStsHits.push_back(h); int sta = th.iStation; - if (fData_->StsHitsStartIndex[sta] == static_cast<THitI>(-1)) - fData_->StsHitsStartIndex[sta] = nEffHits; + if (fData_->StsHitsStartIndex[sta] == static_cast<THitI>(-1)) fData_->StsHitsStartIndex[sta] = nEffHits; nEffHits++; fData_->StsHitsStopIndex[sta] = nEffHits; @@ -1160,12 +1083,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { if (fVerbose >= 10) cout << "ReadEvent: z-pos are saved." << endl; - algo->SetData(fData_->GetStsHits(), - fData_->GetNStsStrips(), - fData_->GetStsZPos(), - fData_->GetSFlag(), - fData_->GetStsHitsStartIndex(), - fData_->GetStsHitsStopIndex()); + algo->SetData(fData_->GetStsHits(), fData_->GetNStsStrips(), fData_->GetStsZPos(), fData_->GetSFlag(), + fData_->GetStsHitsStartIndex(), fData_->GetStsHitsStopIndex(), max(nEntSts, nHits)); if (fPerformance) { @@ -1179,7 +1098,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) { } // void CbmL1::ReadEvent() -void CbmL1::Fill_vMCTracks() { +void CbmL1::Fill_vMCTracks() +{ PrimVtx.MC_ID = 999; { CbmL1Vtx Vtxcurr; @@ -1187,9 +1107,7 @@ void CbmL1::Fill_vMCTracks() { vMCTracks.clear(); - for (DFSET::iterator set_it = vFileEvent.begin(); - set_it != vFileEvent.end(); - ++set_it) { + for (DFSET::iterator set_it = vFileEvent.begin(); set_it != vFileEvent.end(); ++set_it) { Int_t iFile = set_it->first; Int_t iEvent = set_it->second; @@ -1198,8 +1116,7 @@ void CbmL1::Fill_vMCTracks() { for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) { - CbmMCTrack* MCTrack = - L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack)); + CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack)); if (!MCTrack) continue; int mother_ID = MCTrack->GetMotherId(); @@ -1213,11 +1130,11 @@ void CbmL1::Fill_vMCTracks() { Int_t pdg = MCTrack->GetPdgCode(); Double_t q = 1, mass = 0.; - if (pdg < 9999999 - && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { + if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) { q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0; mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); - } else + } + else q = 0; Int_t IND_Track = vMCTracks.size(); //or iMCTrack? CbmL1MCTrack T(mass, q, vr, vp, IND_Track, mother_ID, pdg); @@ -1228,12 +1145,10 @@ void CbmL1::Fill_vMCTracks() { vMCTracks.push_back(T); // Double_t dtrck =dFEI(iFile,iEvent,iMCTrack); - dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack), - vMCTracks.size() - 1)); + dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack), vMCTracks.size() - 1)); - if (T.mother_ID < 0) { // vertex track - if (PrimVtx.MC_ID == 999 - || fabs(T.z - Vtxcurr.MC_z) > 1.e-7) { // new vertex + if (T.mother_ID < 0) { // vertex track + if (PrimVtx.MC_ID == 999 || fabs(T.z - Vtxcurr.MC_z) > 1.e-7) { // new vertex if (nvtrackscurr > nvtracks) { PrimVtx = Vtxcurr; nvtracks = nvtrackscurr; @@ -1243,7 +1158,8 @@ void CbmL1::Fill_vMCTracks() { Vtxcurr.MC_z = T.z; Vtxcurr.MC_ID = T.mother_ID; nvtrackscurr = 1; - } else + } + else nvtrackscurr++; } } //iTracks @@ -1254,17 +1170,13 @@ void CbmL1::Fill_vMCTracks() { if (fVerbose && PrimVtx.MC_ID == 999) cout << "No primary vertex !!!" << endl; } //Fill_vMCTracks -bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, - int iPoint, - int file, - int event, - int MVD) { +bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int MVD) +{ TVector3 xyzI, PI, xyzO, PO; Int_t mcID = -1; Double_t time = 0.f; if (MVD == 1) { - CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>( - fMvdPoints->Get(file, event, iPoint)); // file, event, object + CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(fMvdPoints->Get(file, event, iPoint)); // file, event, object //CbmMvdPoint *pt = L1_DYNAMIC_CAST<CbmMvdPoint*> (Point); if (!pt) return 1; @@ -1276,8 +1188,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, time = pt->GetTime(); } if (MVD == 0) { - CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>( - fStsPoints->Get(file, event, iPoint)); // file, event, object + CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>(fStsPoints->Get(file, event, iPoint)); // file, event, object if (!pt) return 1; // if ( fTimesliceMode ) // { @@ -1301,14 +1212,12 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, if (MVD == 2) { - CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>( - fMuchPoints->Get(file, event, iPoint)); // file, event, object + CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>(fMuchPoints->Get(file, event, iPoint)); // file, event, object if (!pt) return 1; if (fTimesliceMode) { - Double_t StartTime = fTimeSlice->GetStartTime(); - Double_t EndTime = fTimeSlice->GetEndTime(); - Double_t Time_MC_point = - pt->GetTime() + fEventList->GetEventTime(event, file); + Double_t StartTime = fTimeSlice->GetStartTime(); + Double_t EndTime = fTimeSlice->GetEndTime(); + Double_t Time_MC_point = pt->GetTime() + fEventList->GetEventTime(event, file); if (Time_MC_point < StartTime) return 1; if (Time_MC_point > EndTime) return 1; @@ -1324,8 +1233,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, if (MVD == 3) { - CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>( - fTrdPoints->Get(file, event, iPoint)); // file, event, object + CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>(fTrdPoints->Get(file, event, iPoint)); // file, event, object if (!pt) return 1; if (fTimesliceMode) { @@ -1349,14 +1257,12 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, } if (MVD == 4) { - CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>( - fTofPoints->Get(file, event, iPoint)); // file, event, object + CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fTofPoints->Get(file, event, iPoint)); // file, event, object if (!pt) return 1; if (fTimesliceMode) { - Double_t StartTime = fTimeSlice->GetStartTime(); - Double_t EndTime = fTimeSlice->GetEndTime(); - Double_t Time_MC_point = - pt->GetTime() + fEventList->GetEventTime(event, file); + Double_t StartTime = fTimeSlice->GetStartTime(); + Double_t EndTime = fTimeSlice->GetEndTime(); + Double_t Time_MC_point = pt->GetTime() + fEventList->GetEventTime(event, file); if (Time_MC_point < StartTime) return 1; if (Time_MC_point > EndTime) return 1; @@ -1390,17 +1296,16 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, MC->pxOut = PO.X(); MC->pyOut = PO.Y(); MC->pzOut = PO.Z(); - MC->p = sqrt(fabs(MC->px * MC->px + MC->py * MC->py + MC->pz * MC->pz)); - MC->ID = mcID; - MC->pointId = iPoint; - MC->file = file; - MC->event = event; + MC->p = sqrt(fabs(MC->px * MC->px + MC->py * MC->py + MC->pz * MC->pz)); + MC->ID = mcID; + MC->pointId = iPoint; + MC->file = file; + MC->event = event; MC->time = time; if (MC->ID < 0) return 1; - CbmMCTrack* MCTrack = - L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(file, event, MC->ID)); + CbmMCTrack* MCTrack = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(file, event, MC->ID)); if (!MCTrack) return 1; MC->pdg = MCTrack->GetPdgCode(); MC->mother_ID = MCTrack->GetMotherId(); @@ -1414,21 +1319,19 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, return 0; } -bool CbmL1::ReadMCPoint(CbmL1MCPoint* /*MC*/, int /*iPoint*/, int /*MVD*/) { - return 0; -} +bool CbmL1::ReadMCPoint(CbmL1MCPoint* /*MC*/, int /*iPoint*/, int /*MVD*/) { return 0; } /// Procedure for match hits and MCPoints. /// Read information about correspondence between hits and mcpoints and fill CbmL1MCPoint::hitIds and CbmL1StsHit::mcPointIds arrays /// should be called after fill of algo -void CbmL1::HitMatch() { +void CbmL1::HitMatch() +{ const int NHits = vStsHits.size(); for (int iH = 0; iH < NHits; iH++) { CbmL1StsHit& hit = vStsHits[iH]; if (hit.Det == 1) { - CbmStsHit* sh = - L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(vStsHits[iH].extIndex)); + CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(vStsHits[iH].extIndex)); int iP = -1; @@ -1437,38 +1340,32 @@ void CbmL1::HitMatch() { if (listStsClusterMatch) { - const CbmMatch* frontClusterMatch = static_cast<const CbmMatch*>( - listStsClusterMatch->At(sh->GetFrontClusterId())); - const CbmMatch* backClusterMatch = static_cast<const CbmMatch*>( - listStsClusterMatch->At(sh->GetBackClusterId())); + const CbmMatch* frontClusterMatch = + static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetFrontClusterId())); + const CbmMatch* backClusterMatch = + static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetBackClusterId())); CbmMatch stsHitMatch; Float_t Sum_of_front = 0; Float_t Sum_of_back = 0; - for (Int_t iFrontLink = 0; - iFrontLink < frontClusterMatch->GetNofLinks(); - iFrontLink++) { + for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) { const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink); Sum_of_front = Sum_of_front + frontLink.GetWeight(); } - for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); - iBackLink++) { + for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); iBackLink++) { const CbmLink& backLink = backClusterMatch->GetLink(iBackLink); Sum_of_back = Sum_of_back + backLink.GetWeight(); } - for (Int_t iFrontLink = 0; - iFrontLink < frontClusterMatch->GetNofLinks(); - iFrontLink++) { + for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) { const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink); // Float_t Fraction_front = frontLink.GetWeight()/Sum_of_front; - for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); - iBackLink++) { + for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks(); iBackLink++) { const CbmLink& backLink = backClusterMatch->GetLink(iBackLink); // Float_t Fraction_back = backLink.GetWeight()/Sum_of_back; @@ -1514,7 +1411,6 @@ void CbmL1::HitMatch() { if (iP >= 0) { - for (unsigned int i = 0; i < iEvent1.size(); i++) { hit.event = iEvent1[i]; } @@ -1522,7 +1418,8 @@ void CbmL1::HitMatch() { hit.event = vMCPoints[iP].event; hit.mcPointIds.push_back(iP); vMCPoints[iP].hitIds.push_back(iH); - } else { + } + else { hit.event = -1; int idPoint = vHitMCRef[iH]; if (idPoint >= 0) { diff --git a/reco/L1/CbmL1StsHit.h b/reco/L1/CbmL1StsHit.h index 8085f442ad41f7c4748006c332e3d179898c8028..56d42988400d7a46beeb970013ea1c52e1044bf6 100644 --- a/reco/L1/CbmL1StsHit.h +++ b/reco/L1/CbmL1StsHit.h @@ -15,6 +15,9 @@ public: , x(0.) , y(0.) , t(0.) + , f(0.) + , b(0.) + , ID(0.) , file(0) , event(0) {}; CbmL1StsHit(int hitId_, int extIndex_, int Det_) @@ -25,16 +28,19 @@ public: , x(0.) , y(0.) , t(0.) + , f(0.) + , b(0.) + , ID(0.) , file(0) , event(0) {}; - int - hitId; // index of L1StsHit in algo->vStsHits array. Should be equal to index of this in L1->vStsHits + int hitId; // index of L1StsHit in algo->vStsHits array. Should be equal to index of this in L1->vStsHits int extIndex; // index of hit in the TClonesArray array int Det; vector<int> mcPointIds; // indices of CbmL1MCPoint in L1->vMCPoints array float x, y, t; - + int f, b; + int ID; int file; int event; }; diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h index 82de45f4270e2b6d29fdc17056356729b1f76c3d..7c846767898878744782ff02053244d053c0740a 100644 --- a/reco/L1/CbmL1Track.h +++ b/reco/L1/CbmL1Track.h @@ -30,14 +30,7 @@ class CbmL1MCTrack; class CbmL1Track : public CbmL1TrackPar { public: - CbmL1Track() - : StsHits() - , nStations(0) - , index(0) - , fTrackTime(0.) - , hitMap() - , mcTracks() - , maxPurity(-1) {} + CbmL1Track() : StsHits(), nStations(0), index(0), fTrackTime(0.), hitMap(), mcTracks(), maxPurity(-1) {} int GetNOfHits() { return StsHits.size(); } @@ -50,12 +43,8 @@ public: void SetMaxPurity(double maxPurity_) { maxPurity = maxPurity_; } double GetMaxPurity() { return maxPurity; } - static bool compareChi2(const CbmL1Track& a, const CbmL1Track& b) { - return (a.chi2 < b.chi2); - } - static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { - return (a->chi2 < b->chi2); - } + static bool compareChi2(const CbmL1Track& a, const CbmL1Track& b) { return (a.chi2 < b.chi2); } + static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { return (a->chi2 < b->chi2); } double Tpv[7], Cpv[21]; @@ -68,13 +57,11 @@ public: double fTrackTime; - map<int, int> - hitMap; // how many hits from each mcTrack belong to current recoTrack + map<int, int> hitMap; // how many hits from each mcTrack belong to current recoTrack private: // next members filled and used in Performance - vector<CbmL1MCTrack*> - mcTracks; // array of assosiated recoTracks. Should be only one. - double maxPurity; // maximum persent of hits, which belong to one mcTrack. + vector<CbmL1MCTrack*> mcTracks; // array of assosiated recoTracks. Should be only one. + double maxPurity; // maximum persent of hits, which belong to one mcTrack. }; #endif diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx index c1d14ab25fd5a39c9f964a086ac5076da6ad6c7a..9f54f4a6daae4313952b8ee0278be47a51906b80 100644 --- a/reco/L1/L1Algo/L1Algo.cxx +++ b/reco/L1/L1Algo/L1Algo.cxx @@ -1,10 +1,10 @@ #include "L1Algo.h" + #include "L1Grid.h" #include "L1HitPoint.h" -void L1Algo::Init(const vector<fscal>& geo, - const bool UseHitErrors, - const bool mCBMmode) { +void L1Algo::Init(const vector<fscal>& geo, const bool UseHitErrors, const bool mCBMmode) +{ for (int iProc = 0; iProc < 4; iProc++) { for (int i = 0; i < 8; i++) { @@ -13,8 +13,7 @@ void L1Algo::Init(const vector<fscal>& geo, } for (int i = 0; i < 2; i++) { threadNumberToCpuMap[2 * i + 0 + 16 + iProc * 20] = 4 * i + iProc + 64; - threadNumberToCpuMap[2 * i + 1 + 16 + iProc * 20] = - 4 * i + 8 + iProc + 64; + threadNumberToCpuMap[2 * i + 1 + 16 + iProc * 20] = 4 * i + 8 + iProc + 64; } } @@ -38,8 +37,8 @@ void L1Algo::Init(const vector<fscal>& geo, B[i].y = geo[ind++]; B[i].z = geo[ind++]; #ifndef TBB2 - std::cout << "L1Algo Input Magnetic field:" << z[i][0] << " " << B[i].x[0] - << " " << B[i].y[0] << " " << B[i].z[0] << std::endl; + std::cout << "L1Algo Input Magnetic field:" << z[i][0] << " " << B[i].x[0] << " " << B[i].y[0] << " " << B[i].z[0] + << std::endl; #endif // TBB2 } vtxFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]); @@ -91,22 +90,17 @@ void L1Algo::Init(const vector<fscal>& geo, double th = b_phi - f_phi; double det = cos(th); det *= det; - st.XYInfo.C00 = - (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; - st.XYInfo.C10 = - -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; - st.XYInfo.C11 = - (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; + st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; + st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; + st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; //std::cout << "!!! det "<< det<<std::endl; - } else { + } + else { double det = c_f * s_b - s_f * c_b; det *= det; - st.XYInfo.C00 = - (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; - st.XYInfo.C10 = - -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; - st.XYInfo.C11 = - (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; + st.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det; + st.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det; + st.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det; //std::cout << "!!! det "<< det<<std::endl; } //std::cout.precision(10); @@ -137,11 +131,9 @@ void L1Algo::Init(const vector<fscal>& geo, for (int iC = 0; iC < N; iC++) st.fieldSlice.cz[iC] = geo[ind++]; #ifndef TBB2 - std::cout << " " << st.z[0] << " " << st.materialInfo.thick[0] << " " - << st.materialInfo.RL[0] << ", " << N << " field coeff." - << std::endl; - std::cout << " " << f_phi << " " << f_sigma << " " << b_phi << " " - << b_sigma << std::endl; + std::cout << " " << st.z[0] << " " << st.materialInfo.thick[0] << " " << st.materialInfo.RL[0] << ", " << N + << " field coeff." << std::endl; + std::cout << " " << f_phi << " " << f_sigma << " " << b_phi << " " << b_sigma << std::endl; #endif // TBB2 } @@ -176,17 +168,15 @@ void L1Algo::Init(const vector<fscal>& geo, } -void L1Algo::SetData(const vector<L1StsHit>& StsHits_, - int nStsStrips_, - const vector<fscal>& StsZPos_, - const vector<unsigned char>& SFlag_, - const THitI* StsHitsStartIndex_, - const THitI* StsHitsStopIndex_) { +void L1Algo::SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_, + const vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_, + const THitI* StsHitsStopIndex_, const int NhitsGlobal) +{ - vStsHits = &StsHits_; - NStsStrips = nStsStrips_; - vStsZPos = &StsZPos_; - vSFlag = &SFlag_; + vStsHits = &StsHits_; + NStsStrips = nStsStrips_; + vStsZPos = &StsZPos_; + vSFlag = &SFlag_; StsHitsStartIndex = StsHitsStartIndex_; StsHitsStopIndex = StsHitsStopIndex_; @@ -206,11 +196,11 @@ void L1Algo::SetData(const vector<L1StsHit>& StsHits_, RealIHit_v_buf2.resize(nHits); #ifdef _OPENMP - hitToBestTrackF.resize(nHits); - hitToBestTrackB.resize(nHits); + hitToBestTrackF.resize(NhitsGlobal * 2); + hitToBestTrackB.resize(NhitsGlobal * 2); #endif - vStripToTrack.resize(nHits); - vStripToTrackB.resize(nHits); + vStripToTrack.resize(NhitsGlobal * 2); + vStripToTrackB.resize(NhitsGlobal * 2); TripForHit[0].resize(nHits); TripForHit[1].resize(nHits); @@ -245,23 +235,19 @@ void L1Algo::SetData(const vector<L1StsHit>& StsHits_, } -void L1Algo::GetHitCoor(const L1StsHit& _h, fscal& _x, fscal& _y, char iS) { +void L1Algo::GetHitCoor(const L1StsHit& _h, fscal& _x, fscal& _y, char iS) +{ L1Station& sta = vStations[int(iS)]; fscal u = _h.u; fscal v = _h.v; // const fscal &z = (*vStsZPos)[_h.iz]; // fscal x, y; - _x = - (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / (*vStsZPos)[_h.iz]; - _y = - (sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v) / (*vStsZPos)[_h.iz]; + _x = (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / (*vStsZPos)[_h.iz]; + _y = (sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v) / (*vStsZPos)[_h.iz]; } -void L1Algo::GetHitCoor(const L1StsHit& _h, - fscal& _x, - fscal& _y, - fscal& _z, - const L1Station& sta) { +void L1Algo::GetHitCoor(const L1StsHit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta) +{ fscal u = _h.u; fscal v = _h.v; fscal x, y; @@ -271,12 +257,9 @@ void L1Algo::GetHitCoor(const L1StsHit& _h, _z = (*vStsZPos)[_h.iz]; } -void L1Algo::StripsToCoor(const fscal& u, - const fscal& v, - fscal& _x, - fscal& _y, - const L1Station& sta) - const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... +void L1Algo::StripsToCoor( + const fscal& u, const fscal& v, fscal& _x, fscal& _y, + const L1Station& sta) const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... { // fvec x,y; // StripsToCoor(u,v,x,y,sta); @@ -286,12 +269,9 @@ void L1Algo::StripsToCoor(const fscal& u, _y = sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v; } /// convert strip positions to coordinates -void L1Algo::StripsToCoor(const fscal& u, - const fscal& v, - fvec& _x, - fvec& _y, - const L1Station& sta) - const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... +void L1Algo::StripsToCoor( + const fscal& u, const fscal& v, fvec& _x, fvec& _y, + const L1Station& sta) const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... { // fvec x,y; // StripsToCoor(u,v,x,y,sta); @@ -301,36 +281,24 @@ void L1Algo::StripsToCoor(const fscal& u, _y = sta.yInfo.cos_phi * u + sta.yInfo.sin_phi * v; } -void L1Algo::dUdV_to_dY(const fvec& u, - const fvec& v, - fvec& _y, - const L1Station& sta) { - _y = sqrt((sta.yInfo.cos_phi * u) * (sta.yInfo.cos_phi * u) - + (sta.yInfo.sin_phi * v) * (sta.yInfo.sin_phi * v)); +void L1Algo::dUdV_to_dY(const fvec& u, const fvec& v, fvec& _y, const L1Station& sta) +{ + _y = sqrt((sta.yInfo.cos_phi * u) * (sta.yInfo.cos_phi * u) + (sta.yInfo.sin_phi * v) * (sta.yInfo.sin_phi * v)); } -void L1Algo::dUdV_to_dX(const fvec& u, - const fvec& v, - fvec& _x, - const L1Station& sta) { - _x = sqrt((sta.xInfo.sin_phi * u) * (sta.xInfo.sin_phi * u) - + (sta.xInfo.cos_phi * v) * (sta.xInfo.cos_phi * v)); +void L1Algo::dUdV_to_dX(const fvec& u, const fvec& v, fvec& _x, const L1Station& sta) +{ + _x = sqrt((sta.xInfo.sin_phi * u) * (sta.xInfo.sin_phi * u) + (sta.xInfo.cos_phi * v) * (sta.xInfo.cos_phi * v)); } -void L1Algo::dUdV_to_dXdY(const fvec& u, - const fvec& v, - fvec& _xy, - const L1Station& sta) { - _xy = ((sta.xInfo.sin_phi * u) * (sta.yInfo.cos_phi * u) - + (sta.xInfo.cos_phi * v) * (sta.yInfo.sin_phi * v)); +void L1Algo::dUdV_to_dXdY(const fvec& u, const fvec& v, fvec& _xy, const L1Station& sta) +{ + _xy = ((sta.xInfo.sin_phi * u) * (sta.yInfo.cos_phi * u) + (sta.xInfo.cos_phi * v) * (sta.yInfo.sin_phi * v)); } -void L1Algo::StripsToCoor(const fvec& u, - const fvec& v, - fvec& x, - fvec& y, - const L1Station& sta) - const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... +void L1Algo::StripsToCoor( + const fvec& u, const fvec& v, fvec& x, fvec& y, + const L1Station& sta) const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ... { // only for same-z // x = u; @@ -349,9 +317,7 @@ L1HitPoint L1Algo::CreateHitPoint(const L1StsHit& hit, char /*ista*/) return L1HitPoint(z, hit.u, hit.v, hit.du, hit.dv, time, hit.t_er); } -void L1Algo::CreateHitPoint(const L1StsHit& hit, - char /*ista*/, - L1HitPoint& point) +void L1Algo::CreateHitPoint(const L1StsHit& hit, char /*ista*/, L1HitPoint& point) /// hit and station number { // L1Station& sta = vStations[int(ista)]; diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index ed290e6353f0b2fc030546eb9463194eb4ac2e29..89c594f930bf89313ce5ea5212064fe26b6a3156 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -33,8 +33,17 @@ class L1AlgoDraw; //#define MERGE_CLONES +#include <iomanip> +#include <iostream> +#include <map> +#include <vector> + #include "L1Branch.h" #include "L1Field.h" +#include "L1Grid.h" +#include "L1HitPoint.h" +#include "L1HitsSortHelper.h" +#include "L1Portion.h" #include "L1Station.h" #include "L1StsHit.h" #include "L1Track.h" @@ -42,16 +51,6 @@ class L1AlgoDraw; #include "L1TrackParFit.h" #include "L1Triplet.h" -#include "L1Grid.h" -#include "L1HitPoint.h" -#include "L1HitsSortHelper.h" -#include "L1Portion.h" - -#include <iomanip> -#include <iostream> -#include <map> -#include <vector> - #ifdef _OPENMP #include "omp.h" #endif @@ -121,12 +120,12 @@ public: , #ifdef _OPENMP - hitToBestTrackF(TypicalSize) - , hitToBestTrackB(TypicalSize) + hitToBestTrackF(TypicalSize * 2) + , hitToBestTrackB(TypicalSize * 2) , #endif - vStripToTrack(TypicalSize) - , vStripToTrackB(TypicalSize) + vStripToTrack(TypicalSize * 2) + , vStripToTrackB(TypicalSize * 2) , //sh (), fNThreads(nThreads) @@ -281,15 +280,11 @@ public: #endif - void - Init(const vector<fscal>& geo, const bool UseHitErrors, const bool mCBMmode); + void Init(const vector<fscal>& geo, const bool UseHitErrors, const bool mCBMmode); - void SetData(const vector<L1StsHit>& StsHits_, - int nStsStrips_, - const vector<fscal>& StsZPos_, - const vector<unsigned char>& SFlag_, - const THitI* StsHitsStartIndex_, - const THitI* StsHitsStopIndex_ + void SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_, + const vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_, + const int NhitsGlobal ); void PrintHits(); @@ -298,9 +293,8 @@ public: void CATrackFinder(); /// Track fitting procedures - void - KFTrackFitter_simple(); // version, which use procedured used during the reconstruction - void L1KFTrackFitter(); // version from SIMD-KF benchmark + void KFTrackFitter_simple(); // version, which use procedured used during the reconstruction + void L1KFTrackFitter(); // version from SIMD-KF benchmark void L1KFTrackFitterMuch(); @@ -309,24 +303,24 @@ public: void SetNThreads(int n = 1) { fNThreads = n; } - enum { MaxNStations = 25 }; + enum + { + MaxNStations = 25 + }; int NStations, // number of all detector stations NMvdStations, // number of mvd stations NStsStations, NFStations; L1Station vStations[MaxNStations] _fvecalignment; // station info - vector<L1Material> fRadThick; // material for each station - - int NStsStrips; // number of strips - const vector<fscal>* vStsZPos; // all possible z-positions of hits - const vector<L1StsHit>* - vStsHits; // hits as a combination of front-, backstrips and z-position - L1Grid vGrid - [MaxNStations]; // hits as a combination of front-, backstrips and z-position + vector<L1Material> fRadThick; // material for each station + + int NStsStrips; // number of strips + const vector<fscal>* vStsZPos; // all possible z-positions of hits + const vector<L1StsHit>* vStsHits; // hits as a combination of front-, backstrips and z-position + L1Grid vGrid[MaxNStations]; // hits as a combination of front-, backstrips and z-position L1Grid vGridTime[MaxNStations]; - const vector<unsigned char>* - vSFlag; // information of hits station & using hits in tracks; + const vector<unsigned char>* vSFlag; // information of hits station & using hits in tracks; double CATime; // time of trackfinding @@ -374,7 +368,8 @@ public: /// standard sizes of the arrays - enum { + enum + { multiCoeff = 1, // central - 1, mbias @@ -387,9 +382,8 @@ public: MaxNPortion = 40 * coeff / multiCoeff, - MaxArrSize = - MaxNPortion * MaxPortionDoublets - / MaxNStations //200000, // standart size of big arrays // mas be 40000 for normal work in cbmroot! + MaxArrSize = MaxNPortion * MaxPortionDoublets + / MaxNStations //200000, // standart size of big arrays // mas be 40000 for normal work in cbmroot! }; @@ -401,10 +395,8 @@ public: vector<THitI>* RealIHitPBuf; vector<L1HitPoint>* vStsHitPointsUnused; THitI* RealIHit; // index in vStsHits indexed by index in vStsHitsUnused - THitI StsHitsUnusedStartIndex[MaxNStations + 1], - StsHitsUnusedStopIndex[MaxNStations + 1]; - THitI StsHitsUnusedStartIndexEnd[MaxNStations + 1], - StsHitsUnusedStopIndexEnd[MaxNStations + 1]; + THitI StsHitsUnusedStartIndex[MaxNStations + 1], StsHitsUnusedStopIndex[MaxNStations + 1]; + THitI StsHitsUnusedStartIndexEnd[MaxNStations + 1], StsHitsUnusedStopIndexEnd[MaxNStations + 1]; vector<int> TripForHit[2]; @@ -418,8 +410,8 @@ public: vector<THitI> fhitsl_3[nTh], fhitsm_3[nTh], fhitsr_3[nTh]; - nsL1::vector<fvec>::TSimd fu_front3[nTh], fu_back3[nTh], fz_pos3[nTh], - fTimeR[nTh], fTimeER[nTh], dx[nTh], dy[nTh], du[nTh], dv[nTh]; + nsL1::vector<fvec>::TSimd fu_front3[nTh], fu_back3[nTh], fz_pos3[nTh], fTimeR[nTh], fTimeER[nTh], dx[nTh], dy[nTh], + du[nTh], dv[nTh]; // Tindex NHits_l[MaxNStations]; @@ -432,54 +424,30 @@ public: const L1FieldRegion& GetVtxFieldRegion() const { return vtxFieldRegion; } /// ----- Hit-point-strips conversion routines ------ - void GetHitCoor(const L1StsHit& _h, - fscal& _x, - fscal& _y, - fscal& _z, - const L1Station& sta); + void GetHitCoor(const L1StsHit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta); void dUdV_to_dY(const fvec& u, const fvec& v, fvec& _y, const L1Station& sta); void dUdV_to_dX(const fvec& u, const fvec& v, fvec& _x, const L1Station& sta); - void - dUdV_to_dXdY(const fvec& u, const fvec& v, fvec& _xy, const L1Station& sta); + void dUdV_to_dXdY(const fvec& u, const fvec& v, fvec& _xy, const L1Station& sta); void GetHitCoor(const L1StsHit& _h, fscal& _x, fscal& _y, char iS); - void StripsToCoor( - const fscal& u, - const fscal& v, - fscal& x, - fscal& y, - const L1Station& sta) const; // convert strip positions to coordinates - void StripsToCoor( - const fscal& u, - const fscal& v, - fvec& x, - fvec& y, - const L1Station& sta) const; // convert strip positions to coordinates - void StripsToCoor(const fvec& u, - const fvec& v, - fvec& x, - fvec& y, - const L1Station& sta) const; - L1HitPoint - CreateHitPoint(const L1StsHit& hit, - char ista); // full the hit point by hit information. + void StripsToCoor(const fscal& u, const fscal& v, fscal& x, fscal& y, + const L1Station& sta) const; // convert strip positions to coordinates + void StripsToCoor(const fscal& u, const fscal& v, fvec& x, fvec& y, + const L1Station& sta) const; // convert strip positions to coordinates + void StripsToCoor(const fvec& u, const fvec& v, fvec& x, fvec& y, const L1Station& sta) const; + L1HitPoint CreateHitPoint(const L1StsHit& hit, + char ista); // full the hit point by hit information. void CreateHitPoint(const L1StsHit& hit, char ista, L1HitPoint& point); inline int PackIndex(const int& a, const int& b, const int& c); inline int UnPackIndex(const int& i, int& a, int& b, int& c); /// -- Flags routines -- - inline __attribute__((always_inline)) static unsigned char - GetFStation(unsigned char flag) { - return flag / 4; - } - inline __attribute__((always_inline)) static bool - GetFUsed(unsigned char flag) { - return (flag & 0x02) != 0; - } + inline __attribute__((always_inline)) static unsigned char GetFStation(unsigned char flag) { return flag / 4; } + inline __attribute__((always_inline)) static bool GetFUsed(unsigned char flag) { return (flag & 0x02) != 0; } // bool GetFUsedD ( unsigned char flag ){ return (flag&0x01)!=0; } private: @@ -487,15 +455,8 @@ private: /// ----- Subroutines used by L1Algo::CATrackFinder() ------ - void CAFindTrack(int ista, - L1Branch& best_tr, - unsigned char& best_L, - fscal& best_chi2, - const L1Triplet* curr_trip, - L1Branch& curr_tr, - unsigned char& curr_L, - fscal& curr_chi2, - unsigned char min_best_l, + void CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fscal& best_chi2, const L1Triplet* curr_trip, + L1Branch& curr_tr, unsigned char& curr_L, fscal& curr_chi2, unsigned char min_best_l, L1Branch* new_tr); @@ -505,17 +466,11 @@ private: /// dir - 0 - forward, 1 - backward /// qp0 - momentum for extrapolation /// initialize - should be params ititialized. 1 - yes. - void BranchFitterFast(const L1Branch& t, - L1TrackPar& T, - const bool dir, - const fvec qp0 = 0., + void BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0 = 0., const bool initParams = true); /// Fit track. more precise than FitterFast - void BranchFitter(const L1Branch& t, - L1TrackPar& T, - const bool dir, - const fvec qp0 = 0., + void BranchFitter(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0 = 0., const bool initParams = true); /// Find additional hits for existing track @@ -523,10 +478,7 @@ private: /// T - track params /// dir - 0 - forward, 1 - backward /// qp0 - momentum for extrapolation - void FindMoreHits(L1Branch& t, - L1TrackPar& T, - const bool dir, - const fvec qp0 = 0.0); + void FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0 = 0.0); /// Find additional hits for existing track /// return chi2 @@ -537,185 +489,109 @@ private: void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]); void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]); void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]); - void FilterTracks(fvec const r[5], - fvec const C[15], - fvec const m[5], - fvec const V[15], - fvec R[5], - fvec W[15], + void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5], fvec W[15], fvec* chi2); void CAMergeClones(); - inline __attribute__((always_inline)) void - PackLocation(unsigned int& location, - unsigned int& triplet, - unsigned int iStation, - unsigned int& thread) { + inline __attribute__((always_inline)) void PackLocation(unsigned int& location, unsigned int& triplet, + unsigned int iStation, unsigned int& thread) + { location = (triplet << 11) | (thread << 3) | iStation; } - inline __attribute__((always_inline)) void - UnPackStation(unsigned int& location, unsigned int& iStation) { + inline __attribute__((always_inline)) void UnPackStation(unsigned int& location, unsigned int& iStation) + { iStation = location & 0x7; } - inline __attribute__((always_inline)) void - UnPackThread(unsigned int& location, unsigned int& thread) { + inline __attribute__((always_inline)) void UnPackThread(unsigned int& location, unsigned int& thread) + { thread = (location >> 3) & 0xFF; } - inline __attribute__((always_inline)) void - UnPackTriplet(unsigned int& location, unsigned int& triplet) { + inline __attribute__((always_inline)) void UnPackTriplet(unsigned int& location, unsigned int& triplet) + { triplet = (location >> 11); } - inline __attribute__((always_inline)) void - SetFStation(unsigned char& flag, unsigned int iStation) { + inline __attribute__((always_inline)) void SetFStation(unsigned char& flag, unsigned int iStation) + { flag = iStation * 4 + (flag % 4); } - inline __attribute__((always_inline)) void SetFUsed(unsigned char& flag) { - flag |= 0x02; - } + inline __attribute__((always_inline)) void SetFUsed(unsigned char& flag) { flag |= 0x02; } // void SetFUsedD ( unsigned char &flag ){ flag |= 0x01; } - inline __attribute__((always_inline)) void SetFUnUsed(unsigned char& flag) { - flag &= 0xFC; - } + inline __attribute__((always_inline)) void SetFUnUsed(unsigned char& flag) { flag &= 0xFC; } // void SetFUnUsedD ( unsigned char &flag ){ flag &= 0xFE; } /// Prepare the portion of left hits data void f10( // input - Tindex start_lh, - Tindex n1_l, - L1HitPoint* StsHits_l, + Tindex start_lh, Tindex n1_l, L1HitPoint* StsHits_l, // output - fvec* u_front_l, - fvec* u_back_l, - fvec* zPos_l, - THitI* hitsl, - fvec* HitTime_l, - fvec* HitTimeEr, - fvec* Event_l, - fvec* d_x, - fvec* d_y, - fvec* d_xy, - fvec* d_u, - fvec* d_v); + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, THitI* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* Event_l, + fvec* d_x, fvec* d_y, fvec* d_xy, fvec* d_u, fvec* d_v); /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station. void f11( // input - int istal, - int istam, - Tindex n1_V, - - fvec* u_front_l, - fvec* u_back_l, - fvec* zPos_l, - fvec* HitTime_l, - fvec* HitTimeEr, + int istal, int istam, Tindex n1_V, + + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr, // output - L1TrackPar* T_1, - L1FieldRegion* fld_1, - fvec* d_x, - fvec* d_y, - fvec* d_xy, - fvec* d_u, - fvec* d_v); + L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_x, fvec* d_y, fvec* d_xy, fvec* d_u, fvec* d_v); /// Find the doublets. Reformat data in the portion of doublets. void f20( // input - Tindex n1, - L1Station& stam, - L1HitPoint* vStsHits_m, - L1TrackPar* T_1, - THitI* hitsl_1, + Tindex n1, L1Station& stam, L1HitPoint* vStsHits_m, L1TrackPar* T_1, THitI* hitsl_1, // output - Tindex& n2, - vector<THitI>& i1_2, + Tindex& n2, vector<THitI>& i1_2, #ifdef DOUB_PERFORMANCE vector<THitI>& hitsl_2, #endif // DOUB_PERFORMANCE - vector<THitI>& hitsm_2, - fvec* Event, - vector<bool>& lmDuplets); + vector<THitI>& hitsm_2, fvec* Event, vector<bool>& lmDuplets); /// Add the middle hits to parameters estimation. Propagate to right station. /// Find the triplets (right hit). Reformat data in the portion of triplets. void f30( // input - L1HitPoint* vStsHits_r, - L1Station& stam, - L1Station& star, + L1HitPoint* vStsHits_r, L1Station& stam, L1Station& star, - int istam, - int istar, - L1HitPoint* vStsHits_m, - L1TrackPar* T_1, - L1FieldRegion* fld_1, - THitI* hitsl_1, + int istam, int istar, L1HitPoint* vStsHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1, - Tindex n2, - vector<THitI>& hitsm_2, - vector<THitI>& i1_2, + Tindex n2, vector<THitI>& hitsm_2, vector<THitI>& i1_2, const vector<bool>& mrDuplets, // output - Tindex& n3, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, - nsL1::vector<fvec>::TSimd& u_front_3, - nsL1::vector<fvec>::TSimd& u_back_3, + Tindex& n3, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, vector<THitI>& hitsm_3, + vector<THitI>& hitsr_3, nsL1::vector<fvec>::TSimd& u_front_3, nsL1::vector<fvec>::TSimd& u_back_3, nsL1::vector<fvec>::TSimd& z_Pos_3, // nsL1::vector<fvec>::TSimd& dx_, // nsL1::vector<fvec>::TSimd& dy_, - nsL1::vector<fvec>::TSimd& du_, - nsL1::vector<fvec>::TSimd& dv_, - nsL1::vector<fvec>::TSimd& timeR, + nsL1::vector<fvec>::TSimd& du_, nsL1::vector<fvec>::TSimd& dv_, nsL1::vector<fvec>::TSimd& timeR, nsL1::vector<fvec>::TSimd& timeER); /// Add the right hits to parameters estimation. void f31( // input - Tindex n3_V, - L1Station& star, - nsL1::vector<fvec>::TSimd& u_front_3, - nsL1::vector<fvec>::TSimd& u_back_3, + Tindex n3_V, L1Station& star, nsL1::vector<fvec>::TSimd& u_front_3, nsL1::vector<fvec>::TSimd& u_back_3, nsL1::vector<fvec>::TSimd& z_Pos_3, // nsL1::vector<fvec>::TSimd& dx_, // nsL1::vector<fvec>::TSimd& dy_, - nsL1::vector<fvec>::TSimd& du_, - nsL1::vector<fvec>::TSimd& dv_, - nsL1::vector<fvec>::TSimd& timeR, + nsL1::vector<fvec>::TSimd& du_, nsL1::vector<fvec>::TSimd& dv_, nsL1::vector<fvec>::TSimd& timeR, nsL1::vector<fvec>::TSimd& timeER, // output nsL1::vector<L1TrackPar>::TSimd& T_3); /// Refit Triplets. void f32( // input - Tindex n3, - int istal, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, - int nIterations = 0); + Tindex n3, int istal, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, vector<THitI>& hitsm_3, + vector<THitI>& hitsr_3, int nIterations = 0); /// Select triplets. Save them into vTriplets. void f4( // input - Tindex n3, - int istal, - int istam, - int istar, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, + Tindex n3, int istal, int istam, int istar, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, + vector<THitI>& hitsm_3, vector<THitI>& hitsr_3, // output - Tindex& nstaltriplets, - vector<THitI>* hitsn_3 = 0, - vector<THitI>* hitsr_5 = 0 + Tindex& nstaltriplets, vector<THitI>* hitsn_3 = 0, vector<THitI>* hitsr_5 = 0 // #ifdef XXX // ,unsigned int &stat_n_trip @@ -731,37 +607,21 @@ private: /// Find doublets on station void DupletsStaPort( // input - int istal, - int istam, - Tindex ip, - vector<Tindex>& n_g, - Tindex* portionStopIndex_, + int istal, int istam, Tindex ip, vector<Tindex>& n_g, Tindex* portionStopIndex_, // output - L1TrackPar* T_1, - L1FieldRegion* fld_1, - THitI* hitsl_1, + L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1, vector<bool>& lmDuplets, - Tindex& n_2, - vector<THitI>& i1_2, - vector<THitI>& hitsm_2); + Tindex& n_2, vector<THitI>& i1_2, vector<THitI>& hitsm_2); /// Find triplets on station void TripletsStaPort( // input - int istal, - int istam, - int istar, - Tindex& nstaltriplets, - L1TrackPar* T_1, - L1FieldRegion* fld_1, - THitI* hitsl_1, - - Tindex& n_2, - vector<THitI>& i1_2, - vector<THitI>& hitsm_2, + int istal, int istam, int istar, Tindex& nstaltriplets, L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1, + + Tindex& n_2, vector<THitI>& i1_2, vector<THitI>& hitsm_2, const vector<bool>& mrDuplets @@ -773,57 +633,23 @@ private: /// ------ Subroutines used by L1Algo::KFTrackFitter() ------ - void GuessVec(L1TrackPar& t, - fvec* xV, - fvec* yV, - fvec* zV, - fvec* Sy, - fvec* wV, - int NHits, - fvec* zCur = 0); - void GuessVec(L1TrackParFit& t, - fvec* xV, - fvec* yV, - fvec* zV, - fvec* Sy, - fvec* wV, - int NHits, - fvec* zCur = 0, - fvec* timeV = 0, - fvec* w_time = 0); + void GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur = 0); + void GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur = 0, + fvec* timeV = 0, fvec* w_time = 0); void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st); - void - FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st); - void FilterFirst(L1TrackParFit& track, - fvec& x, - fvec& y, - fvec& t, - fvec& t_er, - L1Station& st); - - void FilterFirst(L1TrackParFit& track, - fvec& x, - fvec& y, - fvec& t, - fvec& t_er, - L1Station& st, - fvec& dx, - fvec& dy, + void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st); + void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st); + + void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy, fvec& dxy); - void FilterFirstL(L1TrackParFit& track, - fvec& x, - fvec& y, - fvec& t, - fvec& t_er, - L1Station& st, - fvec& dx, - fvec& dy, + void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy, fvec& dxy); #ifdef TBB - enum { + enum + { nthreads = 3, // number of threads nblocks = 1 // number of stations on one thread }; @@ -845,7 +671,8 @@ private: // fNFindIterations - set number of interation for trackfinding // itetation of finding: #ifdef FIND_GAPED_TRACKS - enum { + enum + { kFastPrimIter, // primary fast tracks kAllPrimIter, // primary all tracks kAllPrimJumpIter, // primary tracks with jumped triplets @@ -858,19 +685,22 @@ private: kAllSecJumpIter // secondary tracks with jumped triplets }; #ifdef TRACKS_FROM_TRIPLETS - enum { + enum + { fNFindIterations = TRACKS_FROM_TRIPLETS_ITERATION + 1 }; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter #else - enum { + enum + { fNFindIterations = 4 }; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter #endif #else - enum { + enum + { kFastPrimIter = 0, // primary fast tracks kAllPrimIter, // primary all tracks kAllSecIter, // secondary all tracks @@ -887,8 +717,7 @@ private: float TRACK_CHI2_CUT; - float - TRIPLET_CHI2_CUT; // = 5.0; // cut for selecting triplets before collecting tracks.per one DoF + float TRIPLET_CHI2_CUT; // = 5.0; // cut for selecting triplets before collecting tracks.per one DoF float DOUBLET_CHI2_CUT; float TIME_CUT1, TIME_CUT2; @@ -897,11 +726,10 @@ private: /// parameters which are different for different iterations. Set in the begin of CAL1TrackFinder - float Pick_gather; // same for attaching additional hits to track - float - PickNeighbour; // (PickNeighbour < dp/dp_error) => triplets are neighbours - fvec MaxInvMom; // max considered q/p for tracks - fvec MaxSlope; // max slope (tx\ty) in prim vertex + float Pick_gather; // same for attaching additional hits to track + float PickNeighbour; // (PickNeighbour < dp/dp_error) => triplets are neighbours + fvec MaxInvMom; // max considered q/p for tracks + fvec MaxSlope; // max slope (tx\ty) in prim vertex fvec targX, targY, targZ; // target coor L1FieldValue targB _fvecalignment; // field in the target point L1XYMeasurementInfo TargetXYInfo _fvecalignment; // target constraint [cm] diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx index 2eea1c3ed42635519419dbfe2b9808d92522cef0..978b851f231ccdaab4c54ff4fc0e134e7a5f0c5b 100644 --- a/reco/L1/L1Algo/L1CATrackFinder.cxx +++ b/reco/L1/L1Algo/L1CATrackFinder.cxx @@ -32,11 +32,10 @@ #include "pthread.h" #endif -#include "L1HitsSortHelper.h" - +#include "TRandom.h" +#include "L1HitsSortHelper.h" #include "L1Timer.h" -#include "TRandom.h" #ifdef DRAW #include "utils/L1AlgoDraw.h" #endif @@ -59,6 +58,7 @@ #include <iostream> #include <list> #include <map> + #include <stdio.h> @@ -71,23 +71,12 @@ using std::vector; /// Prepare the portion of data of left hits of a triplet: all hits except the last and the second last station will be procesesed in the algorythm, /// the data is orginesed in order to be used by SIMD inline void L1Algo::f10( // input - Tindex start_lh, - Tindex n1_l, - L1HitPoint* StsHits_l, + Tindex start_lh, Tindex n1_l, L1HitPoint* StsHits_l, // output - fvec* u_front_l, - fvec* u_back_l, - fvec* zPos_l, - THitI* hitsl, - fvec* HitTime_l, - fvec* HitTimeEr, + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, THitI* hitsl, fvec* HitTime_l, fvec* HitTimeEr, // comment unused parameters, FU, 18.01.21 - fvec* /*Event_l*/, - fvec* /*d_x*/, - fvec* /*d_y*/, - fvec* /*d_xy*/, - fvec* d_u, - fvec* d_v) { + fvec* /*Event_l*/, fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v) +{ const Tindex& end_lh = start_lh + n1_l; @@ -127,21 +116,13 @@ inline void L1Algo::f11( /// input 1st stage of singlet search int istam, /// indexes of left and middle stations of a triplet Tindex n1_V, /// - fvec* u_front_l, - fvec* u_back_l, - fvec* zPos_l, - fvec* HitTime_l, - fvec* HitTimeEr, + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr, // output // L1TrackPar *T_1, - L1TrackPar* T_1, - L1FieldRegion* fld_1, + L1TrackPar* T_1, L1FieldRegion* fld_1, // comment unused parameters, FU, 18.01.21 - fvec* /*d_x*/, - fvec* /*d_y*/, - fvec* /*d_xy*/, - fvec* d_u, - fvec* d_v) { + fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v) +{ L1Station& stal = vStations[istal]; L1Station& stam = vStations[istam]; fvec zstal = stal.z; @@ -162,9 +143,9 @@ inline void L1Algo::f11( /// input 1st stage of singlet search m_B_global _fvecalignment; // field for the next extrapolations for (int i1_V = 0; i1_V < n1_V; i1_V++) { - L1TrackPar& T = T_1[i1_V]; - L1FieldRegion& fld1 = fld_1 - [i1_V]; // by middle hit, left hit and "center" . Will be used for extrapolation to right hit + L1TrackPar& T = T_1[i1_V]; + L1FieldRegion& fld1 = + fld_1[i1_V]; // by middle hit, left hit and "center" . Will be used for extrapolation to right hit // get the magnetic field along the track. // (suppose track is straight line with origin in the target) @@ -199,55 +180,39 @@ inline void L1Algo::f11( /// input 1st stage of singlet search L1Station& stac_global = vStations[istac_global]; fvec zstac_global = stac.z; - stac.fieldSlice.GetFieldValue( - targX + tx * (zstac - targZ), targY + ty * (zstac - targZ), centB); - stal.fieldSlice.GetFieldValue( - targX + tx * (zstal - targZ), targY + ty * (zstal - targZ), l_B); + stac.fieldSlice.GetFieldValue(targX + tx * (zstac - targZ), targY + ty * (zstac - targZ), centB); + stal.fieldSlice.GetFieldValue(targX + tx * (zstal - targZ), targY + ty * (zstal - targZ), l_B); - stam_global.fieldSlice.GetFieldValue(targX + tx * (zstam_global - targZ), - targY + ty * (zstam_global - targZ), + stam_global.fieldSlice.GetFieldValue(targX + tx * (zstam_global - targZ), targY + ty * (zstam_global - targZ), m_B_global); - stal_global.fieldSlice.GetFieldValue(targX + tx * (zstal_global - targZ), - targY + ty * (zstal_global - targZ), + stal_global.fieldSlice.GetFieldValue(targX + tx * (zstal_global - targZ), targY + ty * (zstal_global - targZ), l_B_global); - stac_global.fieldSlice.GetFieldValue(targX + tx * (zstac_global - targZ), - targY + ty * (zstac_global - targZ), + stac_global.fieldSlice.GetFieldValue(targX + tx * (zstac_global - targZ), targY + ty * (zstac_global - targZ), centB_global); - if (istac != istal) - fld0.Set(l_B, stal.z, centB, stac.z, targB, targZ); + if (istac != istal) fld0.Set(l_B, stal.z, centB, stac.z, targB, targZ); else fld0.Set(l_B, zstal, targB, targZ); // estimate field for the next extrapolations - stam.fieldSlice.GetFieldValue( - targX + tx * (zstam - targZ), targY + ty * (zstam - targZ), m_B); + stam.fieldSlice.GetFieldValue(targX + tx * (zstam - targZ), targY + ty * (zstam - targZ), m_B); #define USE_3HITS #ifndef USE_3HITS - if (istac != istal) - fld1.Set(m_B, zstam, l_B, zstal, centB, zstac); + if (istac != istal) fld1.Set(m_B, zstam, l_B, zstal, centB, zstac); else fld1.Set(m_B, zstam, l_B, zstal, targB, targZ); -#else // if USE_3HITS // the best now +#else // if USE_3HITS // the best now L1FieldValue r_B _fvecalignment; L1Station& star = vStations[istam + 1]; fvec zstar = star.z; - star.fieldSlice.GetFieldValue( - targX + tx * (zstar - targZ), targY + ty * (zstar - targZ), r_B); + star.fieldSlice.GetFieldValue(targX + tx * (zstar - targZ), targY + ty * (zstar - targZ), r_B); fld1.Set(r_B, zstar, m_B, zstam, l_B, zstal); if ((istam + 1) >= NFStations) - fld1.Set(m_B_global, - zstam_global, - l_B_global, - zstal_global, - centB_global, - zstac_global); + fld1.Set(m_B_global, zstam_global, l_B_global, zstal_global, centB_global, zstac_global); #endif // USE_3HITS T.chi2 = 0.; T.NDF = 2.; - if ((isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec == kAllSecJumpIter)) - T.NDF = 0; + if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) T.NDF = 0; T.tx = tx; T.ty = ty; @@ -296,7 +261,8 @@ inline void L1Algo::f11( /// input 1st stage of singlet search J[4] = dz; J[5] = J14; L1FilterVtx(T, targX, targY, TargetXYInfo, eX, eY, J); - } else { + } + else { L1ExtrapolateLine(T, targZ); L1FilterXY(T, targX, targY, TargetXYInfo); } @@ -333,11 +299,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search if (fGlobal || fmCBMmode) // extrapolate to left hit { - if (istal <= NFStations) - L1Extrapolate0(T, zl, fld0); + if (istal <= NFStations) L1Extrapolate0(T, zl, fld0); else L1ExtrapolateLine(T, zl); - } else + } + else L1Extrapolate0(T, zl, fld0); @@ -349,20 +315,15 @@ inline void L1Algo::f11( /// input 1st stage of singlet search L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom); #endif if (ista == NMvdStations - 1) L1AddPipeMaterial(T, MaxInvMom); - } else { + } + else { #ifdef USE_RL_TABLE - L1AddMaterial(T, - fRadThick[ista].GetRadThick(T.x, T.y), - MaxInvMom, - 1, - 0.000511f * 0.000511f); + L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f); #else - L1AddMaterial( - T, vStations[ista].materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f); + L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f); #endif - if (ista == NMvdStations - 1) - L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f); + if (ista == NMvdStations - 1) L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f); } } @@ -372,11 +333,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search if (fUseHitErrors) info.sigma2 = du1 * du1; if (fGlobal || fmCBMmode) { - if (istal < NFStations) - L1Filter(T, info, u); + if (istal < NFStations) L1Filter(T, info, u); else L1FilterNoField(T, info, u); - } else + } + else L1Filter(T, info, u); info = stal.backInfo; @@ -384,11 +345,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search if (fUseHitErrors) info.sigma2 = dv1 * dv1; if (fGlobal || fmCBMmode) { - if (istal < NFStations) - L1Filter(T, info, v); + if (istal < NFStations) L1Filter(T, info, v); else L1FilterNoField(T, info, v); - } else + } + else L1Filter(T, info, v); #endif @@ -401,15 +362,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search #else L1AddMaterial(T, stal.materialInfo, MaxInvMom); #endif - if ((istam >= NMvdStations) && (istal <= NMvdStations - 1)) - L1AddPipeMaterial(T, MaxInvMom); - } else { + if ((istam >= NMvdStations) && (istal <= NMvdStations - 1)) L1AddPipeMaterial(T, MaxInvMom); + } + else { #ifdef USE_RL_TABLE - L1AddMaterial(T, - fRadThick[istal].GetRadThick(T.x, T.y), - MaxInvMom, - 1, - 0.000511f * 0.000511f); + L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f); #else L1AddMaterial(T, stal.materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f); #endif @@ -422,11 +379,11 @@ inline void L1Algo::f11( /// input 1st stage of singlet search if (fGlobal || fmCBMmode) // extrapolate to left hit { - if (istam < NFStations) - L1Extrapolate0(T, zstam, fld0); + if (istam < NFStations) L1Extrapolate0(T, zstam, fld0); else L1ExtrapolateLine(T, zstam); // TODO: fld1 doesn't work! - } else + } + else L1Extrapolate0(T, zstam, fld0); // TODO: fld1 doesn't work! } // i1_V @@ -435,23 +392,19 @@ inline void L1Algo::f11( /// input 1st stage of singlet search /// Find the doublets. Reformat data in the portion of doublets. inline void L1Algo::f20( // input - Tindex n1, - L1Station& stam, + Tindex n1, L1Station& stam, - L1HitPoint* vStsHits_m, - L1TrackPar* T_1, - THitI* hitsl_1, + L1HitPoint* vStsHits_m, L1TrackPar* T_1, THitI* hitsl_1, // output - Tindex& n2, - vector<THitI>& i1_2, + Tindex& n2, vector<THitI>& i1_2, #ifdef DOUB_PERFORMANCE vector<THitI>& hitsl_2, #endif // DOUB_PERFORMANCE vector<THitI>& hitsm_2, // comment unused parameters, FU, 18.01.21 - fvec* /*Event*/, - vector<bool>& lmDuplets) { + fvec* /*Event*/, vector<bool>& lmDuplets) +{ n2 = 0; // number of doublet for (Tindex i1 = 0; i1 < n1; ++i1) // for each singlet { @@ -462,8 +415,7 @@ inline void L1Algo::f20( // input const int n2Saved = n2; const fvec& Pick_m22 = - (DOUBLET_CHI2_CUT - - T1.chi2); // if make it bigger the found hits will be rejected later because of the chi2 cut. + (DOUBLET_CHI2_CUT - T1.chi2); // if make it bigger the found hits will be rejected later because of the chi2 cut. // Pick_m22 is not used, search for mean squared, 2nd version // -- collect possible doublets -- @@ -471,32 +423,29 @@ inline void L1Algo::f20( // input const float& timeError = T1.C55[i1_4]; const float& time = T1.t[i1_4]; - L1HitAreaTime areaTime( - vGridTime[&stam - vStations], - T1.x[i1_4] * iz, - T1.y[i1_4] * iz, - (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T1.tx))[i1_4] - * iz, - (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T1.ty))[i1_4] - * iz, - time, - sqrt(timeError) * 5); + L1HitAreaTime areaTime(vGridTime[&stam - vStations], T1.x[i1_4] * iz, T1.y[i1_4] * iz, + (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T1.tx))[i1_4] * iz, + (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T1.ty))[i1_4] * iz, time, + sqrt(timeError) * 5); THitI imh = 0; - - - while (areaTime.GetNext(imh)) - - { + int irm1 = -1; + while (true) { + if (fGlobal || fmCBMmode) { + irm1++; + if (irm1 >= (StsHitsUnusedStopIndex[&stam - vStations] - StsHitsUnusedStartIndex[&stam - vStations])) break; + imh = irm1; + } + else { + if (!areaTime.GetNext(imh)) break; + } const L1HitPoint& hitm = vStsHits_m[imh]; // check y-boundaries - if (fabs(time - hitm.time) - > sqrt(timeError + hitm.timeEr * hitm.timeEr) * 5) - continue; - if (fabs(time - hitm.time) > 40) continue; + if (fabs(time - hitm.time) > sqrt(timeError + hitm.timeEr * hitm.timeEr) * 5) continue; + if (fabs(time - hitm.time) > 100) continue; #ifdef USE_EVENT_NUMBER if ((Event[i1_V][i1_4] != hitm.n)) continue; @@ -588,15 +537,11 @@ inline void L1Algo::f20( // input #endif // DOUB_PERFORMANCE hitsm_2.push_back(imh); - TripForHit[0][hitsl_1[i1] - + StsHitsUnusedStartIndex[&stam - vStations - 1]] = 0; - TripForHit[1][hitsl_1[i1] - + StsHitsUnusedStartIndex[&stam - vStations - 1]] = 0; + TripForHit[0][hitsl_1[i1] + StsHitsUnusedStartIndex[&stam - vStations - 1]] = 0; + TripForHit[1][hitsl_1[i1] + StsHitsUnusedStartIndex[&stam - vStations - 1]] = 0; - TripForHit[0][hitsl_1[i1] - + StsHitsUnusedStartIndex[&stam - vStations - 2]] = 0; - TripForHit[1][hitsl_1[i1] - + StsHitsUnusedStartIndex[&stam - vStations - 2]] = 0; + TripForHit[0][hitsl_1[i1] + StsHitsUnusedStartIndex[&stam - vStations - 2]] = 0; + TripForHit[1][hitsl_1[i1] + StsHitsUnusedStartIndex[&stam - vStations - 2]] = 0; if (n2 > 8000) return; @@ -610,34 +555,18 @@ inline void L1Algo::f20( // input /// Add the middle hits to parameters estimation. Propagate to right station. /// Find the triplets(right hit). Reformat data in the portion of triplets. inline void L1Algo::f30( // input - L1HitPoint* vStsHits_r, - L1Station& stam, - L1Station& star, - int istam, - int istar, - L1HitPoint* vStsHits_m, - L1TrackPar* T_1, - L1FieldRegion* fld_1, - THitI* hitsl_1, - Tindex n2, - vector<THitI>& hitsm_2, - vector<THitI>& i1_2, + L1HitPoint* vStsHits_r, L1Station& stam, L1Station& star, int istam, int istar, L1HitPoint* vStsHits_m, + L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1, Tindex n2, vector<THitI>& hitsm_2, vector<THitI>& i1_2, const vector<bool>& /*mrDuplets*/, // output - Tindex& n3, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, - nsL1::vector<fvec>::TSimd& u_front_3, - nsL1::vector<fvec>::TSimd& u_back_3, + Tindex& n3, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, vector<THitI>& hitsm_3, + vector<THitI>& hitsr_3, nsL1::vector<fvec>::TSimd& u_front_3, nsL1::vector<fvec>::TSimd& u_back_3, nsL1::vector<fvec>::TSimd& z_Pos_3, // nsL1::vector<fvec>::TSimd& dx_, // nsL1::vector<fvec>::TSimd& dy_, - nsL1::vector<fvec>::TSimd& dv_, - nsL1::vector<fvec>::TSimd& du_, - nsL1::vector<fvec>::TSimd& timeR, - nsL1::vector<fvec>::TSimd& timeER) { + nsL1::vector<fvec>::TSimd& dv_, nsL1::vector<fvec>::TSimd& du_, nsL1::vector<fvec>::TSimd& timeR, + nsL1::vector<fvec>::TSimd& timeER) +{ THitI hitsl_2[fvecLen]; THitI hitsm_2_tmp[fvecLen]; fvec fvec_0; @@ -720,22 +649,22 @@ inline void L1Algo::f30( // input if (fUseHitErrors) info.sigma2 = du2 * du2; if (fGlobal || fmCBMmode) { - if (istam < NFStations) - L1Filter(T2, info, u_front_2); + if (istam < NFStations) L1Filter(T2, info, u_front_2); else L1FilterNoField(T2, info, u_front_2); - } else + } + else L1Filter(T2, info, u_front_2); info = stam.backInfo; if (fUseHitErrors) info.sigma2 = dv2 * dv2; if (fGlobal || fmCBMmode) { - if (istam < NFStations) - L1Filter(T2, info, u_back_2); + if (istam < NFStations) L1Filter(T2, info, u_back_2); else L1FilterNoField(T2, info, u_back_2); - } else + } + else L1Filter(T2, info, u_back_2); FilterTime(T2, timeM, timeMEr); @@ -746,15 +675,11 @@ inline void L1Algo::f30( // input #else L1AddMaterial(T2, stam.materialInfo, T2.qp); #endif - if ((istar >= NMvdStations) && (istam <= NMvdStations - 1)) - L1AddPipeMaterial(T2, T2.qp); - } else { + if ((istar >= NMvdStations) && (istam <= NMvdStations - 1)) L1AddPipeMaterial(T2, T2.qp); + } + else { #ifdef USE_RL_TABLE - L1AddMaterial(T2, - fRadThick[istam].GetRadThick(T2.x, T2.y), - T2.qp, - 1, - 0.000511f * 0.000511f); + L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1, 0.000511f * 0.000511f); #else L1AddMaterial(T2, stam.materialInfo, T2.qp, 1, 0.000511f * 0.000511f); #endif @@ -769,21 +694,19 @@ inline void L1Algo::f30( // input if (fGlobal || fmCBMmode) // extrapolate to the right hit station { - if (istar <= NFStations) - L1Extrapolate(T2, star.z, T2.qp, f2); + if (istar <= NFStations) L1Extrapolate(T2, star.z, T2.qp, f2); else L1ExtrapolateLine(T2, star.z); - } else + } + else L1Extrapolate(T2, star.z, T2.qp, f2); // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ---- for (Tindex i2_4 = 0; i2_4 < n2_4; ++i2_4) { - // if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 - // || T2.C33[i2_4] < 0 || T2.C44[i2_4] < 0 || T2.C55[i2_4] < 0) - // continue; - if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 - || T2.C33[i2_4] < 0) + if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0 || T2.C44[i2_4] < 0 + || T2.C55[i2_4] < 0) continue; + if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0) continue; const fvec& Pick_r22 = (TRIPLET_CHI2_CUT - T2.chi2); const float& timeError = T2.C55[i2_4]; @@ -795,21 +718,26 @@ inline void L1Algo::f30( // input if (isec == TRACKS_FROM_TRIPLETS_ITERATION) Pick_r22 = Pick_r2 + 1; #endif const fscal& iz = 1 / T2.z[i2_4]; - L1HitAreaTime area(vGridTime[&star - vStations], - T2.x[i2_4] * iz, - T2.y[i2_4] * iz, - (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) - + MaxDZ * fabs(T2.tx))[i2_4] - * iz, - (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) - + MaxDZ * fabs(T2.ty))[i2_4] - * iz, - time, + L1HitAreaTime area(vGridTime[&star - vStations], T2.x[i2_4] * iz, T2.y[i2_4] * iz, + (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T2.tx))[i2_4] * iz, + (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T2.ty))[i2_4] * iz, time, sqrt(timeError) * 5); THitI irh = 0; + int irh1 = -1; + while (true) { + if (fGlobal || fmCBMmode) { + irh1++; + if (irh1 >= (StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar])) break; + irh = irh1; + } + else { + if (!area.GetNext(irh)) break; + } - while (area.GetNext(irh)) { + + // while (area.GetNext(irh)) { + //for (int irh = 0; irh < ( StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar] ); irh++){ const L1HitPoint& hitr = vStsHits_r[irh]; #ifdef USE_EVENT_NUMBER @@ -826,10 +754,8 @@ inline void L1Algo::f30( // input L1ExtrapolateTime(T2, dz3); - if (fabs(T2.t[i2_4] - hitr.time) - > sqrt(T2.C55[i2_4] + hitr.timeEr) * 5) - continue; - if (fabs(T2.t[i2_4] - hitr.time) > 40) continue; + if (fabs(T2.t[i2_4] - hitr.time) > sqrt(T2.C55[i2_4] + hitr.timeEr) * 5) continue; + if (fabs(T2.t[i2_4] - hitr.time) > 100) continue; // // if (fabs(T2.t[i2_4]-hitr.time)>sqrt(2.9*2.9)*5) continue; @@ -842,8 +768,7 @@ inline void L1Algo::f30( // input (Pick_r22[i2_4] * (fabs( C11[i2_4] - + star.XYInfo.C11 - [i2_4]))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets + + star.XYInfo.C11[i2_4]))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets if (fUseHitErrors) { fvec dyr = 0; @@ -853,18 +778,15 @@ inline void L1Algo::f30( // input const fscal& dY = yr[i2_4] - y[i2_4]; const fscal& dY2 = dY * dY; - if (dY2 > dy_est2 && dY2 < 0) - continue; // if (yr < y_minus_new[i2_4]) continue; + if (dY2 > dy_est2 && dY2 < 0) continue; // if (yr < y_minus_new[i2_4]) continue; // check upper boundary - if (dY2 > dy_est2) - continue; // if (yr > y_plus_new [i2_4] ) continue; + if (dY2 > dy_est2) continue; // if (yr > y_plus_new [i2_4] ) continue; // check x-boundaries fvec x, C00; L1ExtrapolateXC00Line(T2, zr, x, C00); - fscal dx_est2 = - (Pick_r22[i2_4] * (fabs(C00[i2_4] + star.XYInfo.C00[i2_4]))); + fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + star.XYInfo.C00[i2_4]))); if (fUseHitErrors) { fvec dxr = 0; @@ -898,11 +820,8 @@ inline void L1Algo::f30( // input #endif if (fGlobal || fmCBMmode) - if (chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 - || C11[i2_4] < 0) - continue; - else if (chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 - || C11[i2_4] < 0 || T.C55[i2_4] < 0) + if (chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 || C11[i2_4] < 0) continue; + else if (chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 || C11[i2_4] < 0 || T.C55[i2_4] < 0) continue; // chi2_triplet < CHI2_CUT @@ -944,6 +863,8 @@ inline void L1Algo::f30( // input if (n3 > 4000) return; } + + } // i2_4 } // i2_V } // if istar @@ -951,20 +872,16 @@ inline void L1Algo::f30( // input /// Add the right hits to parameters estimation. inline void L1Algo::f31( // input - Tindex n3_V, - L1Station& star, - nsL1::vector<fvec>::TSimd& u_front_, - nsL1::vector<fvec>::TSimd& u_back_, + Tindex n3_V, L1Station& star, nsL1::vector<fvec>::TSimd& u_front_, nsL1::vector<fvec>::TSimd& u_back_, nsL1::vector<fvec>::TSimd& z_Pos, // nsL1::vector<fvec>::TSimd& dx_, // nsL1::vector<fvec>::TSimd& dy_, - nsL1::vector<fvec>::TSimd& dv_, - nsL1::vector<fvec>::TSimd& du_, - nsL1::vector<fvec>::TSimd& timeR, + nsL1::vector<fvec>::TSimd& dv_, nsL1::vector<fvec>::TSimd& du_, nsL1::vector<fvec>::TSimd& timeR, nsL1::vector<fvec>::TSimd& timeER, // output // L1TrackPar *T_3 - nsL1::vector<L1TrackPar>::TSimd& T_3) { + nsL1::vector<L1TrackPar>::TSimd& T_3) +{ for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) { fvec dz = z_Pos[i3_V] - T_3[i3_V].z; @@ -977,12 +894,12 @@ inline void L1Algo::f31( // input if (fUseHitErrors) info.sigma2 = du_[i3_V] * du_[i3_V]; if (fGlobal || fmCBMmode) { - if ((&star - vStations) < NFStations) - L1Filter(T_3[i3_V], info, u_front_[i3_V]); + if ((&star - vStations) < NFStations) L1Filter(T_3[i3_V], info, u_front_[i3_V]); else { L1FilterNoField(T_3[i3_V], info, u_front_[i3_V]); } - } else + } + else L1Filter(T_3[i3_V], info, u_front_[i3_V]); @@ -991,11 +908,11 @@ inline void L1Algo::f31( // input if (fUseHitErrors) info.sigma2 = dv_[i3_V] * dv_[i3_V]; if (fGlobal || fmCBMmode) { - if ((&star - vStations) < NFStations) - L1Filter(T_3[i3_V], info, u_back_[i3_V]); + if ((&star - vStations) < NFStations) L1Filter(T_3[i3_V], info, u_back_[i3_V]); else L1FilterNoField(T_3[i3_V], info, u_back_[i3_V]); - } else + } + else L1Filter(T_3[i3_V], info, u_back_[i3_V]); FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V]); @@ -1005,13 +922,9 @@ inline void L1Algo::f31( // input /// Refit Triplets. inline void L1Algo::f32( // input // TODO not updated after gaps introduction - Tindex n3, - int istal, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, - int nIterations) { + Tindex n3, int istal, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, vector<THitI>& hitsm_3, + vector<THitI>& hitsr_3, int nIterations) +{ const int NHits = 3; // triplets // prepare data @@ -1030,10 +943,9 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction fscal& qp0 = T3.qp[i3_4]; // prepare data - THitI ihit[NHits] = { - (*RealIHitP)[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]], - (*RealIHitP)[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]], - (*RealIHitP)[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]}; + THitI ihit[NHits] = {(*RealIHitP)[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]], + (*RealIHitP)[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]], + (*RealIHitP)[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]}; fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits]; for (int ih = 0; ih < NHits; ++ih) { @@ -1073,16 +985,11 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction L1FieldValue B[3] _fvecalignment; L1FieldRegion fld _fvecalignment; - fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), - (x[2] - x[0]) / (z[2] - z[0]), - (x[2] - x[1]) / (z[2] - z[1])}; - fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), - (y[2] - y[0]) / (z[2] - z[0]), - (y[2] - y[1]) / (z[2] - z[1])}; + fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])}; + fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])}; for (int ih = 0; ih < NHits; ++ih) { fvec dz = (sta[ih].z - z[ih]); - sta[ih].fieldSlice.GetFieldValue( - x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]); + sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]); }; fld.Set(B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z); @@ -1174,19 +1081,12 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction /// Select triplets. Save them into vTriplets. inline void L1Algo::f4( // input - Tindex n3, - int istal, - int istam, - int istar, - nsL1::vector<L1TrackPar>::TSimd& T_3, - vector<THitI>& hitsl_3, - vector<THitI>& hitsm_3, - vector<THitI>& hitsr_3, + Tindex n3, int istal, int istam, int istar, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3, + vector<THitI>& hitsm_3, vector<THitI>& hitsr_3, // output - Tindex& nstaltriplets, - vector<THitI>* /*hitsn_3*/, - vector<THitI>* /*hitsr_5*/ -) { + Tindex& nstaltriplets, vector<THitI>* /*hitsn_3*/, vector<THitI>* /*hitsr_5*/ +) +{ THitI ihitl_priv = 0; unsigned int Station = 0; @@ -1214,16 +1114,13 @@ inline void L1Algo::f4( // input #ifdef _OPENMP Thread = omp_get_thread_num(); #else - Thread = 0; + Thread = 0; #endif - TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetLevel(0); - TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]] - .SetFNeighbour(0); - TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]] - .SetNNeighbours(0); + TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetFNeighbour(0); + TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetNNeighbours(0); Triplet = nTripletsThread[istal][Thread]; @@ -1266,26 +1163,21 @@ inline void L1Algo::f4( // input const THitI& ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal]; const THitI& ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam]; const THitI& ihitr = hitsr_3[i3] + StsHitsUnusedStartIndex[istar]; - L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], - ihitl << " < " << StsHitsUnusedStopIndex[istal]); - L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], - ihitm << " < " << StsHitsUnusedStopIndex[istam]); - L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], - ihitr << " < " << StsHitsUnusedStopIndex[istar]); + L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]); + L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]); + L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]); fscal& time = T3.t[i3_4]; // int n = T3.n[i3_4]; - L1Triplet& tr1 = - TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]]; + L1Triplet& tr1 = TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]]; tr1.SetLevel(0); - tr1.Set( - ihitl, ihitm, ihitr, istal, istam, istar, 0, qp, chi2, time, Cqp, 0); + tr1.Set(ihitl, ihitm, ihitr, istal, istam, istar, 0, qp, chi2, time, Cqp, 0); tr1.tx = sqrt(fabs(T3.tx[i3_4])); tr1.ty = sqrt(fabs(T3.ty[i3_4])); @@ -1341,21 +1233,18 @@ inline void L1Algo::f4( // input inline void L1Algo::f5( // input // output - int* nlevel) { + int* nlevel) +{ #ifdef TRACKS_FROM_TRIPLETS if (isec != TRACKS_FROM_TRIPLETS_ITERATION) #endif for (int istal = NStations - 4; istal >= FIRSTCASTATION; istal--) { - for ( - int tripType = 0; tripType < 3; - tripType++) // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet + for (int tripType = 0; tripType < 3; tripType++) // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet { - if ((((isec != kFastPrimJumpIter) && (isec != kAllPrimJumpIter) - && (isec != kAllSecJumpIter)) + if ((((isec != kFastPrimJumpIter) && (isec != kAllPrimJumpIter) && (isec != kAllSecJumpIter)) && (tripType != 0)) - || (((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) - || (isec == kAllSecJumpIter)) + || (((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter)) && (tripType != 0) && (istal == NStations - 4))) continue; @@ -1382,18 +1271,14 @@ inline void L1Algo::f5( // input THitI ihitl = trip->GetLHit(); THitI ihitm = trip->GetMHit(); THitI ihitr = trip->GetRHit(); - L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], - ihitl << " < " << StsHitsUnusedStopIndex[istal]); - L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], - ihitm << " < " << StsHitsUnusedStopIndex[istam]); - L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], - ihitr << " < " << StsHitsUnusedStopIndex[istar]); + L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]); + L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]); + L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]); vector<unsigned int> neighCands; // save neighbour candidates - neighCands.reserve(8); // ~average is 1-2 for central, up to 5 + neighCands.reserve(8); // ~average is 1-2 for central, up to 5 - unsigned int Neighbours = - TripForHit[1][ihitm] - TripForHit[0][ihitm]; + unsigned int Neighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm]; for (unsigned int iN = 0; iN < Neighbours; ++iN) { // for (iN = first_triplet; iN <= last_triplet; ++iN){ @@ -1414,8 +1299,7 @@ inline void L1Algo::f5( // input fscal qp2 = tripn->GetQp(); fscal Cqp1 = trip->Cqp; fscal Cqp2 = tripn->Cqp; - if (fabs(qp - qp2) > PickNeighbour * (Cqp1 + Cqp2)) - continue; // neighbours should have same qp + if (fabs(qp - qp2) > PickNeighbour * (Cqp1 + Cqp2)) continue; // neighbours should have same qp // calculate level unsigned char jlevel = tripn->GetLevel(); @@ -1447,33 +1331,23 @@ inline void L1Algo::f5( // input inline void L1Algo:: DupletsStaPort( /// creates duplets: input: @istal - start station number, @istam - last station number, @ip - index of portion, @&n_g - number of elements in portion, @*portionStopIndex - int istal, - int istam, - Tindex ip, - vector<Tindex>& n_g, - Tindex* portionStopIndex_, + int istal, int istam, Tindex ip, vector<Tindex>& n_g, Tindex* portionStopIndex_, /// output: L1TrackPar* T_1, /// @*T_1 - singlets perameters, @*fld_1 - field aproximation, @*hitsl_1- left hits of triplets, @&lmDuplets - existance of a doublet starting from the left hit, L1FieldRegion* fld_1, /// @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index - THitI* hitsl_1, - vector<bool>& lmDuplets, - Tindex& n_2, - vector<THitI>& i1_2, - vector<THitI>& hitsm_2) { + THitI* hitsl_1, vector<bool>& lmDuplets, Tindex& n_2, vector<THitI>& i1_2, vector<THitI>& hitsm_2) +{ if (istam < NStations) { L1Station& stam = vStations[istam]; // prepare data - L1HitPoint* vStsHits_l = - &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istal]; - L1HitPoint* vStsHits_m = - &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam]; + L1HitPoint* vStsHits_l = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istal]; + L1HitPoint* vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam]; fvec u_front[Portion / fvecLen], u_back[Portion / fvecLen]; - fvec dx0[Portion / fvecLen], dy0[Portion / fvecLen], - dxy0[Portion / fvecLen]; + fvec dx0[Portion / fvecLen], dy0[Portion / fvecLen], dxy0[Portion / fvecLen]; fvec dv0[Portion / fvecLen], du0[Portion / fvecLen]; fvec zPos[Portion / fvecLen]; fvec HitTime[Portion / fvecLen]; @@ -1484,51 +1358,23 @@ inline void L1Algo:: Tindex& n1 = n_g[ip]; f10( // input - (ip - portionStopIndex_[istal + 1]) * Portion, - n1, - vStsHits_l, + (ip - portionStopIndex_[istal + 1]) * Portion, n1, vStsHits_l, // output - u_front, - u_back, - zPos, - hitsl_1, - HitTime, - HitTimeEr, - Event, - dx0, - dy0, - dxy0, - du0, - dv0); + u_front, u_back, zPos, hitsl_1, HitTime, HitTimeEr, Event, dx0, dy0, dxy0, du0, dv0); for (Tindex i = 0; i < n1; ++i) - L1_ASSERT(hitsl_1[i] < StsHitsUnusedStopIndex[istal] - - StsHitsUnusedStartIndex[istal], - hitsl_1[i] << " < " - << StsHitsUnusedStopIndex[istal] - - StsHitsUnusedStartIndex[istal]); + L1_ASSERT(hitsl_1[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal], + hitsl_1[i] << " < " << StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]); Tindex n1_V = (n1 + fvecLen - 1) / fvecLen; /// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station. - f11(istal, - istam, - n1_V, + f11(istal, istam, n1_V, - u_front, - u_back, - zPos, - HitTime, - HitTimeEr, + u_front, u_back, zPos, HitTime, HitTimeEr, // output - T_1, - fld_1, - dx0, - dy0, - dxy0, - du0, - dv0); + T_1, fld_1, dx0, dy0, dxy0, du0, dv0); /// Find the doublets. Reformat data in the portion of doublets. @@ -1537,27 +1383,17 @@ inline void L1Algo:: #endif // DOUB_PERFORMANCE f20( // input - n1, - stam, - vStsHits_m, - T_1, - hitsl_1, + n1, stam, vStsHits_m, T_1, hitsl_1, // output - n_2, - i1_2, + n_2, i1_2, #ifdef DOUB_PERFORMANCE hitsl_2, #endif // DOUB_PERFORMANCE - hitsm_2, - Event, - lmDuplets); + hitsm_2, Event, lmDuplets); for (Tindex i = 0; i < static_cast<Tindex>(hitsm_2.size()); ++i) - L1_ASSERT(hitsm_2[i] < StsHitsUnusedStopIndex[istam] - - StsHitsUnusedStartIndex[istam], - hitsm_2[i] << " " - << StsHitsUnusedStopIndex[istam] - - StsHitsUnusedStartIndex[istam]); + L1_ASSERT(hitsm_2[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam], + hitsm_2[i] << " " << StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]); #ifdef DOUB_PERFORMANCE THitI* RealIHitL = &((*RealIHitP)[StsHitsUnusedStartIndex[istal]]); @@ -1577,38 +1413,31 @@ inline void L1Algo:: inline void L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station number, @istam - middle station number, @istar - last station number, @ip - index of portion, @&n_g - numer of elements in portion, @*portionStopIndex - int istal, - int istam, - int istar, + int istal, int istam, int istar, ///@nstaltriplets - , @*portionStopIndex, @*T_1 - track parameters for singlets, @*fld_1 - field approximation for singlets, @&n_2 - number of doublets in portion /// @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index - Tindex& nstaltriplets, - L1TrackPar* T_1, - L1FieldRegion* fld_1, - THitI* hitsl_1, + Tindex& nstaltriplets, L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1, - Tindex& n_2, - vector<THitI>& i1_2, - vector<THitI>& hitsm_2, + Tindex& n_2, vector<THitI>& i1_2, vector<THitI>& hitsm_2, const vector<bool>& mrDuplets /// output: @*vTriplets_part - array of triplets, @*TripStartIndexH, @*TripStopIndexH - start/stop index of a triplet in the array -) { +) +{ if (istar < NStations) { // prepare data L1Station& stam = vStations[istam]; L1Station& star = vStations[istar]; - L1HitPoint* vStsHits_m = - &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam]; + L1HitPoint* vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam]; L1HitPoint* vStsHits_r = 0; - vStsHits_r = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istar]; + vStsHits_r = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istar]; Tindex n3 = 0, n3_V; @@ -1651,67 +1480,39 @@ L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station nu /// Find the triplets(right hit). Reformat data in the portion of triplets. + f30( // input - vStsHits_r, - stam, - star, + vStsHits_r, stam, star, - istam, - istar, - vStsHits_m, - T_1, - fld_1, - hitsl_1, + istam, istar, vStsHits_m, T_1, fld_1, hitsl_1, - n_2, - hitsm_2, - i1_2, + n_2, hitsm_2, i1_2, mrDuplets, // output - n3, - T_3, - hitsl_3, - hitsm_3, - hitsr_3, - u_front3, - u_back3, - z_pos3, + n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, // dx3, // dy3, - du3, - dv3, - timeR, - timeER); + du3, dv3, timeR, timeER); n3_V = (n3 + fvecLen - 1) / fvecLen; for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i) - L1_assert(hitsl_3[i] < StsHitsUnusedStopIndex[istal] - - StsHitsUnusedStartIndex[istal]); + L1_assert(hitsl_3[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]); for (Tindex i = 0; i < static_cast<Tindex>(hitsm_3.size()); ++i) - L1_assert(hitsm_3[i] < StsHitsUnusedStopIndex[istam] - - StsHitsUnusedStartIndex[istam]); + L1_assert(hitsm_3[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]); for (Tindex i = 0; i < static_cast<Tindex>(hitsr_3.size()); ++i) - L1_assert(hitsr_3[i] < StsHitsUnusedStopIndex[istar] - - StsHitsUnusedStartIndex[istar]); + L1_assert(hitsr_3[i] < StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar]); // if (n3 >= MaxPortionTriplets) cout << "isec: " << isec << " station: " << istal << " portion number: " << ip << " CATrackFinder: Warning: Too many Triplets created in portion" << endl; /// Add the right hits to parameters estimation. f31( // input - n3_V, - star, - u_front3, - u_back3, - z_pos3, + n3_V, star, u_front3, u_back3, z_pos3, // dx3, // dy3, - du3, - dv3, - timeR, - timeER, + du3, dv3, timeR, timeER, // output T_3); @@ -1726,29 +1527,20 @@ L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station nu for (Tindex i = 0; i < n3; ++i) { Tindex i_4 = i % 4; Tindex i_V = i / 4; - THitI iHits[3] = { - RealIHitL[hitsl_3[i]], RealIHitM[hitsm_3[i]], RealIHitR[hitsr_3[i]]}; + THitI iHits[3] = {RealIHitL[hitsl_3[i]], RealIHitM[hitsm_3[i]], RealIHitR[hitsr_3[i]]}; #ifdef PULLS - if (fL1Eff_triplets->AddOne(iHits)) - fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]); + if (fL1Eff_triplets->AddOne(iHits)) fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]); #else fL1Eff_triplets->AddOne(iHits); #endif - if (T_3[i_V].chi2[i_4] < TRIPLET_CHI2_CUT) - fL1Eff_triplets2->AddOne(iHits); + if (T_3[i_V].chi2[i_4] < TRIPLET_CHI2_CUT) fL1Eff_triplets2->AddOne(iHits); } #endif // TRIP_PERFORMANCE + /// Fill Triplets. f4( // input - n3, - istal, - istam, - istar, - T_3, - hitsl_3, - hitsm_3, - hitsr_3, + n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3, // output nstaltriplets @@ -1786,7 +1578,8 @@ L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station nu ///********************************************************************************************************************** -void L1Algo::CATrackFinder() { +void L1Algo::CATrackFinder() +{ #ifdef _OPENMP omp_set_num_threads(fNThreads); @@ -1799,19 +1592,16 @@ void L1Algo::CATrackFinder() { fL1Pulls->Init(); #endif #ifdef TRIP_PERFORMANCE - static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets_ = - new L1AlgoEfficiencyPerformance<3>(); - fL1Eff_triplets = l1Eff_triplets_; + static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets_ = new L1AlgoEfficiencyPerformance<3>(); + fL1Eff_triplets = l1Eff_triplets_; fL1Eff_triplets->Init(); - static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets2_ = - new L1AlgoEfficiencyPerformance<3>(); - fL1Eff_triplets2 = l1Eff_triplets2_; + static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets2_ = new L1AlgoEfficiencyPerformance<3>(); + fL1Eff_triplets2 = l1Eff_triplets2_; fL1Eff_triplets2->Init(); #endif #ifdef DOUB_PERFORMANCE - static L1AlgoEfficiencyPerformance<2>* l1Eff_doublets_ = - new L1AlgoEfficiencyPerformance<2>(); - fL1Eff_doublets = l1Eff_doublets_; + static L1AlgoEfficiencyPerformance<2>* l1Eff_doublets_ = new L1AlgoEfficiencyPerformance<2>(); + fL1Eff_doublets = l1Eff_doublets_; fL1Eff_doublets->Init(); #endif @@ -1865,19 +1655,16 @@ void L1Algo::CATrackFinder() { static Tindex stat_nTriplets[fNFindIterations] = {0}; static Tindex stat_nLevels[MaxNStations - 2][fNFindIterations] = {{0}}; - static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack - static Tindex stat_nTrCandidates[fNFindIterations] = {0}; + static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack + static Tindex stat_nTrCandidates[fNFindIterations] = {0}; #endif - RealIHitP = &RealIHit_v; - RealIHitPBuf = &RealIHit_v_buf; - vStsHitsUnused = - &vStsDontUsedHits_B; /// array of hits used on current iteration - vector<L1StsHit>* vStsHitsUnused_buf = - &vStsDontUsedHits_A; /// buffer for copy + RealIHitP = &RealIHit_v; + RealIHitPBuf = &RealIHit_v_buf; + vStsHitsUnused = &vStsDontUsedHits_B; /// array of hits used on current iteration + vector<L1StsHit>* vStsHitsUnused_buf = &vStsDontUsedHits_A; /// buffer for copy - vStsHitPointsUnused = - &vStsDontUsedHitsxy_B; /// array of info for hits used on current iteration + vStsHitPointsUnused = &vStsDontUsedHitsxy_B; /// array of info for hits used on current iteration vector<L1HitPoint>* vStsHitPointsUnused_buf = &vStsDontUsedHitsxy_A; NHitsIsecAll = 0; @@ -1892,13 +1679,17 @@ void L1Algo::CATrackFinder() { StsHitsUnusedStopIndex[ista] = StsHitsStopIndex[ista]; } - float lasttime = 0; + float lasttime = 0; + float starttime = std::numeric_limits<float>::max(); + ; for (int ist = 0; ist < NStations; ++ist) for (THitI ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih) { - if ((lasttime < (*vStsHits)[ih].t_reco) - && (!isinf((*vStsHits)[ih].t_reco))) - lasttime = (*vStsHits)[ih].t_reco; + + const float& time = (*vStsHits)[ih].t_reco; + if ((lasttime < time) && (!isinf(time))) lasttime = time; + if ((starttime > time) && (time > 0)) starttime = time; + if (ist < NMvdStations) { L1StsHit& h = *(const_cast<L1StsHit*>(&((*vStsHits)[ih]))); h.t_reco = 0; @@ -1933,17 +1724,13 @@ void L1Algo::CATrackFinder() { for (int iS = 0; iS < NStations; ++iS) { - vGridTime[iS].BuildBins( - -1, 1, -0.6, 0.6, 0, lasttime, xStep, yStep, (lasttime + 1) / 35); - vGridTime[iS].StoreHits( - StsHitsUnusedStopIndex[iS] - StsHitsUnusedStartIndex[iS], - &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), - iS, - *this, - StsHitsUnusedStartIndex[iS], - &(vStsDontUsedHits_Buf[StsHitsUnusedStartIndex[iS]]), - &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), - &(RealIHit_v[StsHitsUnusedStartIndex[iS]])); + + vGridTime[iS].BuildBins(-1, 1, -0.6, 0.6, starttime, lasttime, xStep, yStep, (lasttime - starttime + 1) / 30); + + vGridTime[iS].StoreHits(StsHitsUnusedStopIndex[iS] - StsHitsUnusedStartIndex[iS], + &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), iS, *this, StsHitsUnusedStartIndex[iS], + &(vStsDontUsedHits_Buf[StsHitsUnusedStartIndex[iS]]), + &((*vStsHits)[StsHitsUnusedStartIndex[iS]]), &(RealIHit_v[StsHitsUnusedStartIndex[iS]])); } @@ -2037,30 +1824,22 @@ void L1Algo::CATrackFinder() { break; } - Pick_gather = - 3.0; /// coefficient for size of region for attach new hits to the created track - if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) - || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) + Pick_gather = 3.0; /// coefficient for size of region for attach new hits to the created track + if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) Pick_gather = 4.0; - PickNeighbour = - 1.0; // (PickNeighbour < dp/dp_error) => triplets are neighbours + PickNeighbour = 1.0; // (PickNeighbour < dp/dp_error) => triplets are neighbours // if ( (isec == kFastPrimIter) ) // PickNeighbour = 0.5; // TODO understand why works with 0.2 MaxInvMom = 1.0 / 0.5; // max considered q/p - if ((isec == kAllPrimJumpIter) || (isec == kAllSecIter) - || (isec == kAllSecJumpIter)) - MaxInvMom = 1.0 / 0.1; - if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) - || (isec == kAllSecEIter)) - MaxInvMom = 1. / 0.05; + if ((isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter)) MaxInvMom = 1.0 / 0.1; + if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter)) MaxInvMom = 1. / 0.05; MaxSlope = 1.1; if ( // (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || - (isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec == kAllSecJumpIter)) + (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) MaxSlope = 1.5; // define the target @@ -2069,21 +1848,16 @@ void L1Algo::CATrackFinder() { targZ = 0; // suppose, what target will be at (0,0,0) float SigmaTargetX = 0, SigmaTargetY = 0; // target constraint [cm] - if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) - || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter) - || (isec == kAllPrimEIter) - || (isec == kAllPrimJumpIter)) { // target + if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter) + || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter)) { // target targB = vtxFieldValue; - if ((isec == kFastPrimIter) || (isec == kAllPrimIter) - || (isec == kAllPrimEIter)) + if ((isec == kFastPrimIter) || (isec == kAllPrimIter) || (isec == kAllPrimEIter)) SigmaTargetX = SigmaTargetY = 1; // target else SigmaTargetX = SigmaTargetY = 5; } - if ( - (isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec - == kAllSecJumpIter)) { //use outer radius of the 1st station as a constraint + if ((isec == kAllSecIter) || (isec == kAllSecEIter) + || (isec == kAllSecJumpIter)) { //use outer radius of the 1st station as a constraint L1Station& st = vStations[0]; SigmaTargetX = SigmaTargetY = 10; //st.Rmax[0]; targZ = 0.; //-1.; @@ -2098,24 +1872,19 @@ void L1Algo::CATrackFinder() { /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy /// If sort by y then it is max diff between same station's modules (~0.4cm) MaxDZ = 0; - if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) - || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) + if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) MaxDZ = 0.1; - if (NStations > MaxNStations) - cout << " CATrackFinder: Error: Too many Stations" << endl; + if (NStations > MaxNStations) cout << " CATrackFinder: Error: Too many Stations" << endl; } #ifndef L1_NO_ASSERT for (int istal = NStations - 1; istal >= 0; istal--) { - L1_ASSERT(StsHitsUnusedStopIndex[istal] - >= StsHitsUnusedStartIndex[istal], - StsHitsUnusedStopIndex[istal] - << " >= " << StsHitsUnusedStartIndex[istal]); + L1_ASSERT(StsHitsUnusedStopIndex[istal] >= StsHitsUnusedStartIndex[istal], + StsHitsUnusedStopIndex[istal] << " >= " << StsHitsUnusedStartIndex[istal]); L1_ASSERT(StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(), - StsHitsUnusedStopIndex[istal] << " <= " - << (*vStsHitsUnused).size()); + StsHitsUnusedStopIndex[istal] << " <= " << (*vStsHitsUnused).size()); } #endif // L1_NO_ASSERT } @@ -2126,10 +1895,8 @@ void L1Algo::CATrackFinder() { portionStopIndex[NStations - 1] = 0; unsigned int ip = 0; //index of curent portion - for (int istal = NStations - 2; istal >= FIRSTCASTATION; - istal--) { //start downstream chambers - int NHits_l = - StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]; + for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) { //start downstream chambers + int NHits_l = StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]; int NHits_l_P = NHits_l / Portion; @@ -2194,45 +1961,28 @@ void L1Algo::CATrackFinder() { L1FieldRegion fldG_1[vPortion]; THitI hitslG_1[Portion]; - vector<THitI> - hitsm_2; /// middle hits indexed by number of doublets in portion(i2) - vector<THitI> - i1_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2) + vector<THitI> hitsm_2; /// middle hits indexed by number of doublets in portion(i2) + vector<THitI> i1_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2) - vector<THitI> - hitsmG_2; /// middle hits indexed by number of doublets in portion(i2) - vector<THitI> - i1G_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2) - vector<bool> lmDuplets - [MaxNStations]; // is exist a doublet started from indexed by left hit - vector<bool> lmDupletsG - [MaxNStations]; // is exist a doublet started from indexed by left hit + vector<THitI> hitsmG_2; /// middle hits indexed by number of doublets in portion(i2) + vector<THitI> i1G_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2) + vector<bool> lmDuplets[MaxNStations]; // is exist a doublet started from indexed by left hit + vector<bool> lmDupletsG[MaxNStations]; // is exist a doublet started from indexed by left hit hitsm_2.reserve(3500); i1_2.reserve(3500); hitsmG_2.reserve(800); i1G_2.reserve(800); - for (int istal = NStations - 2; istal >= FIRSTCASTATION; - istal--) // //start downstream chambers + for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) // //start downstream chambers { #ifdef _OPENMP -#pragma omp parallel for firstprivate(T_1, \ - fld_1, \ - hitsl_1, \ - hitsm_2, \ - i1_2, \ - TG_1, \ - fldG_1, \ - hitslG_1, \ - hitsmG_2, \ +#pragma omp parallel for firstprivate(T_1, fld_1, hitsl_1, hitsm_2, i1_2, TG_1, fldG_1, hitslG_1, hitsmG_2, \ i1G_2) //schedule(dynamic, 2) #endif - for (Tindex ip = portionStopIndex[istal + 1]; - ip < portionStopIndex[istal]; - ++ip) { + for (Tindex ip = portionStopIndex[istal + 1]; ip < portionStopIndex[istal]; ++ip) { Tindex n_2; /// number of doublets in portion int NHitsSta = StsHitsStopIndex[istal] - StsHitsStartIndex[istal]; lmDuplets[istal].resize(NHitsSta); @@ -2242,95 +1992,52 @@ void L1Algo::CATrackFinder() { i1_2.clear(); - DupletsStaPort(istal, - istal + 1, - ip, - n_g1, - portionStopIndex, + DupletsStaPort(istal, istal + 1, ip, n_g1, portionStopIndex, // output - T_1, - fld_1, - hitsl_1, + T_1, fld_1, hitsl_1, lmDuplets[istal], - n_2, - i1_2, - hitsm_2); + n_2, i1_2, hitsm_2); Tindex nstaltriplets = 0; TripletsStaPort( // input - istal, - istal + 1, - istal + 2, - nstaltriplets, - T_1, - fld_1, - hitsl_1, - - n_2, - i1_2, - hitsm_2, + istal, istal + 1, istal + 2, nstaltriplets, T_1, fld_1, hitsl_1, + + n_2, i1_2, hitsm_2, lmDuplets[istal + 1] // output ); - if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) - || (isec == kAllSecJumpIter)) { + if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter)) { Tindex nG_2; hitsmG_2.clear(); i1G_2.clear(); DupletsStaPort( // input - istal, - istal + 2, - ip, - n_g1, - portionStopIndex, + istal, istal + 2, ip, n_g1, portionStopIndex, // output - TG_1, - fldG_1, - hitslG_1, + TG_1, fldG_1, hitslG_1, lmDupletsG[istal], - nG_2, - i1G_2, - hitsmG_2); + nG_2, i1G_2, hitsmG_2); TripletsStaPort( // input - istal, - istal + 1, - istal + 3, - nstaltriplets, - T_1, - fld_1, - hitsl_1, - - n_2, - i1_2, - hitsm_2, - lmDupletsG[istal + 1]); + istal, istal + 1, istal + 3, nstaltriplets, T_1, fld_1, hitsl_1, + + n_2, i1_2, hitsm_2, lmDupletsG[istal + 1]); TripletsStaPort( // input - istal, - istal + 2, - istal + 3, - nstaltriplets, - TG_1, - fldG_1, - hitslG_1, - - nG_2, - i1G_2, - hitsmG_2, - lmDuplets[istal + 2] + istal, istal + 2, istal + 3, nstaltriplets, TG_1, fldG_1, hitslG_1, + + nG_2, i1G_2, hitsmG_2, lmDuplets[istal + 2] ); } @@ -2362,11 +2069,9 @@ void L1Algo::CATrackFinder() { // cout<<"---- Collect track candidates. ----"<<endl; // #endif - int min_level = - 0; // min level for start triplet. So min track length = min_level+3. + int min_level = 0; // min level for start triplet. So min track length = min_level+3. // if (isec == kFastPrimJumpIter) min_level = 1; - if ((isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec == kAllSecJumpIter)) + if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) min_level = 1; // only the long low momentum tracks #ifdef TRACKS_FROM_TRIPLETS @@ -2390,8 +2095,9 @@ void L1Algo::CATrackFinder() { unsigned char curr_L = 1; int ndf = 1; - vStripToTrack.assign(StsHitsStopIndex[NStations - 1], -1); - vStripToTrackB.assign(StsHitsStopIndex[NStations - 1], -1); + vStripToTrack.assign(vStripToTrack.size(), -1); + vStripToTrackB.assign(vStripToTrackB.size(), -1); + // collect consequtive: the longest tracks, shorter, more shorter and so on for (int ilev = NStations - 3; ilev >= min_level; ilev--) { @@ -2401,8 +2107,7 @@ void L1Algo::CATrackFinder() { // how many levels to check int nlevel = (NStations - 2) - ilev + 1; - const unsigned char min_best_l = - (ilev > min_level) ? ilev + 2 : min_level + 3; // loose maximum + const unsigned char min_best_l = (ilev > min_level) ? ilev + 2 : min_level + 3; // loose maximum for (int i = 0; i < fNThreads; ++i) numberCandidateThread[i] = 0; @@ -2410,13 +2115,7 @@ void L1Algo::CATrackFinder() { for (int istaF = FIRSTCASTATION; istaF <= NStations - 3 - ilev; ++istaF) { #ifdef _OPENMP -#pragma omp parallel for firstprivate(curr_tr, \ - new_tr, \ - best_tr, \ - curr_chi2, \ - best_chi2, \ - best_L, \ - curr_L, \ +#pragma omp parallel for firstprivate(curr_tr, new_tr, best_tr, curr_chi2, best_chi2, best_L, curr_L, \ ndf) // schedule(dynamic, 10) #endif for (Tindex ip = 0; ip < fNThreads; ++ip) { @@ -2429,38 +2128,30 @@ void L1Algo::CATrackFinder() { #endif L1Triplet& first_trip = (TripletsLocal1[istaF][ip][itrip]); - if (GetFUsed( - (*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f] - | (*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].b])) + if (GetFUsed((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f] + | (*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].b])) continue; // ghost supression !!! #ifndef FIND_GAPED_TRACKS - if (/*(isec == kFastPrimIter) ||*/ (isec == kAllPrimIter) - || (isec == kAllPrimEIter) || (isec == kAllSecIter) - || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) { + if (/*(isec == kFastPrimIter) ||*/ (isec == kAllPrimIter) || (isec == kAllPrimEIter) + || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) { #else - if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) - || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter) - || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) - || (isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec == kAllSecJumpIter)) { + if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) + || (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) + || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) { #endif #ifdef TRACKS_FROM_TRIPLETS if (isec != TRACKS_FROM_TRIPLETS_ITERATION) #endif { // ghost supression !!! - if (isec != kFastPrimIter && isec != kAllPrimIter - && isec != kAllPrimEIter && isec != kAllSecEIter) + if (isec != kFastPrimIter && isec != kAllPrimIter && isec != kAllPrimEIter && isec != kAllSecEIter) if (first_trip.GetLevel() == 0) continue; // ghost suppression // find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet if (!fmCBMmode) - if ((ilev == 0) - && (GetFStation(( - *vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) - != 0)) + if ((ilev == 0) && (GetFStation((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0)) continue; // ghost supression // collect only MAPS tracks-triplets CHECK!!! } if (first_trip.GetLevel() < ilev) @@ -2485,23 +2176,13 @@ void L1Algo::CATrackFinder() { best_chi2 = curr_chi2; best_L = curr_L; - CAFindTrack( - istaF, - best_tr, - best_L, - best_chi2, - &first_trip, - (curr_tr), - curr_L, - curr_chi2, - min_best_l, - new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best + CAFindTrack(istaF, best_tr, best_L, best_chi2, &first_trip, (curr_tr), curr_L, curr_chi2, min_best_l, + new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best // if ( best_L < min_best_l ) continue; if (best_L < ilev + 2) continue; // lose maximum one hit - if (best_L < min_level + 3) - continue; // should find all hits for min_level + if (best_L < min_level + 3) continue; // should find all hits for min_level ndf = best_L * 2 - 5; best_chi2 = best_chi2 / ndf; //normalize @@ -2510,50 +2191,47 @@ void L1Algo::CATrackFinder() { if (fGhostSuppression) { if (best_L == 3) { // if( isec == kAllSecIter ) continue; // too /*short*/ secondary track - if (((isec == kAllSecIter) || (isec == kAllSecEIter) - || (isec == kAllSecJumpIter)) - && (istaF != 0)) + if (((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) && (istaF != 0)) continue; // too /*short*/ non-MAPS track - if ((isec != kAllSecIter) && (isec != kAllSecEIter) - && (isec != kAllSecJumpIter) && (best_chi2 > 5.0)) + if ((isec != kAllSecIter) && (isec != kAllSecEIter) && (isec != kAllSecJumpIter) && (best_chi2 > 5.0)) continue; } } #endif best_tr.Set(istaF, best_L, best_chi2, first_trip.GetQpOrig()); - L1Branch& tr = - CandidatesTrack[thread_num][numberCandidateThread[thread_num]]; - tr = best_tr; - tr.CandIndex = - numberCandidateThread[thread_num] + thread_num * 100000; + L1Branch& tr = CandidatesTrack[thread_num][numberCandidateThread[thread_num]]; + tr = best_tr; + tr.CandIndex = numberCandidateThread[thread_num] + thread_num * 10000000; + bool check = 1; - for (vector<THitI>::iterator phitIt = - tr.StsHits.begin(); /// used strips are marked - phitIt != tr.StsHits.begin() + tr.NHits; - ++phitIt) { + for (vector<THitI>::iterator phitIt = tr.StsHits.begin(); /// used strips are marked + phitIt != tr.StsHits.begin() + tr.NHits; ++phitIt) { const L1StsHit& h = (*vStsHits)[*phitIt]; #ifdef _OPENMP omp_set_lock(&hitToBestTrackB[h.b]); #endif + int& strip1 = (vStripToTrackB)[h.b]; if (strip1 != -1) { - int thread = strip1 / 100000; - int num = (strip1 - thread * 100000); + int thread = strip1 / 10000000; + int num = (strip1 - thread * 10000000); if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])) { CandidatesTrack[thread][num].CandIndex = -1; strip1 = tr.CandIndex; - } else { + } + else { check = 0; #ifdef _OPENMP omp_unset_lock(&hitToBestTrackB[h.b]); #endif break; } - } else + } + else strip1 = tr.CandIndex; #ifdef _OPENMP omp_unset_lock(&hitToBestTrackB[h.b]); @@ -2565,21 +2243,22 @@ void L1Algo::CATrackFinder() { #endif int& strip2 = (vStripToTrack)[h.f]; if (strip2 != -1) { - int thread = strip2 / 100000; - int num = (strip2 - thread * 100000); + int thread = strip2 / 10000000; + int num = (strip2 - thread * 10000000); if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])) { CandidatesTrack[thread][num].CandIndex = -1; strip2 = tr.CandIndex; - - } else { + } + else { check = 0; #ifdef _OPENMP omp_unset_lock(&hitToBestTrackF[h.f]); #endif break; } - } else + } + else strip2 = tr.CandIndex; #ifdef _OPENMP omp_unset_lock(&hitToBestTrackF[h.f]); @@ -2588,6 +2267,7 @@ void L1Algo::CATrackFinder() { } if (check) numberCandidateThread[thread_num]++; + } // itrip } } @@ -2605,20 +2285,16 @@ void L1Algo::CATrackFinder() { #ifdef _OPENMP #pragma omp parallel for schedule(dynamic, 10) firstprivate(t) #endif - for (Tindex iCandidate = 0; iCandidate < numberCandidateThread[i]; - ++iCandidate) { + for (Tindex iCandidate = 0; iCandidate < numberCandidateThread[i]; ++iCandidate) { L1Branch& tr = CandidatesTrack[i][iCandidate]; bool check = 1; if (CandidatesTrack[i][iCandidate].CandIndex != -1) { - for (vector<THitI>::iterator phIt = - tr.StsHits.begin(); /// used strips are marked - phIt != tr.StsHits.begin() + tr.NHits; - ++phIt) { + for (vector<THitI>::iterator phIt = tr.StsHits.begin(); /// used strips are marked + phIt != tr.StsHits.begin() + tr.NHits; ++phIt) { const L1StsHit& h = (((*vStsHits))[*phIt]); - if (((vStripToTrackB)[h.b] != tr.CandIndex) - || ((vStripToTrack)[h.f] != tr.CandIndex)) { + if (((vStripToTrackB)[h.b] != tr.CandIndex) || ((vStripToTrack)[h.f] != tr.CandIndex)) { check = 0; break; } @@ -2644,10 +2320,8 @@ void L1Algo::CATrackFinder() { #endif - for (vector<THitI>::iterator phIt = - tr.StsHits.begin(); /// used strips are marked - phIt != tr.StsHits.begin() + tr.NHits; - ++phIt) { + for (vector<THitI>::iterator phIt = tr.StsHits.begin(); /// used strips are marked + phIt != tr.StsHits.begin() + tr.NHits; ++phIt) { L1StsHit& h = *(const_cast<L1StsHit*>(&(((*vStsHits))[*phIt]))); @@ -2662,17 +2336,14 @@ void L1Algo::CATrackFinder() { const L1StsHit& hit = (*vStsHits)[*phIt]; - L1HitPoint tempPoint = CreateHitPoint( - hit, 0); //TODO take number of station from hit + L1HitPoint tempPoint = CreateHitPoint(hit, 0); //TODO take number of station from hit float xcoor, ycoor = 0; L1Station stah = vStations[0]; StripsToCoor(tempPoint.U(), tempPoint.V(), xcoor, ycoor, stah); float zcoor = tempPoint.Z(); - float timeFlight = - sqrt(xcoor * xcoor + ycoor * ycoor + zcoor * zcoor) - / 30.f; // c = 30[cm/ns] + float timeFlight = sqrt(xcoor * xcoor + ycoor * ycoor + zcoor * zcoor) / 30.f; // c = 30[cm/ns] sumTime += (hit.t_reco - timeFlight); } @@ -2733,21 +2404,12 @@ void L1Algo::CATrackFinder() { int NHitsOnStation = 0; for (int ista = 0; ista < NStations; ++ista) { - int Nelements = - StsHitsUnusedStopIndex[ista] - StsHitsUnusedStartIndex[ista]; + int Nelements = StsHitsUnusedStopIndex[ista] - StsHitsUnusedStartIndex[ista]; int NHitsOnStationTmp = NHitsOnStation; vGridTime[ista].UpdateIterGrid( - Nelements, - &((*vStsHitsUnused)[StsHitsUnusedStartIndex[ista]]), - RealIHitPBuf, - &((*RealIHitP)[StsHitsUnusedStartIndex[ista]]), - vStsHitsUnused_buf, - vStsHitPointsUnused_buf, - &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), - NHitsOnStation, - ista, - *this, - vSFlag); + Nelements, &((*vStsHitsUnused)[StsHitsUnusedStartIndex[ista]]), RealIHitPBuf, + &((*RealIHitP)[StsHitsUnusedStartIndex[ista]]), vStsHitsUnused_buf, vStsHitPointsUnused_buf, + &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), NHitsOnStation, ista, *this, vSFlag); StsHitsUnusedStartIndex[ista] = NHitsOnStationTmp; StsHitsUnusedStopIndex[ista] = NHitsOnStation; } @@ -2901,32 +2563,22 @@ void L1Algo::CATrackFinder() { * * ****************************************************************/ -inline void L1Algo::CAFindTrack(int ista, - L1Branch& best_tr, - unsigned char& best_L, - fscal& best_chi2, - const L1Triplet* curr_trip, - L1Branch& curr_tr, - unsigned char& curr_L, - fscal& curr_chi2, - unsigned char min_best_l, - L1Branch* new_tr) +inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fscal& best_chi2, + const L1Triplet* curr_trip, L1Branch& curr_tr, unsigned char& curr_L, fscal& curr_chi2, + unsigned char min_best_l, L1Branch* new_tr) /// recursive search for tracks /// input: @ista - station index, @&best_tr - best track for the privious call, @&best_L - /// output: @&NCalls - number of function calls { if (curr_trip->GetLevel() == 0) // the end of the track -> check and store { - - // -- finish with current track // add rest of hits const THitI& ihitm = curr_trip->GetMHit(); const THitI& ihitr = curr_trip->GetRHit(); - if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitm].f] - | (*vSFlag)[(*vStsHitsUnused)[ihitm].b])) { + if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitm].f] | (*vSFlag)[(*vStsHitsUnused)[ihitm].b])) { // curr_tr.StsHits.push_back((*RealIHitP)[ihitm]); @@ -2937,8 +2589,7 @@ inline void L1Algo::CAFindTrack(int ista, curr_L++; } - if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitr].f] - | (*vSFlag)[(*vStsHitsUnused)[ihitr].b])) { + if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitr].f] | (*vSFlag)[(*vStsHitsUnused)[ihitr].b])) { //curr_tr.StsHits.push_back((*RealIHitP)[ihitr]); curr_tr.StsHits[curr_tr.NHits] = ((*RealIHitP)[ihitr]); @@ -2974,7 +2625,8 @@ inline void L1Algo::CAFindTrack(int ista, } return; - } else //MEANS level ! = 0 + } + else //MEANS level ! = 0 { unsigned int Station = 0; unsigned int Thread = 0; @@ -3030,17 +2682,16 @@ inline void L1Algo::CAFindTrack(int ista, if (GetFUsed((*vSFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].f] - | (*vSFlag)[(*vStsHitsUnused)[new_trip.GetLHit()] - .b])) { //hits are used + | (*vSFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].b])) { //hits are used // no used hits allowed -> compare and store track - if ((curr_L > best_L) - || ((curr_L == best_L) && (curr_chi2 < best_chi2))) { + if ((curr_L > best_L) || ((curr_L == best_L) && (curr_chi2 < best_chi2))) { best_tr = curr_tr; best_chi2 = curr_chi2; best_L = curr_L; } - } else { // if hits are used add new triplet to the current track + } + else { // if hits are used add new triplet to the current track new_tr[ista] = curr_tr; @@ -3050,8 +2701,7 @@ inline void L1Algo::CAFindTrack(int ista, fscal new_chi2 = curr_chi2; // add new hit - new_tr[ista].StsHits[new_tr[ista].NHits] = - ((*RealIHitP)[new_trip.GetLHit()]); + new_tr[ista].StsHits[new_tr[ista].NHits] = ((*RealIHitP)[new_trip.GetLHit()]); new_tr[ista].NHits++; new_L += 1; dqp = dqp / Cqp * 5.; // CHECKME: understand 5, why no sqrt(5)? @@ -3062,7 +2712,8 @@ inline void L1Algo::CAFindTrack(int ista, if (fGlobal || fmCBMmode) { new_chi2 += dtx * dtx; new_chi2 += dty * dty; - } else + } + else new_chi2 += dqp * dqp; if (new_chi2 > TRACK_CHI2_CUT * new_L) continue; @@ -3070,23 +2721,15 @@ inline void L1Algo::CAFindTrack(int ista, const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta(); - CAFindTrack(new_ista, - best_tr, - best_L, - best_chi2, - &new_trip, - new_tr[ista], - new_L, - new_chi2, - min_best_l, - new_tr); + CAFindTrack(new_ista, best_tr, best_L, best_chi2, &new_trip, new_tr[ista], new_L, new_chi2, min_best_l, new_tr); } // add triplet to track } // for neighbours } // level = 0 } #ifdef DRAW -void L1Algo::DrawRecoTracksTime(const vector<CbmL1Track>& tracks) { +void L1Algo::DrawRecoTracksTime(const vector<CbmL1Track>& tracks) +{ draw->DrawRecoTracksTime(tracks); draw->SaveCanvas(" "); } diff --git a/reco/L1/L1Algo/L1Portion.h b/reco/L1/L1Algo/L1Portion.h index a4e76095c282a4c05f8bf744e2df212783c62823..4f60a6332588dc282d43ae7d15abc52542a41d20 100644 --- a/reco/L1/L1Algo/L1Portion.h +++ b/reco/L1/L1Algo/L1Portion.h @@ -2,6 +2,9 @@ #define L1Portion_H #include <vector> + +#include "L1TrackPar.h" + using std::vector; // class for portions organization diff --git a/reco/L1/L1Algo/L1StsHit.h b/reco/L1/L1Algo/L1StsHit.h index 082eea47afa36395ccc4db53ad12912abacf7df0..f0f78baaf7d4c8e269a0594282a39a142e7d0b0f 100644 --- a/reco/L1/L1Algo/L1StsHit.h +++ b/reco/L1/L1Algo/L1StsHit.h @@ -2,8 +2,8 @@ #define _L1StsHit_h_ //struct L1Branch; -typedef unsigned /*short*/ int THitI; // strip index type -typedef unsigned short int TZPosI; // strip z-coor index type +typedef unsigned /*short*/ int THitI; // strip index type +typedef unsigned short int TZPosI; // strip z-coor index type typedef unsigned /*short*/ int TStripI; // strip index type class L1StsHit { @@ -13,13 +13,14 @@ public: float du, dv; float t_reco; float t_er; + int ID; // int ista; //int track, n; float u = 0, v = 0; TZPosI iz; // index of z coor. in L1Algo::vStsZPos - L1StsHit() : f(0), b(0), du(0.), dv(0.), t_reco(0.f), t_er(0.), iz(0) {} + L1StsHit() : f(0), b(0), du(0.), dv(0.), t_reco(0.f), t_er(0.), ID(0.), iz(0) {} }; #endif diff --git a/reco/L1/L1AlgoInputData.h b/reco/L1/L1AlgoInputData.h index 86739300f5c642c85b8867feaacc2a9500f0dce0..c2d3551d22909f5fdb9d520781c95232c980942e 100644 --- a/reco/L1/L1AlgoInputData.h +++ b/reco/L1/L1AlgoInputData.h @@ -2,23 +2,20 @@ #define _L1AlgoInputData_h #include "CbmL1Def.h" -#include "L1StsHit.h" #include <fstream> #include <iostream> #include <vector> +#include "L1StsHit.h" + using std::istream; using std::vector; class L1AlgoInputData { public: - L1AlgoInputData() - : vStsHits() - , NStsStrips(0) - , vStsZPos() - , vSFlag() + L1AlgoInputData() : vStsHits(), NStsStrips(0), vStsZPos(), vSFlag() // MaxNStations(12) { @@ -37,18 +34,12 @@ public: const THitI* GetStsHitsStopIndex() const { return StsHitsStopIndex; } - bool ReadHitsFromFile(const char work_dir[100], - const int maxNEvent, - const int iVerbose); + bool ReadHitsFromFile(const char work_dir[100], const int maxNEvent, const int iVerbose); // void PrintHits(); /// redefine new\delete for use alignment memmory - void* operator new(size_t size, void* ptr) { - return ::operator new(size, ptr); - } - void* operator new[](size_t size, void* ptr) { - return ::operator new(size, ptr); - } + void* operator new(size_t size, void* ptr) { return ::operator new(size, ptr); } + void* operator new[](size_t size, void* ptr) { return ::operator new(size, ptr); } void* operator new(size_t size) { return _mm_malloc(size, 16); } void* operator new[](size_t size) { return _mm_malloc(size, 16); } void operator delete(void* ptr, size_t) { _mm_free(ptr); } @@ -58,12 +49,13 @@ public: const L1AlgoInputData& operator=(const L1AlgoInputData& a); - void Clear() { + void Clear() + { vStsHits.clear(); NStsStrips = 0; vStsZPos.clear(); - vSFlag.clear(); + //vSFlag.clear(); { for (int i = 0; i < MaxNStations + 1; ++i) @@ -79,14 +71,15 @@ public: static istream& eatwhite(istream& is); // skip spaces /// read data from data_algo.txt // data - enum { MaxNStations = 25 }; - vector<L1StsHit> - vStsHits; // hits as a combination of front-, backstrips and z-position - int NStsStrips; // Number of strips in sts - vector<fscal> vStsZPos; // all possible z-positions of hits - - L1Vector<unsigned char> - vSFlag; // information of hits station & using hits in tracks; + enum + { + MaxNStations = 25 + }; + vector<L1StsHit> vStsHits; // hits as a combination of front-, backstrips and z-position + int NStsStrips; // Number of strips in sts + vector<fscal> vStsZPos; // all possible z-positions of hits + + L1Vector<unsigned char> vSFlag; // information of hits station & using hits in tracks; THitI StsHitsStartIndex[MaxNStations + 1], StsHitsStopIndex[MaxNStations + 1]; // station-bounders in vStsHits array