/** @file CbmTofSimpClusterizer.cxx ** @author Pierre-Alain Loizeau ** @date 23.08.2013 **/ #include "CbmTofSimpClusterizer.h" // TOF Classes and includes #include "CbmTofPoint.h" // in cbmdata/tof #include "CbmTofDigi.h" // in cbmdata/tof #include "CbmTofDigiExp.h" // in cbmdata/tof #include "CbmTofHit.h" // in cbmdata/tof #include "CbmTofGeoHandler.h" // in tof/TofTools #include "CbmTofDetectorId_v12b.h" // in cbmdata/tof #include "CbmTofCell.h" // in tof/TofData #include "CbmTofDigiPar.h" // in tof/TofParam #include "CbmTofDigiBdfPar.h" // in tof/TofParam #include "CbmTofAddress.h" // in cbmdata/tof // CBMroot classes and includes #include "CbmMCTrack.h" // FAIR classes and includes #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairLogger.h" // ROOT Classes and includes #include "TClonesArray.h" #include "TMath.h" #include "TLine.h" #include "TRandom3.h" #include "TF2.h" #include "TVector3.h" #include "TH1.h" #include "TH2.h" #include "TDirectory.h" #include "TROOT.h" // C++ Classes and includes /************************************************************************************/ CbmTofSimpClusterizer::CbmTofSimpClusterizer(): FairTask("CbmTofSimpClusterizer"), fVerbose(1), fGeoHandler(new CbmTofGeoHandler()), fDigiPar(NULL), fChannelInfo(NULL), fTofPointsColl(NULL), fMcTracksColl(NULL), fbWriteHitsInOut(kTRUE), fTofDigisColl(NULL) { } CbmTofSimpClusterizer::CbmTofSimpClusterizer(const char *name, Int_t verbose, Bool_t writeDataInOut): FairTask(TString(name),verbose), fVerbose(verbose), fGeoHandler(new CbmTofGeoHandler()), fDigiPar(NULL), fChannelInfo(NULL), fTofPointsColl(NULL), fMcTracksColl(NULL), fbWriteHitsInOut(writeDataInOut), fTofDigisColl(NULL) { } CbmTofSimpClusterizer::~CbmTofSimpClusterizer() { if( fGeoHandler ) delete fGeoHandler; // DeleteHistos(); // <-- if needed ? } /************************************************************************************/ // FairTasks inherited functions InitStatus CbmTofSimpClusterizer::Init() { if( kFALSE == RegisterInputs() ) return kFATAL; if( kFALSE == RegisterOutputs() ) return kFATAL; if( kFALSE == InitParameters() ) return kFATAL; if( kFALSE == LoadGeometry() ) return kFATAL; if( kFALSE == CreateHistos() ) return kFATAL; return kSUCCESS; } void CbmTofSimpClusterizer::SetParContainers() { LOG(INFO)<<" CbmTofSimpClusterizer => Get the digi parameters for tof"<GetRuntimeDb(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); } void CbmTofSimpClusterizer::Exec(Option_t * option) { fTofHitsColl->Clear("C"); fiNbHits = 0; LOG(DEBUG)<<" CbmTofSimpClusterizer => New event"<GetObject("TofPoint"); if( NULL == fTofPointsColl) { LOG(ERROR)<<"CbmTofSimpClusterizer::RegisterInputs => Could not get the TofPoint TClonesArray!!!"<GetObject("MCTrack"); if( NULL == fMcTracksColl) { LOG(ERROR)<<"CbmTofSimpClusterizer::RegisterInputs => Could not get the MCTrack TClonesArray!!!"<GetObject("TofDigi"); if( NULL == fTofDigisColl) { LOG(ERROR)<<"CbmTofSimpClusterizer::RegisterInputs => Could not get the TofDigi TClonesArray!!!"<Register( "TofHit","Tof", fTofHitsColl, fbWriteHitsInOut); return kTRUE; } Bool_t CbmTofSimpClusterizer::InitParameters() { // Initialize the TOF GeoHandler Bool_t isSimulation=kFALSE; Int_t iGeoVersion = fGeoHandler->Init(isSimulation); if( k12b > iGeoVersion ) { LOG(ERROR)<<"CbmTofSimpClusterizer::InitParameters => Only compatible with geometries after v12b !!!" <GetNbSmTypes(); if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { fStorDigiExp.resize( iNbSmTypes ); fviClusterSize.resize( iNbSmTypes ); fviTrkMul.resize( iNbSmTypes ); fvdX.resize( iNbSmTypes ); fvdY.resize( iNbSmTypes ); fvdDifX.resize( iNbSmTypes ); fvdDifY.resize( iNbSmTypes ); fvdDifCh.resize( iNbSmTypes ); for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); // fStorDigiExp[iSmType].resize( iNbSm ); fStorDigiExp[iSmType].resize( iNbSm*iNbRpc ); fviClusterSize[iSmType].resize( iNbRpc ); fviTrkMul[iSmType].resize( iNbRpc ); fvdX[iSmType].resize( iNbRpc ); fvdY[iSmType].resize( iNbRpc ); fvdDifX[iSmType].resize( iNbRpc ); fvdDifY[iSmType].resize( iNbRpc ); fvdDifCh[iSmType].resize( iNbRpc ); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) { // fStorDigiExp[iSmType][iSm].resize( iNbRpc ); for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { Int_t iNbChan = fDigiBdfPar->GetNbChan( iSmType, iRpc ); // fStorDigiExp[iSmType][iSm][iRpc].resize( iNbChan ); fStorDigiExp[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); } // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) else { fStorDigi.resize( iNbSmTypes ); fviClusterSize.resize( iNbSmTypes ); fviTrkMul.resize( iNbSmTypes ); fvdX.resize( iNbSmTypes ); fvdY.resize( iNbSmTypes ); fvdDifX.resize( iNbSmTypes ); fvdDifY.resize( iNbSmTypes ); fvdDifCh.resize( iNbSmTypes ); for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); // fStorDigi[iSmType].resize( iNbSm ); fStorDigi[iSmType].resize( iNbSm*iNbRpc ); fviClusterSize[iSmType].resize( iNbRpc ); fviTrkMul[iSmType].resize( iNbRpc ); fvdX[iSmType].resize( iNbRpc ); fvdY[iSmType].resize( iNbRpc ); fvdDifX[iSmType].resize( iNbRpc ); fvdDifY[iSmType].resize( iNbRpc ); fvdDifCh[iSmType].resize( iNbRpc ); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) { // fStorDigi[iSmType][iSm].resize( iNbRpc ); for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { Int_t iNbChan = fDigiBdfPar->GetNbChan( iSmType, iRpc ); fStorDigi[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); // fStorDigi[iSmType][iSm][iRpc].resize( iNbChan ); } // for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) return kTRUE; } Bool_t CbmTofSimpClusterizer::DeleteGeometry() { Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes(); if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) { for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) // fStorDigiExp[iSmType][iSm][iRpc].clear(); fStorDigiExp[iSmType][iSm*iNbRpc+iRpc].clear(); // fStorDigiExp[iSmType][iSm].clear(); } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) fStorDigiExp[iSmType].clear(); } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) fStorDigiExp.clear(); } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) else { for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) { for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) // fStorDigi[iSmType][iSm][iRpc].clear(); fStorDigi[iSmType][iSm*iNbRpc+iRpc].clear(); // fStorDigi[iSmType][iSm].clear(); } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) fStorDigi[iSmType].clear(); } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) fStorDigi.clear(); } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) return kTRUE; } /************************************************************************************/ // Histogramming functions Bool_t CbmTofSimpClusterizer::CreateHistos() { TDirectory * oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager! gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager ! fhClustBuildTime = new TH1I( "TofSimpClus_ClustBuildTime", "Time needed to build clusters in each event; Time [s]", 4000, 0.0, 4.0 ); fhHitsPerTracks = new TH1I( "TofSimpClus_TofHitPerTrk", "Mean Number of TofHit per Mc Track; Nb TofHits/Nb MC Tracks []", 2000, 0.0, 20.0 ); if( kFALSE == fDigiBdfPar->ClustUseTrackId() ) fhPtsPerHit = new TH1I( "TofSimpClus_TofPtsPerHit", "Distribution of the Number of MCPoints associated to each TofHit; Nb MCPoint []", 20, 0.0, 20.0 ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { fhTimeResSingHits = new TH1I( "TofSimpClus_TofTimeResClust", "Time resolution for TofHits containing Digis from a single MC Track; t(1st Mc Point) -tTofHit [ns]", 10000, -25.0, 25.0 ); fhTimeResSingHitsB = new TH2I( "TofSimpClus_TofTimeResClustB", "Time resolution for TofHits containing Digis from a single MC Track; (1st Mc Point) -tTofHit [ns]", 5000, -25.0, 25.0, 6, 0, 6); fhTimePtVsHits = new TH2I( "TofSimpClus_TofTimePtVsHit", "Time resolution for TofHits containing Digis from a single MC Track; t(1st Mc Point) -tTofHit [ps]", 2000, 0.0, 50000.0, 2000, 0.0, 50000.0 ); } else { fhTimeResSingHits = new TH1I( "TofSimpClus_TofTimeResClust", "Time resolution for TofHits containing Digis from a single TofPoint; tMcPoint -tTofHit [ns]", 10000, -25.0, 25.0 ); fhTimeResSingHitsB = new TH2I( "TofSimpClus_TofTimeResClustB", "Time resolution for TofHits containing Digis from a single TofPoint; tMcPoint -tTofHit [ns]", 5000, -25.0, 25.0, 6, 0, 6); fhTimePtVsHits = new TH2I( "TofSimpClus_TofTimePtVsHit", "Time resolution for TofHits containing Digis from a single TofPoint; tMcPoint -tTofHit [ps]", 2000, 0.0, 50000.0, 2000, 0.0, 50000.0 ); } // else of if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) fhClusterSize = new TH1I( "TofSimpClus_ClusterSize", "Cluster Size distribution; Cluster Size [Strips]", 100, 0.5, 100.5 ); fhClusterSizeType = new TH2I( "TofSimpClus_ClusterSizeType", "Cluster Size distribution in each (SM type, Rpc) pair; Cluster Size [Strips]; 10*SM Type + Rpc Index []", 100, 0.5, 100.5, 40*fDigiBdfPar->GetNbSmTypes(), 0.0, 40*fDigiBdfPar->GetNbSmTypes() ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { fhTrackMul = new TH1I( "TofSimpClus_TrackMul", "Number of MC tracks generating the cluster; MC Tracks multiplicity []", 100, 0.5, 100.5 ); fhClusterSizeMulti = new TH2I( "TofSimpClus_ClusterSizeMulti", "Cluster Size distribution as function of Number of MC tracks generating the cluster; Cluster Size [Strips]; MC tracks mul. []", 100, 0.5, 100.5, 100, 0.5, 100.5 ); fhTrk1MulPos = new TH2D( "TofSimpClus_Trk1MulPos", "Position of Clusters with only 1 MC tracks generating the cluster; X [cm]; Y [cm]", 1500, -750, 750, 1000, -500, 500 ); fhHiTrkMulPos = new TH2D( "TofSimpClus_HiTrkMulPos", "Position of Clusters with >1 MC tracks generating the cluster; X [cm]; Y [cm]", 1500, -750, 750, 1000, -500, 500 ); fhAllTrkMulPos = new TH2D( "TofSimpClus_AllTrkMulPos", "Position of all clusters generating the cluster; X [cm]; Y [cm]", 1500, -750, 750, 1000, -500, 500 ); fhMultiTrkProbPos = new TH2D( "TofSimpClus_MultiTrkProbPos", "Probability of having a cluster with multiple tracks as function of position; X [cm]; Y [cm]; Prob. [%]", 1500, -750, 750, 1000, -500, 500 ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) fhDigSpacDifClust = new TH1I( "TofSimpClus_DigSpacDifClust", "Space difference along channel direction between Digi pairs on adjacent channels; PosCh i - Pos Ch i+1 [cm]", 5000, -10.0, 10.0 ); fhDigTimeDifClust = new TH1I( "TofSimpClus_DigTimeDifClust", "Time difference between Digi pairs on adjacent channels; 0.5*(tDigiA + tDigiA)chi - 0.5*(tDigiA + tDigiA)chi+1 [ns]", 5000, -5.0, 5.0 ); fhDigDistClust = new TH2I( "TofSimpClus_DigDistClust", "Distance between Digi pairs on adjacent channels; PosCh i - Pos Ch i+1 [cm along ch]; 0.5*(tDigiA + tDigiA)chi - 0.5*(tDigiA + tDigiA)chi+1 [ns]", 5000, -10.0, 10.0, 2000, -5.0, 5.0 ); fhClustSizeDifX = new TH2I( "TofSimpClus_ClustSizeDifX", "Position difference between Point and Hit as function of Cluster Size; Cluster Size [Strips]; dX [cm]", 100, 0.5, 100.5, 500,-50,50); fhClustSizeDifY = new TH2I( "TofSimpClus_ClustSizeDifY", "Position difference between Point and Hit as function of Cluster Size; Cluster Size [Strips]; dY [cm]", 100, 0.5, 100.5, 500,-50,50); fhChDifDifX = new TH2I( "TofSimpClus_ChDifDifX", "Position difference between Point and Hit as function of Channel dif; Ch Dif [Strips]; dX [cm]", 101, -50.5, 50.5, 500,-50,50); fhChDifDifY = new TH2I( "TofSimpClus_ChDifDifY", "Position difference between Point and Hit as function of Channel Dif; Ch Dif [Strips]; dY [cm]", 101, -50.5, 50.5, 500,-50,50); fhNbSameSide = new TH1I( "TofSimpClus_NbSameSide", "How many time per event the 2 digis on a channel were of the same side ; Counts/Event []", 500, 0.0, 500.0 ); fhNbDigiPerChan = new TH1I( "TofSimpClus_NbDigiPerChan", "Nb of Digis per channel; Nb Digis []", 100, 0.0, 100.0 ); gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager! return kTRUE; } Bool_t CbmTofSimpClusterizer::FillHistos() { fhClustBuildTime->Fill( fStop.GetSec() - fStart.GetSec() + (fStop.GetNanoSec() - fStart.GetNanoSec())/1e9 ); Int_t iNbTofHits = fTofHitsColl->GetEntries(); Int_t iNbTracks = fMcTracksColl->GetEntries(); // Trakcs Info Int_t iNbTofTracks = 0; Int_t iNbTofTracksPrim = 0; CbmMCTrack *pMcTrk; for(Int_t iTrkInd = 0; iTrkInd < iNbTracks; iTrkInd++) { pMcTrk = (CbmMCTrack*) fMcTracksColl->At( iTrkInd ); if( 0 < pMcTrk->GetNPoints(kTOF) ) { iNbTofTracks++; } if( 0 < pMcTrk->GetNPoints(kTOF) && -1 == pMcTrk->GetMotherId() ) iNbTofTracksPrim++; } // for(Int_t iTrkInd = 0; iTrkInd < nMcTracks; iTrkInd++) if( 0 < iNbTofTracks ) fhHitsPerTracks->Fill( (Double_t)iNbTofHits/(Double_t)iNbTofTracks ); CbmTofHit *pHit; for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if( kFALSE == fDigiBdfPar->ClustUseTrackId() ) fhPtsPerHit->Fill(pHit->GetFlag()); if( 1 == pHit->GetFlag() ) { // CbmTofPoint* pPt = (CbmTofPoint*)pHit->GetRefId(); // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // CbmTofPoint* pPt = (CbmTofPoint*)(pHit->GetLinks()); // Use instead the index CbmTofPoint* pPt = (CbmTofPoint*)fTofPointsColl->At( pHit->GetRefId() ); fhTimePtVsHits->Fill( pPt->GetTime(), pHit->GetTime() ); fhTimeResSingHits->Fill( pHit->GetTime() - pPt->GetTime() ); fhTimeResSingHitsB->Fill( pHit->GetTime() - pPt->GetTime(), fGeoHandler->GetSMType(pPt->GetDetectorID()) ); } // if( 1 == pHit->GetFlag() ) } // for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) for( Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++ ) for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ ) { for( Int_t iCluster = 0; iCluster < fviClusterSize[iSmType][iRpc].size(); iCluster++ ) { fhClusterSize->Fill( fviClusterSize[iSmType][iRpc][iCluster]); fhClusterSizeType->Fill( fviClusterSize[iSmType][iRpc][iCluster], 40*iSmType + iRpc ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { fhTrackMul->Fill( fviTrkMul[iSmType][iRpc][iCluster] ); fhClusterSizeMulti->Fill( fviClusterSize[iSmType][iRpc][iCluster], fviTrkMul[iSmType][iRpc][iCluster] ); if( 1 == fviTrkMul[iSmType][iRpc][iCluster] ) fhTrk1MulPos->Fill( fvdX[iSmType][iRpc][iCluster], fvdY[iSmType][iRpc][iCluster] ); if( 1 < fviTrkMul[iSmType][iRpc][iCluster] ) fhHiTrkMulPos->Fill(fvdX[iSmType][iRpc][iCluster], fvdY[iSmType][iRpc][iCluster] ); fhAllTrkMulPos->Fill(fvdX[iSmType][iRpc][iCluster], fvdY[iSmType][iRpc][iCluster] ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) if( 1 == fviTrkMul[iSmType][iRpc][iCluster] ) { fhClustSizeDifX->Fill( fviClusterSize[iSmType][iRpc][iCluster], fvdDifX[iSmType][iRpc][iCluster]); fhClustSizeDifY->Fill( fviClusterSize[iSmType][iRpc][iCluster], fvdDifY[iSmType][iRpc][iCluster]); if( 1 == fviClusterSize[iSmType][iRpc][iCluster] ) { fhChDifDifX->Fill( fvdDifCh[iSmType][iRpc][iCluster], fvdDifX[iSmType][iRpc][iCluster]); fhChDifDifY->Fill( fvdDifCh[iSmType][iRpc][iCluster], fvdDifY[iSmType][iRpc][iCluster]); } } } // for( Int_t iCluster = 0; iCluster < fviClusterSize[iSmType][iRpc].size(); iCluster++ ) fviClusterSize[iSmType][iRpc].clear(); fviTrkMul[iSmType][iRpc].clear(); fvdX[iSmType][iRpc].clear(); fvdY[iSmType][iRpc].clear(); fvdDifX[iSmType][iRpc].clear(); fvdDifY[iSmType][iRpc].clear(); fvdDifCh[iSmType][iRpc].clear(); } // for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ ) fhNbSameSide->Fill(fiNbSameSide); return kTRUE; } Bool_t CbmTofSimpClusterizer::WriteHistos() { TDirectory * oldir = gDirectory; TFile *fHist = new TFile("./tofSimpClust.hst.root","RECREATE"); fHist->cd(); fhClustBuildTime->Write(); fhHitsPerTracks->Write(); if( kFALSE == fDigiBdfPar->ClustUseTrackId() ) fhPtsPerHit->Write(); fhTimeResSingHits->Write(); fhTimeResSingHitsB->Write(); fhTimePtVsHits->Write(); fhClusterSize->Write(); fhClusterSizeType->Write(); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { fhTrackMul->Write(); fhClusterSizeMulti->Write(); fhTrk1MulPos->Write(); fhHiTrkMulPos->Write(); fhAllTrkMulPos->Write(); fhMultiTrkProbPos->Divide( fhHiTrkMulPos, fhAllTrkMulPos); fhMultiTrkProbPos->Scale( 100.0 ); fhMultiTrkProbPos->Write(); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) fhDigSpacDifClust->Write(); fhDigTimeDifClust->Write(); fhDigDistClust->Write(); fhClustSizeDifX->Write(); fhClustSizeDifY->Write(); fhChDifDifX->Write(); fhChDifDifY->Write(); fhNbSameSide->Write(); fhNbDigiPerChan->Write(); gDirectory->cd( oldir->GetPath() ); fHist->Close(); return kTRUE; } Bool_t CbmTofSimpClusterizer::DeleteHistos() { delete fhClustBuildTime; delete fhHitsPerTracks; delete fhPtsPerHit; delete fhTimeResSingHits; delete fhTimeResSingHitsB; delete fhTimePtVsHits; delete fhClusterSize; delete fhClusterSizeType; if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { delete fhTrackMul; delete fhClusterSizeMulti; delete fhTrk1MulPos; delete fhHiTrkMulPos; delete fhAllTrkMulPos; delete fhMultiTrkProbPos; } delete fhDigSpacDifClust; delete fhDigTimeDifClust; delete fhDigDistClust; delete fhClustSizeDifX; delete fhClustSizeDifY; delete fhChDifDifX; delete fhChDifDifY; delete fhNbSameSide; delete fhNbDigiPerChan; return kTRUE; } /************************************************************************************/ Bool_t CbmTofSimpClusterizer::BuildClusters() { /* * FIXME: maybe use the 2D distance (X/Y) as criteria instead of the distance long channel * direction */ Double_t dMaxTimeDist = fDigiBdfPar->GetMaxTimeDist(); Double_t dMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh(); // First Time order the Digis array // fTofDigisColl->Sort(); // Then loop over the digis array and store the Digis in separate vectors for // each RPC modules Int_t iNbTofDigi = fTofDigisColl->GetEntries(); if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { CbmTofDigiExp *pDigi; for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd ); // cout<GetType()<<" " // <GetSm()<<" " // <GetRpc()<<" " // <GetChannel()<<" " // <GetTime()<<" " // <GetTot() // <GetType()][pDigi->GetSm()][pDigi->GetRpc()][pDigi->GetChannel()].push_back(pDigi); fStorDigiExp[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(pDigi); } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ ) } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) else { CbmTofDigi *pDigi; for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { pDigi = (CbmTofDigi*) fTofDigisColl->At( iDigInd ); // fStorDigi[pDigi->GetType()][pDigi->GetSm()][pDigi->GetRpc()][pDigi->GetChannel()].push_back(pDigi); fStorDigi[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(pDigi); } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ ) } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) // Then build clusters inside each RPC module // Assume only 0 or 1 Digi per channel/side in each event // Use simplest method possible, scan direction independent: // a) Loop over channels in the RPC starting from 0 // * If strips // i) Loop over Digis to check if both ends of the channel have a Digi // ii) Reconstruct a mean channel time and a mean position // + If a Hit is currently filled & the mean position (space, time) is less than XXX from last channel position // iii) Add the mean channel time and the mean position to the ones of the hit // + else // iii) Use nb of strips in cluster to cal. the hit mean time and pos (charge/tot weighting) // iv) Save the hit // v) Start a new hit with current channel // * else (pads) // i) Loop over Digis to find if this channel fired // ii) FIXME: either scan all other channels to check for matching Digis or have more than 1 hit open Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes(); // Hit variables Double_t dWeightedTime = 0.0; Double_t dWeightedPosX = 0.0; Double_t dWeightedPosY = 0.0; Double_t dWeightedPosZ = 0.0; Double_t dWeightsSum = 0.0; std::vector< CbmTofPoint* > vPtsRef; Int_t iNbChanInHit = 0; // Last Channel Temp variables Int_t iLastChan = -1; Double_t dLastPosX = 0.0; Double_t dLastPosY = 0.0; Double_t dLastTime = 0.0; // Channel Temp variables Double_t dPosX = 0.0; Double_t dPosY = 0.0; Double_t dPosZ = 0.0; Double_t dTime = 0.0; Double_t dTotS = 0.0; fiNbSameSide = 0; if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); if( 0 == iChType ) { Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); // Don't spread clusters over RPCs!!! dWeightedTime = 0.0; dWeightedPosX = 0.0; dWeightedPosY = 0.0; dWeightedPosZ = 0.0; dWeightsSum = 0.0; iNbChanInHit = 0; vPtsRef.clear(); // For safety reinitialize everything iLastChan = -1; dLastPosX = 0.0; dLastPosY = 0.0; dLastTime = 0.0; if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) { // Horizontal strips => X comes from left right time difference for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { if( 0 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); // if( 2 == fStorDigiExp[iSmType][iSm][iRpc][iCh].size() ) if( 2 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { if( (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; continue; } // 2 Digis = both sides present CbmTofDetectorInfo xDetInfo(kTOF, iSmType, iSm, iRpc, 0, iCh + 1); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); // CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm][iRpc][iCh][0]; // CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm][iRpc][iCh][1]; CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]; CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1]; // The "Strip" time is the mean time between each end dTime = xDigiA->GetTime() + xDigiB->GetTime(); dTime /= 2.0; // Weight is the total charge => sum of both ends ToT dTotS = xDigiA->GetTot() + xDigiB->GetTot(); // Y position is just the center of the channel dPosY = fChannelInfo->GetY(); // X position form strip ends time difference dPosX = fChannelInfo->GetX(); if( 1 == xDigiA->GetSide() ) // 0 is the right side, 1 is the left side dPosX += fDigiBdfPar->GetSignalSpeed() * ( xDigiA->GetTime() - xDigiB->GetTime() )/2.0; else // 0 is the left side, 1 is the right side dPosX += fDigiBdfPar->GetSignalSpeed() * ( xDigiB->GetTime() - xDigiA->GetTime() )/2.0; // For Z always just take the one of the channel itself( in fact its gap one) dPosZ = fChannelInfo->GetZ(); // Now check if a hit/cluster is already started if( 0 < iNbChanInHit) { if( iLastChan == iCh - 1 ) { fhDigTimeDifClust->Fill( dTime - dLastTime ); fhDigSpacDifClust->Fill( dPosX - dLastPosX ); fhDigDistClust->Fill( dPosX - dLastPosX, dTime - dLastTime ); } // if( iLastChan == iCh - 1 ) // a cluster is already started => check distance in space/time // For simplicity, just check along strip direction for now // and break cluster when a not fired strip is found if( TMath::Abs( dTime - dLastTime) < dMaxTimeDist && iLastChan == iCh - 1 && TMath::Abs( dPosX - dLastPosX) < dMaxSpaceDist ) { // Add to cluster/hit dWeightedTime += dTime*dTotS; dWeightedPosX += dPosX*dTotS; dWeightedPosY += dPosY*dTotS; dWeightedPosZ += dPosZ*dTotS; dWeightsSum += dTotS; iNbChanInHit += 1; // if one of the side digi comes from a CbmTofPoint not already found // in this cluster, save its pointer Bool_t bFoundA = kFALSE; Bool_t bFoundB = kFALSE; if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() ) bFoundA = kTRUE; if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) bFoundB = kTRUE; } // for( Int else for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiA->GetLinks() ) bFoundA = kTRUE; if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiB->GetLinks() ) bFoundB = kTRUE; } // for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) // CbmTofPoint pointer for 1st digi not found => save it if( kFALSE == bFoundA ) vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // CbmTofPoint pointer for 2nd digi not found => save it if( kFALSE == bFoundB ) vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if current Digis compatible with last fired chan else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fChannelInfo->GetSizey()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); // Start a new hit dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if current Digis compatible with last fired chan } // if( 0 < iNbChanInHit) else { // first fired strip in this RPC dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if( 0 < iNbChanInHit) iLastChan = iCh; dLastPosX = dPosX; dLastPosY = dPosY; dLastTime = dTime; } // if( 2 <= fStorDigiExp[iSmType][iSm][iRpc][iCh].size() ) // fStorDigiExp[iSmType][iSm][iRpc][iCh].clear(); fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { // Vertical strips => Y comes from bottom top time difference for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { if( 0 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); // if( 2 == fStorDigiExp[iSmType][iSm][iRpc][iCh].size() ) if( 2 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { if( (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; continue; } // 2 Digis = both sides present CbmTofDetectorInfo xDetInfo(kTOF, iSmType, iSm, iRpc, 0, iCh + 1); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); // CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm][iRpc][iCh][0]; // CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm][iRpc][iCh][1]; CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]; CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1]; // The "Strip" time is the mean time between each end dTime = xDigiA->GetTime() + xDigiB->GetTime(); dTime /= 2.0; // Weight is the total charge => sum of both ends ToT dTotS = xDigiA->GetTot() + xDigiB->GetTot(); // X position is just the center of the channel dPosX = fChannelInfo->GetX(); // Y position form strip ends time difference dPosY = fChannelInfo->GetY(); if( 1 == xDigiA->GetSide() ) // 0 is the top side, 1 is the bottom side dPosY += fDigiBdfPar->GetSignalSpeed() * ( xDigiA->GetTime() - xDigiB->GetTime() )/2.0; else // 0 is the bottom side, 1 is the top side dPosY += fDigiBdfPar->GetSignalSpeed() * ( xDigiB->GetTime() - xDigiA->GetTime() )/2.0; // For Z always just take the one of the channel itself( in fact its gap one) dPosZ = fChannelInfo->GetZ(); // Now check if a hit/cluster is already started if( 0 < iNbChanInHit) { if( iLastChan == iCh - 1 ) { fhDigTimeDifClust->Fill( dTime - dLastTime ); fhDigSpacDifClust->Fill( dPosY - dLastPosY ); fhDigDistClust->Fill( dPosY - dLastPosY, dTime - dLastTime ); } // if( iLastChan == iCh - 1 ) // a cluster is already started => check distance in space/time // For simplicity, just check along strip direction for now // and break cluster when a not fired strip is found if( TMath::Abs( dTime - dLastTime) < dMaxTimeDist && iLastChan == iCh - 1 && TMath::Abs( dPosY - dLastPosY) < dMaxSpaceDist ) { // Add to cluster/hit dWeightedTime += dTime*dTotS; dWeightedPosX += dPosX*dTotS; dWeightedPosY += dPosY*dTotS; dWeightedPosZ += dPosZ*dTotS; dWeightsSum += dTotS; iNbChanInHit += 1; // if one of the side digi comes from a CbmTofPoint not already found // in this cluster, save its pointer Bool_t bFoundA = kFALSE; Bool_t bFoundB = kFALSE; if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() ) bFoundA = kTRUE; if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) bFoundB = kTRUE; } // for( Int else for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiA->GetLinks() ) bFoundA = kTRUE; if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiB->GetLinks() ) bFoundB = kTRUE; } // for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) // CbmTofPoint pointer for 1st digi not found => save it if( kFALSE == bFoundA ) vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // CbmTofPoint pointer for 2nd digi not found => save it if( kFALSE == bFoundB ) vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if current Digis compatible with last fired chan else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); // Start a new hit dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if current Digis compatible with last fired chan } // if( 0 < iNbChanInHit) else { // first fired strip in this RPC dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if( 0 < iNbChanInHit) iLastChan = iCh; dLastPosX = dPosX; dLastPosY = dPosY; dLastTime = dTime; } // if( 2 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) // fStorDigiExp[iSmType][iSm][iRpc][iCh].clear(); fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 == iChType) else { LOG(ERROR)<<"CbmTofSimpClusterizer::BuildClusters => Cluster building " <<"from digis to hits not implemented for pads, Sm type " <GetChanOrient( iSmType, iRpc ) ) { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fChannelInfo->GetSizey()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size // cout<<"a "<GetDetectorID()<GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 < iNbChanInHit) } // for each sm/rpc pair } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) else { for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); if( 0 == iChType ) { Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); // Don't spread clusters over RPCs!!! dWeightedTime = 0.0; dWeightedPosX = 0.0; dWeightedPosY = 0.0; dWeightedPosZ = 0.0; dWeightsSum = 0.0; iNbChanInHit = 0; vPtsRef.clear(); // For safety reinitialize everything iLastChan = -1; dLastPosX = 0.0; dLastPosY = 0.0; dLastTime = 0.0; if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) { // Horizontal strips => X comes from left right time difference for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { if( 0 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); // if( 2 == fStorDigi[iSmType][iSm][iRpc][iCh].size() ) if( 2 == fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { if( (fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; continue; } CbmTofDetectorInfo xDetInfo(kTOF, iSmType, iSm, iRpc, 0, iCh + 1); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); // CbmTofDigi * xDigiA = fStorDigi[iSmType][iSm][iRpc][iCh][0]; // CbmTofDigi * xDigiB = fStorDigi[iSmType][iSm][iRpc][iCh][1]; CbmTofDigi * xDigiA = fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][0]; CbmTofDigi * xDigiB = fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][1]; // The "Strip" time is the mean time between each end dTime = xDigiA->GetTime() + xDigiB->GetTime(); dTime /= 2.0; // Weight is the total charge => sum of both ends ToT dTotS = xDigiA->GetTot() + xDigiB->GetTot(); // Y position is just the center of the channel dPosY = fChannelInfo->GetY(); // X position form strip ends time difference dPosX = fChannelInfo->GetX(); if( 1 == xDigiA->GetSide() ) // 0 is the right side, 1 is the left side dPosX += fDigiBdfPar->GetSignalSpeed() * ( xDigiA->GetTime() - xDigiB->GetTime() )/2.0; else // 0 is the left side, 1 is the right side dPosX += fDigiBdfPar->GetSignalSpeed() * ( xDigiB->GetTime() - xDigiA->GetTime() )/2.0; // For Z always just take the one of the channel itself( in fact its gap one) dPosZ = fChannelInfo->GetZ(); // Now check if a hit/cluster is already started if( 0 < iNbChanInHit) { if( iLastChan == iCh - 1 ) { fhDigTimeDifClust->Fill( dTime - dLastTime ); fhDigSpacDifClust->Fill( dPosX - dLastPosX ); fhDigDistClust->Fill( dPosX - dLastPosX, dTime - dLastTime ); } // if( iLastChan == iCh - 1 ) // a cluster is already started => check distance in space/time // For simplicity, just check along strip direction for now if( TMath::Abs( dTime - dLastTime) < dMaxTimeDist && iLastChan == iCh - 1 && TMath::Abs( dPosX - dLastPosX) < dMaxSpaceDist ) { // Add to cluster/hit dWeightedTime += dTime*dTotS; dWeightedPosX += dPosX*dTotS; dWeightedPosY += dPosY*dTotS; dWeightedPosZ += dPosZ*dTotS; dWeightsSum += dTotS; iNbChanInHit += 1; // if one of the side digi comes from a CbmTofPoint not already found // in this cluster, save its pointer Bool_t bFoundA = kFALSE; Bool_t bFoundB = kFALSE; if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() ) bFoundA = kTRUE; if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) bFoundB = kTRUE; } // for( Int else for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiA->GetLinks() ) bFoundA = kTRUE; if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiB->GetLinks() ) bFoundB = kTRUE; } // for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) // CbmTofPoint pointer for 1st digi not found => save it if( kFALSE == bFoundA ) vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // CbmTofPoint pointer for 2nd digi not found => save it if( kFALSE == bFoundB ) vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if current Digis compatible with last fired chan else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fChannelInfo->GetSizey()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); // Start a new hit dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if current Digis compatible with last fired chan } // if( 0 < iNbChanInHit) else { // first fired strip in this RPC dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if( 0 < iNbChanInHit) iLastChan = iCh; dLastPosX = dPosX; dLastPosY = dPosY; dLastTime = dTime; } // if( 2 == fStorDigi[iSmType][iSm][iRpc][iCh].size() ) // fStorDigi[iSmType][iSm][iRpc][iCh].clear(); fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { // Vertical strips => Y comes from bottom top time difference for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { if( 0 < fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); // if( 2 == fStorDigi[iSmType][iSm][iRpc][iCh].size() ) if( 2 == fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { if( (fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; continue; } // 2 Digis = both sides present CbmTofDetectorInfo xDetInfo(kTOF, iSmType, iSm, iRpc, 0, iCh + 1); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); // CbmTofDigi * xDigiA = fStorDigi[iSmType][iSm][iRpc][iCh][0]; // CbmTofDigi * xDigiB = fStorDigi[iSmType][iSm][iRpc][iCh][1]; CbmTofDigi * xDigiA = fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][0]; CbmTofDigi * xDigiB = fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh][1]; // The "Strip" time is the mean time between each end dTime = xDigiA->GetTime() + xDigiB->GetTime(); dTime /= 2.0; // Weight is the total charge => sum of both ends ToT dTotS = xDigiA->GetTot() + xDigiB->GetTot(); // X position is just the center of the channel dPosX = fChannelInfo->GetX(); // Y position form strip ends time difference dPosY = fChannelInfo->GetY(); if( 1 == xDigiA->GetSide() ) // 0 is the top side, 1 is the bottom side dPosY += fDigiBdfPar->GetSignalSpeed() * ( xDigiA->GetTime() - xDigiB->GetTime() )/2.0; else // 0 is the bottom side, 1 is the top side dPosY += fDigiBdfPar->GetSignalSpeed() * ( xDigiB->GetTime() - xDigiA->GetTime() )/2.0; // For Z always just take the one of the channel itself( in fact its gap one) dPosZ = fChannelInfo->GetZ(); // Now check if a hit/cluster is already started if( 0 < iNbChanInHit) { if( iLastChan == iCh - 1 ) { fhDigTimeDifClust->Fill( dTime - dLastTime ); fhDigSpacDifClust->Fill( dPosY - dLastPosY ); fhDigDistClust->Fill( dPosY - dLastPosY, dTime - dLastTime ); } // if( iLastChan == iCh - 1 ) // a cluster is already started => check distance in space/time // For simplicity, just check along strip direction for now if( TMath::Abs( dTime - dLastTime) < dMaxTimeDist && iLastChan == iCh - 1 && TMath::Abs( dPosY - dLastPosY) < dMaxSpaceDist ) { // Add to cluster/hit dWeightedTime += dTime*dTotS; dWeightedPosX += dPosX*dTotS; dWeightedPosY += dPosY*dTotS; dWeightedPosZ += dPosZ*dTotS; dWeightsSum += dTotS; iNbChanInHit += 1; // if one of the side digi comes from a CbmTofPoint not already found // in this cluster, save its pointer Bool_t bFoundA = kFALSE; Bool_t bFoundB = kFALSE; if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() ) bFoundA = kTRUE; if( ((CbmTofPoint*)(vPtsRef[iPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) bFoundB = kTRUE; } // for( Int else for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) { if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiA->GetLinks() ) bFoundA = kTRUE; if( vPtsRef[iPtRef] == (CbmTofPoint*)xDigiB->GetLinks() ) bFoundB = kTRUE; } // for( Int_t iPtRef = 0; iPtRef < vPtsRef.size(); iPtRef++) // CbmTofPoint pointer for 1st digi not found => save it if( kFALSE == bFoundA ) vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // CbmTofPoint pointer for 2nd digi not found => save it if( kFALSE == bFoundB ) vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if current Digis compatible with last fired chan else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fChannelInfo->GetSizey()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); // Start a new hit dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if current Digis compatible with last fired chan } // if( 0 < iNbChanInHit) else { // first fired strip in this RPC dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { if( ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() != ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // if other side come from a different Track, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) else if( xDigiA->GetLinks() != xDigiB->GetLinks() ) // if other side come from a different TofPoint, // also save the pointer on CbmTofPoint vPtsRef.push_back( (CbmTofPoint*)(xDigiB->GetLinks()) ); } // else of if( 0 < iNbChanInHit) iLastChan = iCh; dLastPosX = dPosX; dLastPosY = dPosY; dLastTime = dTime; } // if( 2 <= fStorDigi[iSmType][iSm][iRpc][iCh].size() ) // fStorDigi[iSmType][iSm][iRpc][iCh].clear(); fStorDigi[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 == iChType) else { LOG(ERROR)<<"CbmTofSimpClusterizer::BuildClusters => Cluster building " <<"from digis to hits not implemented for pads, Sm type " <GetChanOrient( iSmType, iRpc ) ) { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fChannelInfo->GetSizey()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); // Simple errors, not properly done at all for now // Right way of doing it should take into account the weight distribution // and real system time resolution TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSignalSpeed(), // Use the electronics resolution fDigiBdfPar->GetNbGaps( iSmType, iRpc)* fDigiBdfPar->GetGapSize( iSmType, iRpc)/10.0 / // Change gap size in cm TMath::Sqrt(12.0) ); // Use full RPC thickness as "Channel" Z size Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = 0; Int_t iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); // Int_t iRefId = 0; // Index of the correspondng TofPoint new((*fTofHitsColl)[fiNbHits]) CbmTofHit( iDetId, hitPos, hitPosErr, iRefId, dWeightedTime, vPtsRef.size() ); // flag = number of TofPoints generating the cluster // Using the SetLinks/GetLinks of the TofHit class seems to conflict // with something in littrack QA // ((CbmTofHit*)fTofHitsColl->At(fiNbHits))->SetLinks(vPtsRef[0]); fiNbHits++; // For Histogramming fviClusterSize[iSmType][iRpc].push_back(iNbChanInHit); fviTrkMul[iSmType][iRpc].push_back( vPtsRef.size() ); fvdX[iSmType][iRpc].push_back(dWeightedPosX); fvdY[iSmType][iRpc].push_back(dWeightedPosY); fvdDifX[iSmType][iRpc].push_back( vPtsRef[0]->GetX() - dWeightedPosX); fvdDifY[iSmType][iRpc].push_back( vPtsRef[0]->GetY() - dWeightedPosY); fvdDifCh[iSmType][iRpc].push_back( fGeoHandler->GetCell( vPtsRef[0]->GetDetectorID() ) -1 -iLastChan ); vPtsRef.clear(); } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 < iNbChanInHit) } // for each sm/rpc pair } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) //cout<<"test"<