/** @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 "CbmTofDetectorId_v14a.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" #include "CbmMatch.h" // FAIR classes and includes #include "FairRootManager.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairLogger.h" #include "FairMCEventHeader.h" // from CbmStsDigitize, for GetEventInfo #include "FairEventHeader.h" // from CbmStsDigitize, for GetEventInfo #include "FairRunSim.h" // from CbmStsDigitize, for GetEventInfo // 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 "TProfile.h" #include "TDirectory.h" #include "TROOT.h" #include "TGeoManager.h" // C++ Classes and includes // const Int_t DetMask = 4194303; (VF) not used const Int_t nbClWalkBinX=20; //const Int_t nbClWalkBinY=41; (VF) not used // choose odd number to have central bin symmetric around 0 //const Double_t WalkNHmin=100; (VF) not used // minimal number of hits in bin for walk correction Double_t TOTMax=5.E4; Double_t TOTMin=2.E4; const Double_t TTotMean=2.E4; const Int_t nbClDelTofBinX=50; //const Int_t nbClDelTofBinY=49; (VF) not used //const Double_t DelTofMax=5000.; (VF) not used //const Int_t nbCldXdYBinX=49; (VF) not used //const Int_t nbCldXdYBinY=49; (VF) not used //const Double_t dXdYMax=10.; (VF) not used const Int_t iNTrg=1; //const Double_t Zref = 200.; (VF) not used // distance of projection plane to target // C++ Classes and includes /************************************************************************************/ CbmTofSimpClusterizer::CbmTofSimpClusterizer(): FairTask("CbmTofSimpClusterizer"), fGeoHandler(new CbmTofGeoHandler()), fTofId(NULL), fDigiPar(NULL), fChannelInfo(NULL), fDigiBdfPar(NULL), fdParFeeTimeRes(0.0), fdParSystTimeRes(0.0), fTofPointsColl(NULL), fMcTracksColl(NULL), fTofDigisColl(NULL), fTofHitsColl(NULL), fTofDigiMatchColl(NULL), fiNbHits(0), fVerbose(1), fStorDigi(), fStorDigiExp(), fStorDigiInd(), fviClusterMul(), fviClusterSize(), fviTrkMul(), fvdX(), fvdY(), fvdDifX(), fvdDifY(), fvdDifCh(), fsHistoOutFilename(""), fhClustBuildTime(NULL), fhHitsPerTracks(NULL), fhPtsPerHit(NULL), fhTimeResSingHits(NULL), fhTimeResSingHitsB(NULL), fhTimePtVsHits(NULL), fhClusterSize(NULL), fhClusterSizeType(NULL), fhTrackMul(NULL), fhClusterSizeMulti(NULL), fhTrk1MulPos(NULL), fhHiTrkMulPos(NULL), fhAllTrkMulPos(NULL), fhMultiTrkProbPos(NULL), fhDigSpacDifClust(NULL), fhDigTimeDifClust(NULL), fhDigDistClust(NULL), fhClustSizeDifX(NULL), fhClustSizeDifY(NULL), fhChDifDifX(NULL), fhChDifDifY(NULL), fhRpcDigiCor(), fhRpcCluMul(), fhRpcSigPropSpeed(), fhRpcCluPosition(), fhRpcCluTOff(), fhRpcCluTrms(), fhRpcCluTot(), fhRpcCluSize(), fhRpcCluAvWalk(), fhRpcCluWalk(), fhTRpcCluMul(), fhTRpcCluPosition(), fhTRpcCluTOff(), fhTRpcCluTot(), fhTRpcCluSize(), fhTRpcCluAvWalk(), fhTRpcCluDelTof(), fhTRpcCludXdY(), fhTRpcCluWalk(), fhTrgdT(), fvCPSigPropSpeed(), fvCPDelTof(), fvCPTOff(), fvCPTotGain(), fvCPWalk(), fiNbSameSide(0), fhNbSameSide(NULL), fhNbDigiPerChan(NULL), fStart(), fStop(), fTimer(), fiNofEvents(0.), fdNofDigisTot(0.), fdNofHitsTot(0.), fdTimeTot(0.), dTRef(0.), fdTRefMax(0.), fCalMode(0), fCalTrg(0), fCalSmType(0), fdCaldXdYMax(0.), fTRefMode(0), fTRefHits(0), fPosYMaxScal(0.), fTRefDifMax(0.), fTotMax(0.), fTotMin(0.), fOutTimeFactor(1.), fCalParFileName(""), fCalParFile(NULL), fbMcTrkMonitor(kFALSE) { } CbmTofSimpClusterizer::CbmTofSimpClusterizer(const char *name, Int_t verbose): FairTask(TString(name),verbose), fGeoHandler(new CbmTofGeoHandler()), fTofId(NULL), fDigiPar(NULL), fChannelInfo(NULL), fDigiBdfPar(NULL), fdParFeeTimeRes(0.0), fdParSystTimeRes(0.0), fTofPointsColl(NULL), fMcTracksColl(NULL), fTofDigisColl(NULL), fTofHitsColl(NULL), fTofDigiMatchColl(NULL), fiNbHits(0), fVerbose(verbose), fStorDigi(), fStorDigiExp(), fStorDigiInd(), fviClusterMul(), fviClusterSize(), fviTrkMul(), fvdX(), fvdY(), fvdDifX(), fvdDifY(), fvdDifCh(), fsHistoOutFilename(""), fhClustBuildTime(NULL), fhHitsPerTracks(NULL), fhPtsPerHit(NULL), fhTimeResSingHits(NULL), fhTimeResSingHitsB(NULL), fhTimePtVsHits(NULL), fhClusterSize(NULL), fhClusterSizeType(NULL), fhTrackMul(NULL), fhClusterSizeMulti(NULL), fhTrk1MulPos(NULL), fhHiTrkMulPos(NULL), fhAllTrkMulPos(NULL), fhMultiTrkProbPos(NULL), fhDigSpacDifClust(NULL), fhDigTimeDifClust(NULL), fhDigDistClust(NULL), fhClustSizeDifX(NULL), fhClustSizeDifY(NULL), fhChDifDifX(NULL), fhChDifDifY(NULL), fhRpcDigiCor(), fhRpcCluMul(), fhRpcSigPropSpeed(), fhRpcCluPosition(), fhRpcCluTOff(), fhRpcCluTrms(), fhRpcCluTot(), fhRpcCluSize(), fhRpcCluAvWalk(), fhRpcCluWalk(), fhTRpcCluMul(), fhTRpcCluPosition(), fhTRpcCluTOff(), fhTRpcCluTot(), fhTRpcCluSize(), fhTRpcCluAvWalk(), fhTRpcCluDelTof(), fhTRpcCludXdY(), fhTRpcCluWalk(), fhTrgdT(), fvCPSigPropSpeed(), fvCPDelTof(), fvCPTOff(), fvCPTotGain(), fvCPWalk(), fiNbSameSide(0), fhNbSameSide(NULL), fhNbDigiPerChan(NULL), fStart(), fStop(), fTimer(), fiNofEvents(0.), fdNofDigisTot(0.), fdNofHitsTot(0.), fdTimeTot(0.), dTRef(0.), fdTRefMax(0.), fCalMode(0), fCalTrg(0), fCalSmType(0), fdCaldXdYMax(0.), fTRefMode(0), fTRefHits(0), fPosYMaxScal(0.), fTRefDifMax(0.), fTotMax(0.), fTotMin(0.), fOutTimeFactor(1.), fCalParFileName(""), fCalParFile(NULL), fbMcTrkMonitor(kFALSE) { } 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 == InitCalibParameter() ) return kFATAL; if( kFALSE == CreateHistos() ) return kFATAL; return kSUCCESS; } void CbmTofSimpClusterizer::SetParContainers() { LOG(info)<<" CbmTofSimpClusterizer => Get the digi parameters for tof"; // Get Base Container FairRunAna* ana = FairRunAna::Instance(); FairRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); LOG(info)<<" CbmTofSimpClusterizer::SetParContainers found " << fDigiPar->GetNrOfModules() << " cells " ; fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); } void CbmTofSimpClusterizer::Exec(Option_t* /*option*/) { // Start timer counter fTimer.Start(); fTofHitsColl->Clear("C"); // fTofDigiMatchColl->Clear("C"); // Not enough => CbmMatch has no Clear functions!! fTofDigiMatchColl->Delete(); fiNbHits = 0; Int_t iNbTofDigi = fTofDigisColl->GetEntries(); LOG(debug)<<" CbmTofSimpClusterizer => New event with " << iNbTofDigi <<" digis "; fStart.Set(); BuildClusters(); fStop.Set(); FillHistos(); // --- Update Counters fTimer.Stop(); fiNofEvents++; fdNofDigisTot += iNbTofDigi; fdNofHitsTot += fiNbHits; fdTimeTot += fTimer.RealTime(); } void CbmTofSimpClusterizer::Finish() { WriteHistos(); // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method DeleteHistos(); // Final printout for reference LOG(info) << "=====================================" ; LOG(info) << GetName() << ": Run summary (Time includes Hist filling)" ; LOG(info) << "Events processed : " << fiNofEvents ; LOG(info) << "Digis / event : " << fdNofDigisTot / static_cast(fiNofEvents) ; LOG(info) << "Hits / event : " << fdNofHitsTot / static_cast(fiNofEvents) ; LOG(info) << "Digis per Hits : " << fdNofDigisTot / fdNofHitsTot ; LOG(info) << "Time per event : " << fdTimeTot / static_cast(fiNofEvents) << " " ; LOG(info) << "=====================================" ; } /************************************************************************************/ // Functions common for all clusters approximations Bool_t CbmTofSimpClusterizer::RegisterInputs() { FairRootManager *fManager = FairRootManager::Instance(); /** VF: The task should run without MC input fTofPointsColl = (TClonesArray *) fManager->GetObject("TofPoint"); if( NULL == fTofPointsColl) { LOG(error)<<"CbmTofSimpClusterizer::RegisterInputs => Could not get the TofPoint TClonesArray!!!"; return kFALSE; } // if( NULL == fTofPointsColl) **/ fTofPointsColl = nullptr; fMcTracksColl = (TClonesArray *) fManager->GetObject("MCTrack"); if( NULL == fMcTracksColl) { LOG(info)<<"CbmTofSimpClusterizer: No MCTrack array."; } // if( NULL == fMcTracksColl) fTofDigisColl = (TClonesArray *) fManager->GetObject("TofDigi"); if( NULL == fTofDigisColl) { LOG(error)<<"CbmTofSimpClusterizer::RegisterInputs => Could not get the TofDigi TClonesArray!!!"; return kFALSE; } // if( NULL == fTofDigisColl) return kTRUE; } Bool_t CbmTofSimpClusterizer::RegisterOutputs() { FairRootManager* rootMgr = FairRootManager::Instance(); fTofHitsColl = new TClonesArray("CbmTofHit"); // Flag check to control whether digis are written in output root file rootMgr->Register( "TofHit","Tof", fTofHitsColl, IsOutputBranchPersistent("TofHit")); fTofDigiMatchColl = new TClonesArray("CbmMatch",100); rootMgr->Register( "TofHitDigiMatch", "Tof", fTofDigiMatchColl, IsOutputBranchPersistent("TofHitDigiMatch")); return kTRUE; } Bool_t CbmTofSimpClusterizer::InitParameters() { // Initialize the TOF GeoHandler Bool_t isSimulation=kFALSE; Int_t iGeoVersion = fGeoHandler->Init(isSimulation); LOG(info)<<"CbmTofSimpClusterizer::InitParameters with GeoVersion "< iGeoVersion ) { LOG(error)<<"CbmTofSimpClusterizer::InitParameters => Only compatible with geometries after v12b !!!" ; return kFALSE; } if(NULL != fTofId) LOG(info)<<"CbmTofSimpClusterizer::InitParameters with GeoVersion "<GetGeoVersion(); else{ switch(iGeoVersion){ case k12b: fTofId = new CbmTofDetectorId_v12b(); break; case k14a: fTofId = new CbmTofDetectorId_v14a(); break; default: LOG(error)<<"CbmTofSimpClusterizer::InitParameters => Invalid geometry!!!"<GetFeeTimeRes(); fdParSystTimeRes = 0.080; LOG(info)<<" CbmTofSimpClusterizer::InitParameters found Tres FEE = " << fdParFeeTimeRes << " ns " ; LOG(info)<<" CbmTofSimpClusterizer::InitParameters found Tres Syst = " << fdParSystTimeRes << " ns " ; return kTRUE; } /************************************************************************************/ Bool_t CbmTofSimpClusterizer::InitCalibParameter() { // dimension and initialize calib parameter // Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes(); fvCPSigPropSpeed.resize( iNbSmTypes ); for (Int_t iT=0; iTGetNbRpc( iT ); fvCPSigPropSpeed[iT].resize(iNbRpc); for (Int_t iRpc=0; iRpcGetSigVel( iT, 0, iRpc ) ) fvCPSigPropSpeed[iT][iRpc] = fDigiBdfPar->GetSigVel( iT, 0, iRpc ); else fvCPSigPropSpeed[iT][iRpc] = fDigiBdfPar->GetSignalSpeed(); } // for (Int_t iT=0; iTGetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); fvCPTOff[iSmType].resize( iNbSm*iNbRpc ); fvCPTotGain[iSmType].resize( iNbSm*iNbRpc ); fvCPWalk[iSmType].resize( iNbSm*iNbRpc ); fvCPDelTof[iSmType].resize( iNbSm*iNbRpc ); for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) { for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { // LOG(info)<GetNbChan( iSmType, iRpc ); fvCPTOff[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); fvCPTotGain[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); fvCPWalk[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); Int_t nbSide =2 - fDigiBdfPar->GetChanType( iSmType, iRpc ); for (Int_t iCh=0; iChcd(); // <= To prevent histos from being sucked in by the param file of the TRootManager ! */ if(0Print(); fCalParFile->cd(); fCalParFile->ls(); */ 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++ ) { TH2F *htempPos_pfx =(TH2F*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx",iSmType,iSm,iRpc)); TH2F *htempTOff_pfx=(TH2F*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx",iSmType,iSm,iRpc)); TH2F *htempTot_pfx =(TH2F*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_pfx",iSmType,iSm,iRpc)); if(NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_pfx) { Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); Int_t iNbinTot = htempTot_pfx->GetNbinsX(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { Double_t YMean=((TProfile *)htempPos_pfx)->GetBinContent(iCh+1); //nh +1 empirical(?) Double_t TMean=((TProfile *)htempTOff_pfx)->GetBinContent(iCh+1); Double_t dTYOff=YMean/fvCPSigPropSpeed[iSmType][iRpc] ; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] += -dTYOff + TMean ; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] += +dTYOff + TMean ; Double_t TotMean=((TProfile *)htempTot_pfx)->GetBinContent(iCh+1); //nh +1 empirical(?) if(1 " << fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] <<", " << fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] <<", " << fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0] <<", NbinTot "<< iNbinTot ; TH1D *htempWalk0=(TH1D*)gDirectory->FindObjectAny( Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); TH1D *htempWalk1=(TH1D*)gDirectory->FindObjectAny( Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); if(NULL != htempWalk0 && NULL != htempWalk1 ) { // reinitialize Walk array LOG(info)<<"Initialize Walk correction for " <GetNbinsX() != nbClWalkBinX) LOG(error)<<"CbmTofSimpClusterizer::InitCalibParameter: Inconsistent Walk histograms" ; for( Int_t iBin = 0; iBin < nbClWalkBinX; iBin++ ) { fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iBin]=htempWalk0->GetBinContent(iBin+1); fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iBin]=htempWalk1->GetBinContent(iBin+1); LOG(debug1)<FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof",iSmType,iSm,iRpc,iTrg)); if (NULL==htmpDelTof) { LOG(info)<<" Histos " << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof", iSmType, iSm, iRpc, iTrg) << " not found. " ; continue; } LOG(info)<<" Load DelTof from histos "<< Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof",iSmType,iSm,iRpc,iTrg)<<"." ; for(Int_t iBx=0; iBxGetBinContent(iBx+1); } // copy Histo to memory TDirectory * curdir = gDirectory; gDirectory->cd( oldir->GetPath() ); TH1D *h1DelTof=(TH1D *)htmpDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Trg%02d_DelTof",iSmType,iSm,iRpc,iTrg)); LOG(info)<<" copy histo " <GetName() <<" to directory " <GetName() ; gDirectory->cd( curdir->GetPath() ); } } } } // fCalParFile->Delete(); gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager! LOG(info)<<"CbmTofSimpClusterizer::InitCalibParameter: initialization done" ; return kTRUE; } /************************************************************************************/ Bool_t CbmTofSimpClusterizer::LoadGeometry() { LOG(debug)<<"CbmTofSimpClusterizer::LoadGeometry starting for " <GetNrOfModules() << " geometrically known modules " ; Int_t iNrOfCells = fDigiPar->GetNrOfModules(); LOG(debug)<<"Digi Parameter container contains "<GetGeoVersion(); fGeoHandler->CheckGeometryVersion(); for (Int_t icell = 0; icell < iNrOfCells; ++icell) { Int_t cellId = fDigiPar->GetCellId(icell); // cellId is assigned in CbmTofCreateDigiPar fChannelInfo = fDigiPar->GetCell(cellId); Int_t smtype = fGeoHandler->GetSMType(cellId); Int_t smodule = fGeoHandler->GetSModule(cellId); Int_t module = fGeoHandler->GetCounter(cellId); Int_t cell = fGeoHandler->GetCell(cellId); Double_t x = fChannelInfo->GetX(); Double_t y = fChannelInfo->GetY(); Double_t z = fChannelInfo->GetZ(); Double_t dx = fChannelInfo->GetSizex(); Double_t dy = fChannelInfo->GetSizey(); LOG(debug1) << "-I- InitPar "<UseExpandedDigi() ) else { fStorDigi.resize( iNbSmTypes ); fStorDigiInd.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*iNbRpc ); fStorDigiInd[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++ ) { for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { Int_t iNbChan = fDigiBdfPar->GetNbChan( iSmType, iRpc ); fStorDigi[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); fStorDigiInd[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++ ) } // else of if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) LOG(debug)<<"CbmTofSimpClusterizer::LoadGeometry: Done!" ; 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*iNbRpc+iRpc].clear(); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc].clear(); } } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) fStorDigiExp[iSmType].clear(); fStorDigiInd[iSmType].clear(); } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) fStorDigiExp.clear(); fStorDigiInd.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*iNbRpc+iRpc].clear(); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc].clear(); } } // for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) fStorDigi[iSmType].clear(); fStorDigiInd[iSmType].clear(); } // for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) fStorDigi.clear(); fStorDigiInd.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(); if( fbMcTrkMonitor && fMcTracksColl ) { 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 ); } // if( fbMcTrkMonitor ) 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 => WRONG and not consistent with hit creation // Need to loop on digi match object, check the digi to pnt match object for each, etc.... // ====> Just comment this code, if anybody needs it, have to implement proper solution 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( UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++ ) { fhClusterSize->Fill( fviClusterSize[iSmType][iRpc][uCluster]); fhClusterSizeType->Fill( fviClusterSize[iSmType][iRpc][uCluster], 40*iSmType + iRpc ); if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) { fhTrackMul->Fill( fviTrkMul[iSmType][iRpc][uCluster] ); fhClusterSizeMulti->Fill( fviClusterSize[iSmType][iRpc][uCluster], fviTrkMul[iSmType][iRpc][uCluster] ); if( 1 == fviTrkMul[iSmType][iRpc][uCluster] ) fhTrk1MulPos->Fill( fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster] ); if( 1 < fviTrkMul[iSmType][iRpc][uCluster] ) fhHiTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster] ); fhAllTrkMulPos->Fill(fvdX[iSmType][iRpc][uCluster], fvdY[iSmType][iRpc][uCluster] ); } // if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) if( 1 == fviTrkMul[iSmType][iRpc][uCluster] ) { fhClustSizeDifX->Fill( fviClusterSize[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]); fhClustSizeDifY->Fill( fviClusterSize[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]); if( 1 == fviClusterSize[iSmType][iRpc][uCluster] ) { fhChDifDifX->Fill( fvdDifCh[iSmType][iRpc][uCluster], fvdDifX[iSmType][iRpc][uCluster]); fhChDifDifY->Fill( fvdDifCh[iSmType][iRpc][uCluster], fvdDifY[iSmType][iRpc][uCluster]); } } } // for( UInt_t uCluster = 0; uCluster < fviClusterSize[iSmType][iRpc].size(); uCluster++ ) 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() { if( "" == fsHistoOutFilename ) { LOG(info) << "CbmTofSimpClusterizer::WriteHistos => Control histograms will not be written to disk!" ; LOG(info) << "CbmTofSimpClusterizer::WriteHistos => To change this behavior please provide a full path " << "with the SetHistoFileName method" ; return kTRUE; } // if( "" == fsHistoOutFilename ) TDirectory * oldir = gDirectory; TFile *fHist = new TFile(fsHistoOutFilename,"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::SetHistoFileName( TString sFilenameIn ) { fsHistoOutFilename = sFilenameIn; 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() { // --- MC Event info (input file, entry number, start time) Int_t iInputNr = 0; Int_t iEventNr = 0; Double_t dEventTime = 0.; GetEventInfo(iInputNr, iEventNr, dEventTime); /* * 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(); LOG(debug)<<" BuildCluster with MaxTimeDist "<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 ); LOG(debug1)<<"CbmTofSimpClusterizer::BuildClusters: " <GetType()<<" " <GetSm()<<" " <GetRpc()<<" " <GetChannel()<<" " <GetTime(),pDigi->GetTot()) ; if( fDigiBdfPar->GetNbSmTypes() > pDigi->GetType() // prevent crash due to misconfiguration && fDigiBdfPar->GetNbSm( pDigi->GetType()) > pDigi->GetSm() && fDigiBdfPar->GetNbRpc( pDigi->GetType()) > pDigi->GetRpc() && fDigiBdfPar->GetNbChan(pDigi->GetType(),0) >pDigi->GetChannel() ){ fStorDigiExp[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(pDigi); fStorDigiInd[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(iDigInd); // apply calibration vectors pDigi->SetTime(pDigi->GetTime()- // calibrate Digi Time fvCPTOff[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()]); pDigi->SetTot(pDigi->GetTot() * // calibrate Digi ToT fvCPTotGain[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()]); // walk correction Double_t dTotBinSize = (TOTMax-TOTMin)/ 2. / nbClWalkBinX; Int_t iWx = (Int_t)((pDigi->GetTot()-TOTMin/2.)/dTotBinSize); if (0>iWx) iWx=0; if (iWx > (nbClWalkBinX - 1) ) iWx=nbClWalkBinX-1; Double_t dDTot = (pDigi->GetTot()-TOTMin/2.)/dTotBinSize-(Double_t)iWx-0.5; Double_t dWT = fvCPWalk[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx]; if(dDTot > 0) { if(iWx < nbClWalkBinX -1) { // linear interpolation to next bin dWT += dDTot * (fvCPWalk[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx+1] -fvCPWalk[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx]); } // if(iWx < nbClWalkBinX -1) }else // dDTot < 0, { if(0 < iWx) { // linear interpolation to next bin dWT -= dDTot * (fvCPWalk[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx-1] -fvCPWalk[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx]); } // if(0 < iWx) } pDigi->SetTime(pDigi->GetTime() - dWT); // calibrate Digi Time if(0) {//pDigi->GetType()==2 && pDigi->GetSm()==0){ LOG(info)<<"CbmTofSimpClusterizer::BuildClusters: CalDigi " <GetType()<<", Sm " <GetSm()<<", R " <GetRpc()<<", Ch " <GetChannel()<<", S " <GetSide()<<", T " <GetTime()<<", Tot " <GetTot() <<", TotGain "<< fvCPTotGain[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] <<", Walk "<GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx] ; LOG(info)<<" dDTot "<GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx-1] <<", "<GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx] <<", "<GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()] [iWx+1] <<" -> dWT = "<< dWT ; } } else { LOG(debug)<<"CbmTofSimpBeamClusterizer::BuildClusters: Skip Digi " <<" Type "<GetType()<<" "<< fDigiBdfPar->GetNbSmTypes() <<" Sm " <GetSm()<<" " << fDigiBdfPar->GetNbSm(pDigi->GetType()) <<" Rpc " <GetRpc()<<" "<< fDigiBdfPar->GetNbRpc(pDigi->GetType()) <<" Ch " <GetChannel()<<" "<GetNbChan(pDigi->GetType(),0) ; } } // 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()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(pDigi); fStorDigiInd[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()].push_back(iDigInd); } // 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 dWeightedTimeErr = 0.0; Double_t dWeightsSum = 0.0; std::vector< CbmTofPoint* > vPtsRef; std::vector< Int_t > vDigiIndRef; CbmTofCell *fTrafoCell=NULL; Int_t iTrafoCell=-1; Int_t iNbChanInHit = 0; // Last Channel Temp variables Int_t iLastChan = -1; // Double_t dLastPosX = 0.0; // -> Comment to remove warning because set but never used 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; /// Go to Top volume of the geometry in the GeoManager to make sure /// our nodes are found gGeoManager->CdTop(); 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 iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); LOG(debug2)<<"CbmTofSimpClusterizer::BuildClusters: RPC - Loop " << Form(" %3d %3d %3d %3d ",iSmType,iSm,iRpc,iChType) ; fviClusterMul[iSmType][iSm][iRpc]=0; Int_t iChId = 0; if( 0 == iChType ) { // Don't spread clusters over RPCs!!! dWeightedTime = 0.0; dWeightedPosX = 0.0; dWeightedPosY = 0.0; dWeightedPosZ = 0.0; dWeightedTimeErr = 0.0; dWeightsSum = 0.0; iNbChanInHit = 0; vPtsRef.clear(); // For safety reinitialize everything iLastChan = -1; // dLastPosX = 0.0; // -> Comment to remove warning because set but never used dLastPosY = 0.0; dLastTime = 0.0; fTrafoCell=NULL; iTrafoCell=-1; LOG(debug2)<<"CbmTofSimpClusterizer::BuildClusters: ChanOrient " << Form(" %3d %3d %3d %3d %3d ",iSmType,iSm,iRpc,fDigiBdfPar->GetChanOrient( iSmType, iRpc ),iNbCh) ; if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) { // Horizontal strips => X comes from left right time difference } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { // Vertical strips => Y comes from bottom top time difference for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { LOG(debug3)<<"CbmTofSimpClusterizer::BuildClusters: VDigisize" << Form(" T %3d Sm %3d R %3d Ch %3d Size %3zu ",iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; if( 0 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); while( 1 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: digis processing for " << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3zu ",iSmType,iSm,iRpc,iCh, fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; /* if (3 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) // needs more work! LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: digis ignored for " << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3d ",iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; */ while( (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; LOG(debug) << "CbmTofSimpClusterizer::BuildClusters: SameSide Hits! Times: " << (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetTime() << ", "<<(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetTime() <<", DeltaT " <<(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetTime() - (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetTime() ; LOG(debug2)<<" SameSide Erase fStor entries(d) "< fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) break; continue; } LOG(debug2) << "CbmTofSimpClusterizer::BuildClusters: digis processing for " << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3zu ",iSmType,iSm,iRpc,iCh, fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; if(2 > fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) break; /* Int_t iLastChId = iChId; // Save Last hit channel */ // 2 Digis = both sides present Int_t iCh1=iCh; if(fGeoHandler->GetGeoVersion() < k14a) iCh1= iCh1+1; //FIXME CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, iCh1); iChId = fTofId->SetDetectorInfo( xDetInfo ); Int_t iUCellId=CbmTofAddress::GetUniqueAddress(iSm,iRpc,iCh,0,iSmType); fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(error)<<"CbmTofSimpClusterizer::BuildClusters: no geometry info! " << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ",iSmType, iSm, iRpc, iCh, iChId,iUCellId) ; break; } LOG(debug1)<<"CbmTofSimpClusterizer::BuildClusters:" << Form(" T %3d Sm %3d R %3d Ch %3d size %3zu ", iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) << Form(" ChId: 0x%08x 0x%08x %p",iChId,iUCellId,fChannelInfo) ; 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 =0.5 * ( xDigiA->GetTime() + xDigiB->GetTime() ) ; // 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 from strip ends time difference //dPosY = fChannelInfo->GetY(); // For Z always just take the one of the channel itself( in fact its gap one) //dPosZ = fChannelInfo->GetZ(); TGeoNode *fNode= // prepare local->global trafo gGeoManager->FindNode(fChannelInfo->GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ()); if(NULL == fTrafoCell) { fTrafoCell=fChannelInfo; iTrafoCell=iCh; } //fNode->Print(); LOG(debug1)<GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ(),fNode,fChannelInfo,fTrafoCell) ; // fNode->Print(); // switch to local coordinates, (0,0,0) is in the center of counter ? dPosX=((Double_t)(-iTrafoCell + iCh))*fChannelInfo->GetSizex(); dPosY=0.; dPosZ=0.; if( 1 == xDigiA->GetSide() ) // 0 is the top side, 1 is the bottom side dPosY += fvCPSigPropSpeed[iSmType][iRpc] * ( xDigiA->GetTime() - xDigiB->GetTime() )/2.0; else // 0 is the bottom side, 1 is the top side dPosY += fvCPSigPropSpeed[iSmType][iRpc] * ( xDigiB->GetTime() - xDigiA->GetTime() )/2.0; if( fChannelInfo->GetSizey()/2.0 < dPosY || -1*fChannelInfo->GetSizey()/2.0 > dPosY ) { // if outside of strip limits, the pair is bad => try to remove one end and check the next pair // (if another possibility exist) fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].begin()); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin()); continue; } // Pair leads to hit oustide of strip limits LOG(debug1) <<"CbmTofSimpClusterizer::BuildClusters: NbChanInHit" << Form(" %3d %3d %3d %p %2.0f Time %6.1f PosY %5.1f Svel %5.1f", iNbChanInHit,iCh,iLastChan,xDigiA,xDigiA->GetSide(), dTime,dPosY,fvCPSigPropSpeed[iSmType][iRpc]) << Form( ", Offs %5.1f, %5.1f",fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]) ; // 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; dWeightedTimeErr += dTotS * 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; // -> Comment to remove warning because set but never used // Bool_t bFoundB = kFALSE; // -> Comment to remove warning because set but never used if(NULL != fTofPointsColl){ // MC /** **Comment the full Block as not used anymore** **/ /* if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) for( UInt_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++) { //if( ((CbmTofPoint*)(vPtsRef[uPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiA->GetLinks()))->GetTrackID() ) // bFoundA = kTRUE; // -> Comment to remove warning because set but never used //if( ((CbmTofPoint*)(vPtsRef[uPtRef]))->GetTrackID() == ((CbmTofPoint*)(xDigiB->GetLinks()))->GetTrackID() ) // bFoundB = kTRUE; // -> Comment to remove warning because set but never used } // for( Int else for( UInt_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++) { // if( vPtsRef[uPtRef] == (CbmTofPoint*)xDigiA->GetLinks() ) // bFoundA = kTRUE; // -> Comment to remove warning because set but never used // if( vPtsRef[uPtRef] == (CbmTofPoint*)xDigiB->GetLinks() ) // bFoundB = kTRUE; // -> Comment to remove warning because set but never used } // for( Int_t uPtRef = 0; uPtRef < vPtsRef.size(); uPtRef++) */ // 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()) ); } // MC end vDigiIndRef.push_back( (Int_t )(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][0])); vDigiIndRef.push_back( (Int_t )(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][1])); LOG(debug1)<<" Add Digi and erase fStor entries(a): NbChanInHit "<< iNbChanInHit<<", " <global trafo gGeoManager->FindNode(fTrafoCell->GetX(),fTrafoCell->GetY(),fTrafoCell->GetZ()); TGeoNode *cNode = gGeoManager->GetCurrentNode(); TGeoRotation rotMatrix ( *gGeoManager->GetCurrentMatrix() ); hitpos[0]=fTrafoCell->GetX(); hitpos[1]=fTrafoCell->GetY(); hitpos[2]=fTrafoCell->GetZ(); gGeoManager->MasterToLocal(hitpos, hitpos_local); LOG(debug1)<GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ(),hitpos_local[0], hitpos_local[1], hitpos_local[2], tNode) ; hitpos_local[0] += dWeightedPosX; hitpos_local[1] += dWeightedPosY; hitpos_local[2] += dWeightedPosZ; gGeoManager->LocalToMaster(hitpos_local, hitpos); LOG(debug1)<< Form(" LTM for node %p, info %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)", cNode, fChannelInfo, hitpos_local[0], hitpos_local[1], hitpos_local[2], hitpos[0], hitpos[1], hitpos[2]) ; TVector3 hitPos(hitpos[0],hitpos[1],hitpos[2]); // 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 Double_t hiterr_local[3]; Double_t hiterr[3]; hiterr_local[0] = fChannelInfo->GetSizex()/TMath::Sqrt(12.0); // Single strips approximation hiterr_local[1] = fdParFeeTimeRes * fvCPSigPropSpeed[iSmType][iRpc]; // Use the electronics resolution hiterr_local[2] = 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 rotMatrix.LocalToMaster(hiterr_local, hiterr); // TVector3 hitPosErr( hiterr_local[0], hiterr_local[1], hiterr_local[2] ); TVector3 hitPosErr( hiterr[0], hiterr[1], hiterr[2] ); // Int_t iDetId = vPtsRef[0]->GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // calc mean ch from dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex(); Int_t iChm=iTrafoCell+floor(dWeightedPosX/fChannelInfo->GetSizex()); if(iChm<0 || iChm >iNbCh){ LOG(debug)<<"CbmTofSimpClusterizer::BuildClusters: Invalid mean channel"; } Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm,iRpc,iChm,0,iSmType); Int_t iRefId = 0; // Index of the correspondng TofPoint // if(NULL != fTofPointsColl) { iRefId = fTofPointsColl->IndexOf( vPtsRef[0] ); } LOG(debug)<<"CbmTofSimpClusterizer::BuildClusters: Save Hit " << Form(" %3d %3d 0x%08x %3d %3d %3d %f %f", fiNbHits,iNbChanInHit,iDetId,iCh,iLastChan,iRefId, dWeightedTime,dWeightedPosY) <<", DigiSize: "<AddLink(CbmLink(0.,vDigiIndRef.at(i), iEventNr, iInputNr)); } new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(*digiMatch); delete digiMatch; // 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); /* no TofPoint available for data! 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(); vDigiIndRef.clear(); // Start a new hit dWeightedTime = dTime*dTotS; dWeightedPosX = dPosX*dTotS; dWeightedPosY = dPosY*dTotS; dWeightedPosZ = dPosZ*dTotS; dWeightedTimeErr = dTotS * dTotS; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // Save next digi address LOG(DEBUG4)<<" Next fStor Digi "<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 { LOG(debug1)<<"CbmTofSimpClusterizer::BuildClusters: 1.Hit on channel " <GetLinks()) ); vDigiIndRef.push_back( (Int_t )(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][0])); vDigiIndRef.push_back( (Int_t )(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][1])); LOG(debug2)<<" Erase fStor entries(c) "<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; // -> Comment to remove warning because set but never used dLastPosY = dPosY; dLastTime = dTime; } // if( 2 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].clear(); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) LOG(debug2)<<"CbmTofSimpClusterizer::BuildClusters: finished V-RPC" << Form(" %3d %3d %3d %d",iSmType,iSm,iRpc,fTofHitsColl->GetEntries()) ; } // 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 ) ) { LOG(debug1)<<"CbmTofSimpClusterizer::BuildClusters: H-Hit " ; } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { LOG(debug2)<<"CbmTofSimpClusterizer::BuildClusters: V-Hit " ; // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; // simple error scaling with TOT // dWeightedTimeErr = TMath::Sqrt( dWeightedTimeErr ) * fdParSystTimeRes / dWeightsSum; // OR harcoded value dWeightedTimeErr = fdParSystTimeRes; Double_t hitpos_local[3]; Double_t hitpos[3]; TGeoNode *tNode= // prepare local->global trafo gGeoManager->FindNode(fTrafoCell->GetX(),fTrafoCell->GetY(),fTrafoCell->GetZ()); //TGeoNode *cNode = gGeoManager->GetCurrentNode(); TGeoRotation rotMatrix ( *gGeoManager->GetCurrentMatrix() ); hitpos[0]=fTrafoCell->GetX(); hitpos[1]=fTrafoCell->GetY(); hitpos[2]=fTrafoCell->GetZ(); gGeoManager->MasterToLocal(hitpos, hitpos_local); LOG(debug1)<GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ(),hitpos_local[0], hitpos_local[1], hitpos_local[2]) ; hitpos_local[0] += dWeightedPosX; hitpos_local[1] += dWeightedPosY; hitpos_local[2] += dWeightedPosZ; gGeoManager->LocalToMaster(hitpos_local, hitpos); LOG(debug1)<< Form(" LTM for V-node %p, info %p, tra %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f) [(%6.1f,%6.1f,%6.1f)]", tNode, fChannelInfo, fTrafoCell, hitpos_local[0], hitpos_local[1], hitpos_local[2], hitpos[0], hitpos[1], hitpos[2], fTrafoCell->GetX(),fTrafoCell->GetY(),fTrafoCell->GetZ()) ; TVector3 hitPos(hitpos[0],hitpos[1],hitpos[2]); // TestBeam 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 Double_t hiterr_local[3]; Double_t hiterr[3]; hiterr_local[0] = fTrafoCell->GetSizex()/TMath::Sqrt(12.0); // Single strips approximation hiterr_local[1] = fdParFeeTimeRes * fvCPSigPropSpeed[iSmType][iRpc]; // Use the electronics resolution hiterr_local[2] = 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 rotMatrix.LocalToMaster(hiterr_local, hiterr); // TVector3 hitPosErr( hiterr_local[0], hiterr_local[1], hiterr_local[2] ); TVector3 hitPosErr( hiterr[0], hiterr[1], hiterr[2] ); /* LOG(info)<< " Size X " << fTrafoCell->GetSizex() ; LOG(info)<< " Nb gaps " << fDigiBdfPar->GetNbGaps( iSmType, iRpc) << " gap size " << fDigiBdfPar->GetGapSize( iSmType, iRpc) ; LOG(info)<<"CbmTofSimpClusterizer::BuildClusters: Errors in local " << hiterr_local[0] << " "<< hiterr_local[1] << " " << hiterr_local[2] ; LOG(info)<<"CbmTofSimpClusterizer::BuildClusters: Errors in master " << hiterr[0] << " "<< hiterr[1] << " " << hiterr[2] ; */ // cout<<"a "<GetDetectorID()<GetDetectorID();// detID = pt->GetDetectorID() <= from TofPoint // Int_t iDetId = iChId; Int_t iChm=iTrafoCell+floor(dWeightedPosX/fChannelInfo->GetSizex()); if(iChm<0 || iChm >iNbCh){ LOG(debug)<<"CbmTofSimpClusterizer::BuildClusters: Invalid mean channel " <IndexOf( vPtsRef[0] ); LOG(debug)<<"CbmTofSimpClusterizer::BuildClusters: Save V-Hit " // << Form(" %3d %3d 0x%08x %3d %3d %3d 0x%08x", << Form(" %3d %3d %3d %3d %3d ", fiNbHits,iNbChanInHit,iDetId,iLastChan,iRefId) //vPtsRef.size(),vPtsRef[0]) // dWeightedTime,dWeightedPosY) <<", DigiSize: "<AddLink(CbmLink(0.,vDigiIndRef.at(i), iEventNr, iInputNr)); } new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(*digiMatch); delete digiMatch; // 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(); vDigiIndRef.clear(); } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 < iNbChanInHit) LOG(DEBUG4)<<" Fini-A "<UseExpandedDigi() ) else { LOG(error)<<" Compressed Digis not implemented ... "; } return kTRUE; } void CbmTofSimpClusterizer::GetEventInfo(Int_t& inputNr, Int_t& eventNr, Double_t& eventTime) { // --- In a FairRunAna, take the information from FairEventHeader if ( FairRunAna::Instance() ) { FairEventHeader* event = FairRunAna::Instance()->GetEventHeader(); inputNr = event->GetInputFileId(); eventNr = event->GetMCEntryNumber(); eventTime = event->GetEventTime(); } // --- In a FairRunSim, the input number and event time are always zero; // --- only the event number is retrieved. else { if ( ! FairRunSim::Instance() ) LOG(fatal) << GetName() << ": neither SIM nor ANA run." ; FairMCEventHeader* event = FairRunSim::Instance()->GetMCEventHeader(); inputNr = 0; eventNr = event->GetEventID(); eventTime = 0.; } }