/** @file CbmTofEventClusterizer.cxx ** @author nh ** @date 13.07.2018 ** adopted from ** @file CbmTofTestBeamClusterizer.cxx ** @file CbmTofSimpClusterizer.cxx ** @author Pierre-Alain Loizeau ** @date 23.08.2013 **/ #include "CbmTofEventClusterizer.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 #include "CbmMatch.h" #include "CbmEvent.h" #include "CbmVertex.h" #include "TTrbHeader.h" // CBMroot classes and includes #include "CbmMCTrack.h" // FAIR classes and includes #include "FairRootManager.h" #include "FairRootFileSink.h" #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairLogger.h" #include "FairEventHeader.h" // ROOT Classes and includes #include "TClonesArray.h" #include "TMath.h" #include "TLine.h" #include "TRandom3.h" #include "TF1.h" #include "TF2.h" #include "TVector3.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" #include "TProfile.h" #include "TDirectory.h" #include "TROOT.h" #include "TGeoManager.h" #include "TMinuit.h" #include "TRandom.h" // Constants definitions #include "CbmTofClusterizersDef.h" // C++ Classes and includes // Globals #include static Int_t iIndexDut = 0; static Double_t StartAnalysisTime = 0.; // const Double_t cLight=29.9792; // in cm/ns (VF) not used static Int_t SelMask=DetMask; static Double_t fdStartAna10s=0.; static Double_t dTLEvt=0.; static Int_t iNSpill=0; const Double_t fdSpillDuration = 4.; // in seconds const Double_t fdSpillBreak = 0.3; // in seconds static Bool_t bIsMC=kFALSE; // std::vector< CbmTofPoint* > vPtsRef; CbmTofEventClusterizer *CbmTofEventClusterizer::fInstance = 0; /************************************************************************************/ CbmTofEventClusterizer::CbmTofEventClusterizer() : CbmTofEventClusterizer("TestbeamClusterizer",0,0) { if ( !fInstance ) fInstance = this; } CbmTofEventClusterizer::CbmTofEventClusterizer(const char *name, Int_t verbose, Bool_t writeDataInOut): FairTask(TString(name),verbose), fGeoHandler(new CbmTofGeoHandler()), fTofId(NULL), fDigiPar(NULL), fChannelInfo(NULL), fDigiBdfPar(NULL), fTrbHeader(NULL), fTofPointsColl(NULL), fMcTracksColl(NULL), fTofDigisColl(NULL), fTofRawDigisColl(nullptr), fEventsColl(nullptr), fbWriteHitsInOut(writeDataInOut), fbWriteDigisInOut(writeDataInOut), fTofCalDigisColl(NULL), fTofHitsColl(NULL), fTofDigiMatchColl(NULL), fTofCalDigisCollOut(NULL), fTofHitsCollOut(NULL), fTofDigiMatchCollOut(NULL), fiNbHits(0), fVerbose(verbose), fStorDigi(), fStorDigiExp(), fStorDigiInd(), vDigiIndRef(), fviClusterMul(), fviClusterSize(), fviTrkMul(), fvdX(), fvdY(), fvdDifX(), fvdDifY(), fvdDifCh(), 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), fhCluMulCorDutSel(NULL), fhEvCluMul(NULL), fhRpcDigiCor(), fhRpcDigiMul(), fhRpcDigiStatus(), fhRpcDigiDTLD(), fhRpcCluMul(), fhRpcCluRate(), fhRpcCluRate10s(), fhRpcCluPosition(), fhRpcCluPositionEvol(), fhRpcCluTimeEvol(), fhRpcCluDelPos(), fhRpcCluDelMatPos(), fhRpcCluTOff(), fhRpcCluDelTOff(), fhRpcCluDelMatTOff(), fhRpcCluTrms(), fhRpcCluTot(), fhRpcCluSize(), fhRpcCluAvWalk(), fhRpcCluAvLnWalk(), fhRpcCluWalk(), fhSmCluPosition(), fhSmCluTOff(), fhSmCluSvel(), fhSmCluFpar(), fhRpcDTLastHits(), fhRpcDTLastHits_Tot(), fhRpcDTLastHits_CluSize(), fhTRpcCluMul(), fhTRpcCluPosition(), fhTRpcCluTOff(), fhTRpcCluTofOff(), fhTRpcCluTot(), fhTRpcCluSize(), fhTRpcCluAvWalk(), fhTRpcCluDelTof(), fhTRpcCludXdY(), fhTRpcCluWalk(), fhTRpcCluWalk2(), fhTSmCluPosition(), fhTSmCluTOff(), fhTSmCluTRun(), fhTRpcCluTOffDTLastHits(), fhTRpcCluTotDTLastHits(), fhTRpcCluSizeDTLastHits(), fhTRpcCluMemMulDTLastHits(), fhSeldT(), fvCPDelTof(), fvCPTOff(), fvCPTotGain(), fvCPTotOff(), fvCPWalk(), fvLastHits(), fvDeadStrips(), fvTimeLastDigi(), fiNbSameSide(0), fhNbSameSide(NULL), fhNbDigiPerChan(NULL), fStart(), fStop(), dTRef(0.), fdTRefMax(0.), fCalMode(0), fCalSel(0), fCalSmAddr(0), fdCaldXdYMax(0.), fiCluMulMax(0), fTRefMode(0), fTRefHits(0), fIdMode(0), fDutId(0), fDutSm(0), fDutRpc(0), fDutAddr(0), fSelId(0), fSelSm(0), fSelRpc(0), fSelAddr(0), fiBeamRefType(0), fiBeamRefSm(0), fiBeamRefDet(0), fiBeamRefAddr(0), fiBeamRefMulMax(1), fiBeamAddRefMul(0), fSel2Id(0), fSel2Sm(0), fSel2Rpc(0), fSel2Addr(0), fSel2MulMax(1), fDetIdIndexMap(), fviDetId(), fPosYMaxScal(0.), fTRefDifMax(0.), fTotMax(0.), fTotMin(0.), fTotOff(0.), fTotMean(0.), fdDelTofMax(60.), fTotPreRange(0.), fMaxTimeDist(0.), fdChannelDeadtime(0.), fdMemoryTime(0.), fdYFitMin(1.E6), fdToDAv(0.033), // in ns/cm fEnableMatchPosScaling(kTRUE), fEnableAvWalk(kFALSE), fbPs2Ns(kFALSE), fCalParFileName(""), fOutHstFileName(""), fCalParFile(NULL), fiNevtBuild(0), fiMsgCnt(100), fdTOTMax(50.), fdTOTMin(0.), fdTTotMean(2.), fdMaxTimeDist(0.), fdMaxSpaceDist(0.), fdEvent(0), fbSwapChannelSides(kFALSE), fiOutputTreeEntry(0), fiFileIndex(0), fbAlternativeBranchNames(kFALSE) { if ( !fInstance ) fInstance = this; } CbmTofEventClusterizer::~CbmTofEventClusterizer() { if( fGeoHandler ) delete fGeoHandler; if( fInstance==this ) fInstance = 0; // DeleteHistos(); // <-- if needed ? } /************************************************************************************/ // FairTasks inherited functions InitStatus CbmTofEventClusterizer::Init() { LOG(info) << "CbmTofEventClusterizer initializing... expect Digis in ns units! "; 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; switch(fIdMode){ case 0: fDutAddr =CbmTofAddress::GetUniqueAddress(fDutSm,fDutRpc,0,0,fDutId); fSelAddr =CbmTofAddress::GetUniqueAddress(fSelSm,fSelRpc,0,0,fSelId); fSel2Addr=CbmTofAddress::GetUniqueAddress(fSel2Sm,fSel2Rpc,0,0,fSel2Id); fiBeamRefAddr=CbmTofAddress::GetUniqueAddress(fiBeamRefSm,fiBeamRefDet,0,0,fiBeamRefType); break; case 1: SelMask=ModMask; fDutAddr =CbmTofAddress::GetUniqueAddress(fDutSm,0,0,0,fDutId); fSelAddr =CbmTofAddress::GetUniqueAddress(fSelSm,0,0,0,fSelId); fSel2Addr=CbmTofAddress::GetUniqueAddress(fSel2Sm,0,0,0,fSel2Id); fiBeamRefAddr=CbmTofAddress::GetUniqueAddress(fiBeamRefSm,0,0,0,fiBeamRefType); break; } iIndexDut=fDigiBdfPar->GetDetInd(fDutAddr); return kSUCCESS; } void CbmTofEventClusterizer::SetParContainers() { LOG(info)<<"=> Get the digi parameters for tof"; //LOG(warning)<<"Return without action"; //return; // Get Base Container FairRunAna* ana = FairRunAna::Instance(); FairRuntimeDb* rtdb=ana->GetRuntimeDb(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); LOG(info)<<"found " << fDigiPar->GetNrOfModules() << " cells " ; fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); } void CbmTofEventClusterizer::Exec(Option_t* option) { FairEventHeader* eventHeader = NULL; if ( FairRunAna::Instance() ) eventHeader = FairRunAna::Instance()->GetEventHeader(); if( NULL != eventHeader) { LOG(debug) << "MC event entry " << eventHeader->GetMCEntryNumber() <<" at t = " << eventHeader->GetEventTime(); if(eventHeader->GetMCEntryNumber() > 0) { // first MC event will not be recognized bIsMC=kTRUE; } } if(fTofCalDigisCollOut) fTofCalDigisCollOut->Delete(); if(fEventsColl) { LOG(debug)<<"CbmTofEventClusterizer::Exec => New time slice with " << fEventsColl->GetEntriesFast() << " events, " << fTofRawDigisColl->GetEntriesFast() << " TOF digis " ; Int_t iNbHits=0; Int_t iNbCalDigis = 0; fTofDigiMatchCollOut->Delete(); // costly, FIXME for(Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++) { CbmEvent* tEvent = dynamic_cast(fEventsColl->At(iEvent)); if (fTofDigisColl) fTofDigisColl->Clear("C"); Int_t iNbDigis=0; for (Int_t iDigi = 0; iDigi < tEvent->GetNofData(kTofDigi); iDigi++) { Int_t iDigiIndex = static_cast(tEvent->GetIndex(kTofDigi, iDigi)); CbmTofDigiExp* tDigi = dynamic_cast(fTofRawDigisColl->At(iDigiIndex)); new((*fTofDigisColl)[iNbDigis++]) CbmTofDigiExp(*tDigi); } ExecEvent(option); // --- In event-by-event mode: copy caldigis, hits and matches to output array and register them to event // Int_t iDigi0=iNbCalDigis; //starting index of current event (VF) not used for (Int_t index = 0; index < fTofCalDigisColl->GetEntriesFast(); index++){ CbmTofDigiExp* tDigi = dynamic_cast(fTofCalDigisColl->At(index)); tEvent->AddData(kTofCalDigi, iNbCalDigis); new((*fTofCalDigisCollOut)[iNbCalDigis++]) CbmTofDigiExp(*tDigi); } for (Int_t index = 0; index < fTofHitsColl->GetEntriesFast(); index++){ CbmTofHit* pHit = (CbmTofHit *)fTofHitsColl->At(index); new((*fTofHitsCollOut)[iNbHits]) CbmTofHit(*pHit); tEvent->AddData(kTofHit, iNbHits); CbmMatch* pDigiMatch=(CbmMatch *)fTofDigiMatchColl->At( index ); // update content of match object, not necessary if event definition is kept ! /* for (Int_t iLink=0; iLinkGetNofLinks(); iLink++) { // loop over digis CbmLink Link = pDigiMatch->GetLink(iLink); Link.SetIndex(Link.GetIndex()+iDigi0); } */ new((*fTofDigiMatchCollOut)[iNbHits]) CbmMatch(*pDigiMatch); iNbHits++; } fTofDigisColl->Delete(); fTofCalDigisColl->Delete();//Clear("C"); //otherwise memoryleak! FIXME fTofHitsColl->Clear("C"); fTofDigiMatchColl->Delete();//Clear("C"); } } else { fTofDigisColl=fTofRawDigisColl; if(bIsMC) { // add fake diamond digis, this part should be moved to digitizer Int_t iNbDigis=fTofDigisColl->GetEntriesFast(); UInt_t uChanUId=0x00005006; Double_t dHitTime=gRandom->Gaus(0.,0.04); const Double_t dHitTot=2.; CbmTofDigiExp* tDigi = new CbmTofDigiExp(uChanUId, dHitTime, dHitTot); new((*fTofDigisColl)[iNbDigis++]) CbmTofDigiExp(*tDigi); uChanUId=0x00805006; CbmTofDigiExp* tDigi1 = new CbmTofDigiExp(uChanUId, dHitTime, dHitTot); new((*fTofDigisColl)[iNbDigis++]) CbmTofDigiExp(*tDigi1); LOG(debug) << Form("Add fake diamond digis 0x%08x in event mode with t = %7.3f",uChanUId,dHitTime); } ExecEvent(option); } } void CbmTofEventClusterizer::ExecEvent(Option_t* /*option*/) { // Clear output arrays fTofCalDigisColl->Delete(); //otherwise memoryleak if 'CbmDigi::fMatch' points to valid MC match objects (simulation)! FIXME fTofHitsColl->Clear("C"); //fTofHitsColl->Delete(); // Computationally costly!, but hopefully safe //for (Int_t i=0; iGetEntries(); i++) ((CbmMatch *)(fTofDigiMatchColl->At(i)))->ClearLinks(); // FIXME, try to tamper memory leak (did not help) //fTofDigiMatchColl->Clear("C+L"); // leads to memory leak fTofDigiMatchColl->Delete(); FairRootFileSink* bla =(FairRootFileSink *)FairRootManager::Instance()->GetSink() ; if (bla) fiOutputTreeEntry = ((FairRootFileSink *)FairRootManager::Instance()->GetSink())->GetOutTree()->GetEntries(); fiNbHits = 0; fStart.Set(); BuildClusters(); MergeClusters(); fStop.Set(); fdEvent++; FillHistos(); // fTofDigisColl->RemoveAll(); } /************************************************************************************/ void CbmTofEventClusterizer::Finish() { if (fdEvent<100 ) return; // don't save histos with insufficient statistics WriteHistos(); // Prevent them from being sucked in by the CbmHadronAnalysis WriteHistograms method // DeleteHistos(); if(fdMemoryTime>0.) CleanLHMemory(); } void CbmTofEventClusterizer::Finish(Double_t calMode) { if (fdEvent<100 ) return; // don't save histos with insufficient statistics SetCalMode(calMode); WriteHistos(); } /************************************************************************************/ // Functions common for all clusters approximations Bool_t CbmTofEventClusterizer::RegisterInputs() { FairRootManager *fManager = FairRootManager::Instance(); if( NULL == fManager) { LOG(error)<<"CbmTofEventClusterizer::RegisterInputs => Could not find FairRootManager!!!"; return kFALSE; } // if( NULL == fTofDigisColl) /* fTofPointsColl = (TClonesArray *) fManager->GetObject("TofPoint"); if( NULL == fTofPointsColl) { LOG(error)<<"CbmTofEventClusterizer::RegisterInputs => Could not get the TofPoint TClonesArray!!!"; return kFALSE; } // if( NULL == fTofPointsColl) fMcTracksColl = (TClonesArray *) fManager->GetObject("MCTrack"); if( NULL == fMcTracksColl) { LOG(error)<<"CbmTofEventClusterizer::RegisterInputs => Could not get the MCTrack TClonesArray!!!"; return kFALSE; } // if( NULL == fMcTracksColl) */ fEventsColl = dynamic_cast(fManager->GetObject("Event")); if(NULL == fEventsColl) fEventsColl = dynamic_cast(fManager->GetObject("CbmEvent")); if(NULL == fEventsColl) LOG(info) << "CbmEvent not found in input file, assume eventwise input" ; fTofRawDigisColl = (TClonesArray *) fManager->GetObject("CbmTofDigiExp"); if( NULL == fTofRawDigisColl) fTofRawDigisColl = (TClonesArray *) fManager->GetObject("CbmTofDigi"); if( NULL == fTofRawDigisColl) fTofRawDigisColl = (TClonesArray *) fManager->GetObject("TofDigiExp"); if( NULL == fTofRawDigisColl) fTofRawDigisColl = (TClonesArray *) fManager->GetObject("TofDigi"); if( NULL == fTofRawDigisColl) { LOG(error)<<"CbmTofEventClusterizer::RegisterInputs => Could not get the CbmTofDigi TClonesArray!!!"; return kFALSE; } // if( NULL == fTofDigisColl) fTrbHeader = (TTrbHeader *) fManager->GetObject("TofTrbHeader."); if( NULL == fTrbHeader) { LOG(info)<<"CbmTofEventClusterizer::RegisterInputs => Could not get TofTrbHeader Object"; } if (NULL == fEventsColl) { } else { // time based input LOG(info) << "CbmEvent found in input file, assume time based input" ; fTofDigisColl = new TClonesArray("CbmTofDigiExp"); } return kTRUE; } Bool_t CbmTofEventClusterizer::RegisterOutputs() { FairRootManager* rootMgr = FairRootManager::Instance(); // FairRunAna* ana = FairRunAna::Instance(); (VF) not used rootMgr->InitSink(); fTofCalDigisColl = new TClonesArray("CbmTofDigiExp"); fTofHitsColl = new TClonesArray("CbmTofHit"); fTofDigiMatchColl = new TClonesArray("CbmMatch",100); TString tHitBranchName; TString tHitDigiMatchBranchName; if(fbAlternativeBranchNames) { tHitBranchName = "ATofHit"; tHitDigiMatchBranchName = "ATofDigiMatch"; } else { tHitBranchName = "TofHit"; tHitDigiMatchBranchName = "TofDigiMatch"; } if (NULL == fEventsColl) { // Flag check to control whether digis are written in output root file rootMgr->Register( "TofCalDigi","Tof", fTofCalDigisColl, fbWriteDigisInOut); // Flag check to control whether digis are written in output root file rootMgr->Register(tHitBranchName,"Tof", fTofHitsColl, fbWriteHitsInOut); rootMgr->Register(tHitDigiMatchBranchName,"Tof", fTofDigiMatchColl, fbWriteHitsInOut); } else { // CbmEvent - mode fTofCalDigisCollOut = new TClonesArray("CbmTofDigiExp"); fTofHitsCollOut = new TClonesArray("CbmTofHit"); fTofDigiMatchCollOut = new TClonesArray("CbmMatch",100); rootMgr->Register( "TofCalDigi","Tof", fTofCalDigisCollOut, fbWriteDigisInOut); rootMgr->Register(tHitBranchName,"Tof", fTofHitsCollOut, fbWriteHitsInOut); rootMgr->Register(tHitDigiMatchBranchName,"Tof", fTofDigiMatchCollOut, fbWriteHitsInOut); } return kTRUE; } Bool_t CbmTofEventClusterizer::InitParameters() { // Initialize the TOF GeoHandler Bool_t isSimulation=kFALSE; LOG(info)<<"CbmTofEventClusterizer::InitParameters - Geometry, Mapping, ... ??" ; // Get Base Container FairRun* ana = FairRun::Instance(); FairRuntimeDb* rtdb=ana->GetRuntimeDb(); Int_t iGeoVersion = fGeoHandler->Init(isSimulation); if( k14a > iGeoVersion ) { LOG(error)<<"CbmTofEventClusterizer::InitParameters => Only compatible with geometries after v14a !!!" ; return kFALSE; } fTofId = new CbmTofDetectorId_v14a(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); if( 0 == fDigiPar ) { LOG(error)<<"CbmTofEventClusterizer::InitParameters => Could not obtain the CbmTofDigiPar "; return kFALSE; } fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); if( 0 == fDigiBdfPar ) { LOG(error)<<"CbmTofEventClusterizer::InitParameters => Could not obtain the CbmTofDigiBdfPar "; return kFALSE; } rtdb->initContainers( ana->GetRunId() ); LOG(info)<<"CbmTofEventClusterizer::InitParameter: currently " << fDigiPar->GetNrOfModules() << " digi cells " ; fdMaxTimeDist = fDigiBdfPar->GetMaxTimeDist(); // in ns fdMaxSpaceDist = fDigiBdfPar->GetMaxDistAlongCh(); // in cm if(fMaxTimeDist!=fdMaxTimeDist) { fdMaxTimeDist=fMaxTimeDist; // modify default fdMaxSpaceDist=fdMaxTimeDist*fDigiBdfPar->GetSignalSpeed()*0.5; // cut consistently on positions (with default signal velocity) } LOG(info)<<" BuildCluster with MaxTimeDist "< BeamRefType = "<GetNbDet(); (VF) not used Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes(); if (fTotMean !=0.) fdTTotMean=fTotMean; // adjust target mean for TOT fvCPTOff.resize( iNbSmTypes ); fvCPTotGain.resize( iNbSmTypes ); fvCPTotOff.resize( iNbSmTypes ); fvCPWalk.resize( iNbSmTypes ); fvCPDelTof.resize( iNbSmTypes ); for( Int_t iSmType = 0; iSmType < iNbSmTypes; iSmType++ ) { Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); fvCPTOff[iSmType].resize( iNbSm*iNbRpc ); fvCPTotGain[iSmType].resize( iNbSm*iNbRpc ); fvCPTotOff[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 ); fvCPTotOff[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 ); TProfile *hSvel = (TProfile *) gDirectory->FindObjectAny( Form("cl_SmT%01d_Svel", iSmType) ); // copy Histo to memory TDirectory * curdir = gDirectory; if(NULL != hSvel) { gDirectory->cd( oldir->GetPath() ); // TProfile *hSvelmem = (TProfile *)hSvel->Clone(); (VF) not used gDirectory->cd( curdir->GetPath() ); }else{ LOG(info)<<"Svel histogram not found for module type "<FindObjectAny( Form("cl_SmT%01d_Fpar%1d", iSmType, iPar) ); if(NULL != hFparcur) { gDirectory->cd( oldir->GetPath() ); // TProfile *hFparmem = (TProfile *)hFparcur->Clone(); (VF) not used gDirectory->cd( curdir->GetPath() ); } } for( Int_t iSm = 0; iSm < iNbSm; iSm++ ) for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { // update default parameter if(NULL != hSvel){ Double_t Vscal=1.;//hSvel->GetBinContent(iSm*iNbRpc+iRpc+1); if (Vscal==0.) Vscal=1.; fDigiBdfPar->SetSigVel(iSmType,iSm,iRpc,fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*Vscal); LOG(info)<<"Modify "<GetSigVel(iSmType,iSm,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)); TH1D *htempTot_Mean =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean",iSmType,iSm,iRpc)); TH1D *htempTot_Off =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off",iSmType,iSm,iRpc)); if(NULL != htempPos_pfx && NULL != htempTOff_pfx && NULL != htempTot_Mean && NULL != htempTot_Off) { Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); Int_t iNbinTot = htempTot_Mean->GetNbinsX(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { for(Int_t iSide=0; iSide<2; iSide++){ Double_t TotMean=htempTot_Mean->GetBinContent(iCh*2+1+iSide); //nh +1 empirical(?) if(0.001 < TotMean){ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] *= fdTTotMean / TotMean; } fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]=htempTot_Off->GetBinContent(iCh*2+1+iSide); } 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/fDigiBdfPar->GetSignalSpeed() ; Double_t dTYOff=YMean/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) ; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] += -dTYOff + TMean ; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] += +dTYOff + TMean ; /* if (iSmType==6 && iSm==0 && iRpc==1) { LOG(info) << "Skip loading other calib parameters for TSR "< " << Form(" %f, %f, %f, %f ",fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1], fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1]) <<", 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(debug)<<"Initialize Walk correction for " <GetNbinsX() != nbClWalkBinX) LOG(error)<<"CbmTofEventClusterizer::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_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); if (NULL==htmpDelTof) { LOG(debug)<<" Histos " << Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc, iSel) << " not found. " ; continue; } LOG(debug)<<" Load DelTof from histos "<< Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)<<"." ; 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_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); LOG(debug)<<" 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)<<"CbmTofEventClusterizer::InitCalibParameter: initialization done" ; return kTRUE; } /************************************************************************************/ Bool_t CbmTofEventClusterizer::LoadGeometry() { LOG(info)<<"CbmTofEventClusterizer::LoadGeometry starting for " <GetNbDet() << " described detectors, " <GetNrOfModules() << " geometrically known modules " ; Int_t iNrOfCells = fDigiPar->GetNrOfModules(); LOG(info)<<"Digi Parameter container contains "<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(debug) << "-I- InitPar "<GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ(),fNode) ; if(icell==0) { TGeoHMatrix* cMatrix = gGeoManager->GetCurrentMatrix(); fNode->Print(); cMatrix->Print(); } } Int_t iNbDet=fDigiBdfPar->GetNbDet(); fvDeadStrips.resize( iNbDet ); fvTimeLastDigi.resize( iNbDet ); for(Int_t iDetIndx=0; iDetIndxGetDetUId( iDetIndx ); Int_t iSmType = CbmTofAddress::GetSmType( iUniqueId ); Int_t iSmId = CbmTofAddress::GetSmId( iUniqueId ); Int_t iRpcId = CbmTofAddress::GetRpcId( iUniqueId ); LOG(info) << " DetIndx "< UniqueId "<GetSigVel(iSmType,iSmId,iRpcId), fvDeadStrips[iDetIndx]) ; Int_t iNbChan = fDigiBdfPar->GetNbChan( iSmType, iRpcId ); fvTimeLastDigi[iDetIndx].resize( iNbChan*2); for (Int_t iCh = 0; iChGetCell(iUCellId); if (NULL == fChannelInfo) break; LOG(debug3) << " Cell " << iCell << Form(" 0x%08x ",iUCellId) << Form(", fCh 0x%p ",fChannelInfo) << ", x: " << fChannelInfo->GetX() << ", y: " << fChannelInfo->GetY() << ", z: " << fChannelInfo->GetZ() ; } } // return kTRUE; Int_t iNbSmTypes = fDigiBdfPar->GetNbSmTypes(); if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { fStorDigiExp.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 ); fviClusterMul.resize( iNbSmTypes ); fvLastHits.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*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++ ) { fviClusterMul[iSmType].resize( iNbSm ); fvLastHits[iSmType].resize( iNbSm ); for( Int_t iRpc = 0; iRpc < iNbRpc; iRpc++ ) { fviClusterMul[iSmType][iSm].resize( iNbRpc ); fvLastHits[iSmType][iSm].resize( iNbRpc ); Int_t iNbChan = fDigiBdfPar->GetNbChan( iSmType, iRpc ); if(iNbChan == 0) { LOG(warning)<<"CbmTofEventClusterizer::LoadGeometry: StoreDigi without channels " << Form("SmTy %3d, Sm %3d, NbRpc %3d, Rpc, %3d ",iSmType,iSm,iNbRpc,iRpc) ; } LOG(debug1)<<"CbmTofEventClusterizer::LoadGeometry: StoreDigi with " << Form(" %3d %3d %3d %3d %5d ",iSmType,iSm,iNbRpc,iRpc,iNbChan) ; fStorDigiExp[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc].resize( iNbChan ); fvLastHits[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++ ) } // if( kTRUE == fDigiBdfPar->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() ) return kTRUE; } Bool_t CbmTofEventClusterizer::DeleteGeometry() { LOG(info)<<"CbmTofEventClusterizer::DeleteGeometry starting" ; return kTRUE; 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 CbmTofEventClusterizer::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( "TofEventClus_ClustBuildTime", "Time needed to build clusters in each event; Time [s]",4000, 0.0, 4.0 ); /* TH2* hTemp = 0;*/ // Sm related distributions fhSmCluPosition.resize( fDigiBdfPar->GetNbSmTypes() ); fhSmCluTOff.resize( fDigiBdfPar->GetNbSmTypes() ); fhSmCluSvel.resize( fDigiBdfPar->GetNbSmTypes() ); fhSmCluFpar.resize( fDigiBdfPar->GetNbSmTypes() ); fhTSmCluPosition.resize( fDigiBdfPar->GetNbSmTypes() ); fhTSmCluTOff.resize( fDigiBdfPar->GetNbSmTypes() ); fhTSmCluTRun.resize( fDigiBdfPar->GetNbSmTypes() ); for (Int_t iS=0; iS< fDigiBdfPar->GetNbSmTypes(); iS++){ Double_t YSCAL=50.; if (fPosYMaxScal !=0.) YSCAL=fPosYMaxScal; /* Int_t iUCellId = CbmTofAddress::GetUniqueAddress(0,0,0,0,iS); fChannelInfo = fDigiPar->GetCell(iUCellId); */ Int_t iUCellId(0); fChannelInfo = NULL; // Cover the case that the geometry file does not contain the module // indexed with 0 of a certain module type BUT modules with higher indices. for(Int_t iSM = 0; iSM < fDigiBdfPar->GetNbSm(iS); iSM++) { iUCellId = CbmTofAddress::GetUniqueAddress(iSM,0,0,0,iS); fChannelInfo = fDigiPar->GetCell(iUCellId); // Retrieve geometry information from the first module of a certain // module type that is found in the geometry file. if(NULL != fChannelInfo) { break; } } if(NULL == fChannelInfo){ LOG(warning)<<"No DigiPar for SmType " <GetSizey())*YSCAL; fhSmCluPosition[iS] = new TH2F( Form("cl_SmT%01d_Pos", iS), Form("Clu position of SmType %d; Sm+Rpc# []; ypos [cm]", iS ), fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 99, -YDMAX,YDMAX); Double_t TSumMax=1.E3; if (fTRefDifMax !=0.) TSumMax=fTRefDifMax; fhSmCluTOff[iS] = new TH2F( Form("cl_SmT%01d_TOff", iS), Form("Clu TimeZero in SmType %d; Sm+Rpc# []; TOff [ns]", iS ), fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 99, -TSumMax*2.,TSumMax*2. ); TProfile *hSvelcur = (TProfile *) gDirectory->FindObjectAny( Form("cl_SmT%01d_Svel", iS) ); if(NULL == hSvelcur) { fhSmCluSvel[iS] = new TProfile( Form("cl_SmT%01d_Svel", iS), Form("Clu Svel in SmType %d; Sm+Rpc# []; v/v_{nominal} ", iS ), fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0.8,1.2 ); fhSmCluSvel[iS]->Sumw2(); } else fhSmCluSvel[iS]=(TProfile *) hSvelcur->Clone(); fhSmCluFpar[iS].resize( 4 ); for (Int_t iPar =0; iPar<4; iPar++) { TProfile *hFparcur = (TProfile *) gDirectory->FindObjectAny( Form("cl_SmT%01d_Fpar%1d", iS, iPar) ); if(NULL == hFparcur){ LOG(info)<<"Histo "<< Form("cl_SmT%01d_Fpar%1d", iS, iPar) <<" not found, recreate ..."; fhSmCluFpar[iS][iPar] = new TProfile( Form("cl_SmT%01d_Fpar%1d", iS, iPar), Form("Clu Fpar %d in SmType %d; Sm+Rpc# []; value ", iPar, iS ), fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), -100.,100. ); } else fhSmCluFpar[iS][iPar]=(TProfile *) hFparcur->Clone(); } fhTSmCluPosition[iS].resize( iNSel ); fhTSmCluTOff[iS].resize( iNSel ); fhTSmCluTRun[iS].resize( iNSel ); for (Int_t iSel=0; iSelGetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 99, -YDMAX,YDMAX); fhTSmCluTOff[iS][iSel] = new TH2F( Form("cl_TSmT%01d_Sel%02d_TOff", iS, iSel), Form("Clu TimeZero in SmType %d under Selector %02d; Sm+Rpc# []; TOff [ns]", iS, iSel ), fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 0, fDigiBdfPar->GetNbSm(iS)*fDigiBdfPar->GetNbRpc(iS), 99, -TSumMax,TSumMax ); fhTSmCluTRun[iS][iSel] = new TH2F( Form("cl_TSmT%01d_Sel%02d_TRun", iS, iSel), Form("Clu TimeZero in SmType %d under Selector %02d; Event# []; TMean [ns]", iS, iSel ), 100, 0, MaxNbEvent, 99, -TSumMax,TSumMax ); } } // RPC related distributions Int_t iNbDet=fDigiBdfPar->GetNbDet(); LOG(info)<<" Define Clusterizer histos for "<GetDetUId( iDetIndx ); fDetIdIndexMap[iUniqueId]=iDetIndx; fviDetId[iDetIndx]=iUniqueId; Int_t iSmType = CbmTofAddress::GetSmType( iUniqueId ); Int_t iSmId = CbmTofAddress::GetSmId( iUniqueId ); Int_t iRpcId = CbmTofAddress::GetRpcId( iUniqueId ); Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSmId,iRpcId,0,0,iSmType); fChannelInfo = fDigiPar->GetCell(iUCellId); if (NULL==fChannelInfo) { LOG(warning)<<"No DigiPar for Det " < UniqueId "<GetSizex() <<", dy "<GetSizey() <GetNbChan( iSmType, 0 ) ; // check access to all channel infos for (Int_t iCh=0; iChGetNbChan( iSmType, iRpcId ); iCh++){ Int_t iCCellId = CbmTofAddress::GetUniqueAddress(iSmId,iRpcId,iCh,0,iSmType); fChannelInfo = fDigiPar->GetCell(iCCellId); if(NULL == fChannelInfo) LOG(warning)<GetCell(iUCellId); fhRpcDigiCor[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DigiCor", iSmType, iSmId, iRpcId ), Form("Digi Correlation of Rpc #%03d in Sm %03d of type %d; digi 0; digi 1", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId)); fhRpcDigiMul[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DigiMul", iSmType, iSmId, iRpcId ), Form("Digi Multiplicity of Rpc #%03d in Sm %03d of type %d; strip #; digi mul", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 15,0,15.); fhRpcDigiStatus[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DigiStatus", iSmType, iSmId, iRpcId ), Form("Digi Status of Rpc #%03d in Sm %03d of type %d; strip #; digi status", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 10,0,10.); fhRpcDigiDTLD[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DigiDTLD", iSmType, iSmId, iRpcId ), Form("Time distance to last digi of Rpc #%03d in Sm %03d of type %d; channel; t_{digi} - t_{previous digi} (s)", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId)*2,0,fDigiBdfPar->GetNbChan(iSmType,iRpcId)*2, 1000., 0., 5.); fhRpcCluMul[iDetIndx] = new TH1F( Form("cl_SmT%01d_sm%03d_rpc%03d_Mul", iSmType, iSmId, iRpcId ), Form("Clu multiplicity of Rpc #%03d in Sm %03d of type %d; M []; Cnts", iRpcId, iSmId, iSmType ), 2.+2.*fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,2.+2.*fDigiBdfPar->GetNbChan(iSmType,iRpcId)); fhRpcCluRate[iDetIndx] = new TH1D( Form("cl_SmT%01d_sm%03d_rpc%03d_rate", iSmType, iSmId, iRpcId ), Form("Clu rate of Rpc #%03d in Sm %03d of type %d; Time (s); Rate (Hz)", iRpcId, iSmId, iSmType ), 36000.,0.,3600.); fhRpcCluRate10s[iDetIndx] = new TH1D( Form("cl_SmT%01d_sm%03d_rpc%03d_rate10s", iSmType, iSmId, iRpcId ), Form(" Clu rate of Rpc #%03d in Sm %03d of type %d in last 10s; Time (s); Rate (Hz)", iRpcId, iSmId, iSmType ), 10000,0.,10.); fhRpcDTLastHits[iDetIndx] = new TH1F( Form("cl_SmT%01d_sm%03d_rpc%03d_DTLastHits", iSmType, iSmId, iRpcId ), Form("Clu #DeltaT to last hits of Rpc #%03d in Sm %03d of type %d; log( #DeltaT (ns)); counts", iRpcId, iSmId, iSmType ), 100.,0.,10.); fhRpcDTLastHits_Tot[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Tot_DTLH", iSmType, iSmId, iRpcId), Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); TOT [a.u.]", iRpcId, iSmId, iSmType), 100, 0., 10., 100, fdTOTMin, 4.*fdTOTMax); fhRpcDTLastHits_CluSize[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_CluSize_DTLH", iSmType, iSmId, iRpcId), Form("Clu Size of Rpc #%03d in Sm %03d of type %d; log(#DeltaT (ns)); CluSize []", iRpcId, iSmId, iSmType), 100, 0., 10., 16, 0.5, 16.5); Double_t YSCAL=50.; if (fPosYMaxScal !=0.) YSCAL=fPosYMaxScal; Double_t YDMAX=TMath::Max(2.,fChannelInfo->GetSizey())*YSCAL; fhRpcCluPosition[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Pos", iSmType, iSmId, iRpcId ), Form("Clu position of Rpc #%03d in Sm %03d of type %d; Strip []; ypos [cm]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -YDMAX,YDMAX); fhRpcCluPositionEvol[iDetIndx] = new TProfile( Form("cl_SmT%01d_sm%03d_rpc%03d_PosEvol", iSmType, iSmId, iRpcId ), Form("Clu position of Rpc #%03d in Sm %03d of type %d; Analysis Time [s]; ypos [cm]", iRpcId, iSmId, iSmType ), 1000,-1.,1.E4,-100.,100.); fhRpcCluTimeEvol[iDetIndx] = new TProfile( Form("cl_SmT%01d_sm%03d_rpc%03d_TimeEvol", iSmType, iSmId, iRpcId ), Form("Clu time of Rpc #%03d in Sm %03d of type %d; Analysis Time [s]; dT [ns]", iRpcId, iSmId, iSmType ), 1000,-1.,1.E4,-10.,10.); fhRpcCluDelPos[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DelPos", iSmType, iSmId, iRpcId ), Form("Clu position difference of Rpc #%03d in Sm %03d of type %d; Strip []; #Deltaypos(clu) [cm]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -10.,10.); fhRpcCluDelMatPos[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DelMatPos", iSmType, iSmId, iRpcId ), Form("Matched Clu position difference of Rpc #%03d in Sm %03d of type %d; Strip []; #Deltaypos(mat) [cm]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -5.,5.); Double_t TSumMax=1.E3; if (fTRefDifMax !=0.) TSumMax=fTRefDifMax; fhRpcCluTOff[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_TOff", iSmType, iSmId, iRpcId ), Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -TSumMax,TSumMax ); fhRpcCluDelTOff[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DelTOff", iSmType, iSmId, iRpcId ), Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip []; #DeltaTOff(clu) [ns]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -0.6,0.6); fhRpcCluDelMatTOff[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_DelMatTOff", iSmType, iSmId, iRpcId ), Form("Clu TimeZero Difference of Rpc #%03d in Sm %03d of type %d; Strip []; #DeltaTOff(mat) [ns]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, -0.6,0.6); fhRpcCluTrms[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Trms", iSmType, iSmId, iRpcId ), Form("Clu Time RMS of Rpc #%03d in Sm %03d of type %d; Strip []; Trms [ns]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 99, 0., 0.5 ); fhRpcCluTot[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Tot", iSmType, iSmId, iRpcId ), Form("Clu Tot of Rpc #%03d in Sm %03d of type %d; StripSide []; TOT [a.u.]", iRpcId, iSmId, iSmType ), 2*fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,2*fDigiBdfPar->GetNbChan(iSmType,iRpcId), 100, fdTOTMin, fdTOTMax); fhRpcCluSize[iDetIndx] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Size", iSmType, iSmId, iRpcId ), Form("Clu size of Rpc #%03d in Sm %03d of type %d; Strip []; size [strips]", iRpcId, iSmId, iSmType ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 16, 0.5, 16.5); // Walk histos fhRpcCluAvWalk[iDetIndx] = new TH2D( Form("cl_SmT%01d_sm%03d_rpc%03d_AvWalk", iSmType, iSmId, iRpcId), Form("Walk in SmT%01d_sm%03d_rpc%03d_AvWalk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId), nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinY,-TSumMax,TSumMax); fhRpcCluAvLnWalk[iDetIndx] = new TH2D( Form("cl_SmT%01d_sm%03d_rpc%03d_AvLnWalk", iSmType, iSmId, iRpcId), Form("Walk in SmT%01d_sm%03d_rpc%03d_AvLnWalk; log_{10}Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId), nbClWalkBinX,TMath::Log10(fdTOTMax/50.),TMath::Log10(fdTOTMax),nbClWalkBinY,-TSumMax,TSumMax); fhRpcCluWalk[iDetIndx].resize( fDigiBdfPar->GetNbChan(iSmType,iRpcId) ); for( Int_t iCh=0; iChGetNbChan(iSmType,iRpcId); iCh++){ fhRpcCluWalk[iDetIndx][iCh].resize( 2 ); for (Int_t iSide=0; iSide<2; iSide++) { fhRpcCluWalk[iDetIndx][iCh][iSide]= new TH2D( Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk", iSmType, iSmId, iRpcId, iCh, iSide ), Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Walk; Tot [a.u.]; #DeltaT [ns]", iSmType, iSmId, iRpcId, iCh, iSide ), nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinY,-TSumMax,TSumMax ); } /* (fhRpcCluWalk[iDetIndx]).push_back( hTemp ); */ } } // Clusterizing monitoring LOG(info)<<" Define Clusterizer monitoring histos "; fhCluMulCorDutSel = new TH2F(Form("hCluMulCorDutSel"), Form("Cluster Multiplicity Correlation of Dut %d%d%d with Sel %d%d%d; Mul(Sel); Mul(Dut)", fDutId, fDutSm, fDutRpc, fSelId, fSelSm, fSelRpc ), 32,0,32,32,0,32); fhEvCluMul = new TH2F(Form("hEvCluMul"), Form("Time Evolution of Cluster Multiplicity in Event; Time (s); Mul"), 1800,0,1800,100,0,100); fhDigSpacDifClust = new TH1I( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_NbDigiPerChan", "Nb of Digis per channel; Nb Digis []", 100, 0.0, 100.0 ); // Trigger selected histograms if (0GetDetUId( iDetIndx ); Int_t iSmType = CbmTofAddress::GetSmType( iUniqueId ); Int_t iSmId = CbmTofAddress::GetSmId( iUniqueId ); Int_t iRpcId = CbmTofAddress::GetRpcId( iUniqueId ); Int_t iUCellId = CbmTofAddress::GetUniqueAddress(iSmId,iRpcId,0,0,iSmType); fChannelInfo = fDigiPar->GetCell(iUCellId); if (NULL==fChannelInfo) { LOG(warning)<<"No DigiPar for Det " < UniqueId "<GetSizex() <<", dy "<GetSizey() <GetNbChan(iSmType,iRpcId) ; fhTRpcCluMul[iDetIndx].resize( iNSel ); fhTRpcCluPosition[iDetIndx].resize( iNSel ); fhTRpcCluTOff[iDetIndx].resize( iNSel ); fhTRpcCluTofOff[iDetIndx].resize( iNSel ); fhTRpcCluTot[iDetIndx].resize( iNSel ); fhTRpcCluSize[iDetIndx].resize( iNSel ); fhTRpcCluAvWalk[iDetIndx].resize( iNSel ); fhTRpcCluDelTof[iDetIndx].resize( iNSel ); fhTRpcCludXdY[iDetIndx].resize( iNSel ); fhTRpcCluWalk[iDetIndx].resize( iNSel ); fhTRpcCluWalk2[iDetIndx].resize( iNSel ); fhTRpcCluTOffDTLastHits[iDetIndx].resize( iNSel ); fhTRpcCluTotDTLastHits[iDetIndx].resize( iNSel ); fhTRpcCluSizeDTLastHits[iDetIndx].resize( iNSel ); fhTRpcCluMemMulDTLastHits[iDetIndx].resize( iNSel ); for (Int_t iSel=0; iSelGetNbChan(iSmType,iRpcId)+2,0.,fDigiBdfPar->GetNbChan(iSmType,iRpcId)+2); if (NULL == fhTRpcCluMul[iDetIndx][iSel]) LOG(fatal)<<" Histo not generated !"; Double_t YSCAL=50.; if (fPosYMaxScal !=0.) YSCAL=fPosYMaxScal; Double_t YDMAX=TMath::Max(2.,fChannelInfo->GetSizey())*YSCAL; fhTRpcCluPosition[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Pos", iSmType, iSmId, iRpcId, iSel ), Form("Clu position of Rpc #%03d in Sm %03d of type %d under Selector %02d; Strip []; ypos [cm]", iRpcId, iSmId, iSmType, iSel ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 100, -YDMAX,YDMAX ); Double_t TSumMax=1.E4; if (fTRefDifMax !=0.) TSumMax=fTRefDifMax; fhTRpcCluTOff[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff", iSmType, iSmId, iRpcId, iSel ), Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d under Selector %02d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType, iSel ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 999, -TSumMax,TSumMax ); fhTRpcCluTofOff[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TofOff", iSmType, iSmId, iRpcId, iSel ), Form("Clu TimeDeviation of Rpc #%03d in Sm %03d of type %d under Selector %02d; Strip []; TOff [ns]", iRpcId, iSmId, iSmType, iSel ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 999, -TSumMax*4.,TSumMax*4.); if (fTotMax !=0.) fdTOTMax=fTotMax; fhTRpcCluTot[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot", iSmType, iSmId, iRpcId, iSel ), Form("Clu Tot of Rpc #%03d in Sm %03d of type %d under Selector %02d; StripSide []; TOT [a.u.]", iRpcId, iSmId, iSmType, iSel ), fDigiBdfPar->GetNbChan(iSmType,iRpcId)*2, 0, fDigiBdfPar->GetNbChan(iSmType,iRpcId)*2, 100, fdTOTMin, fdTOTMax); fhTRpcCluSize[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size", iSmType, iSmId, iRpcId, iSel ), Form("Clu size of Rpc #%03d in Sm %03d of type %d under Selector %02d; Strip []; size [strips]", iRpcId, iSmId, iSmType, iSel ), fDigiBdfPar->GetNbChan(iSmType,iRpcId),0,fDigiBdfPar->GetNbChan(iSmType,iRpcId), 16, 0.5, 16.5); // Walk histos fhTRpcCluAvWalk[iDetIndx][iSel] = new TH2D( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk", iSmType, iSmId, iRpcId, iSel), Form("Walk in SmT%01d_sm%03d_rpc%03d_Sel%02d_AvWalk; TOT; T-TSel", iSmType, iSmId, iRpcId, iSel), nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinY,-TSumMax,TSumMax); // Tof Histos fhTRpcCluDelTof[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSmId, iRpcId, iSel), Form("SmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof; TRef-TSel; T-TSel", iSmType, iSmId, iRpcId, iSel), nbClDelTofBinX,-fdDelTofMax,fdDelTofMax,nbClDelTofBinY,-TSumMax,TSumMax); // Position deviation histos fhTRpcCludXdY[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY", iSmType, iSmId, iRpcId, iSel), Form("SmT%01d_sm%03d_rpc%03d_Sel%02d_dXdY; #Delta x [cm]; #Delta y [cm];", iSmType, iSmId, iRpcId, iSel), nbCldXdYBinX,-dXdYMax,dXdYMax,nbCldXdYBinY,-dXdYMax,dXdYMax); fhTRpcCluWalk2[iDetIndx][iSel]= new TH3F ( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Walk2", iSmType, iSmId, iRpcId, iSel ), Form("Walk in SmT%01d_sm%03d_rpc%03d_Sel%02d_Walk2", iSmType, iSmId, iRpcId, iSel ), nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinY,-TSumMax,TSumMax); fhTRpcCluWalk[iDetIndx][iSel].resize( fDigiBdfPar->GetNbChan(iSmType,iRpcId) ); for( Int_t iCh=0; iChGetNbChan(iSmType,iRpcId); iCh++){ fhTRpcCluWalk[iDetIndx][iSel][iCh].resize( 2 ); for (Int_t iSide=0; iSide<2; iSide++) { fhTRpcCluWalk[iDetIndx][iSel][iCh][iSide]= new TH2D( Form("cl_SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk", iSmType, iSmId, iRpcId, iCh, iSide, iSel ), Form("Walk in SmT%01d_sm%03d_rpc%03d_Ch%03d_S%01d_Sel%02d_Walk", iSmType, iSmId, iRpcId, iCh, iSide, iSel ), nbClWalkBinX,fdTOTMin,fdTOTMax,nbClWalkBinY,-TSumMax,TSumMax); } } fhTRpcCluTOffDTLastHits[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_TOff_DTLH", iSmType, iSmId, iRpcId, iSel ), Form("Clu TimeZero of Rpc #%03d in Sm %03d of type %d under Selector %02d; log(#DeltaT (ns)); TOff [ns]", iRpcId, iSmId, iSmType, iSel ), 100, 0., 10., 99, -TSumMax,TSumMax ); fhTRpcCluTotDTLastHits[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Tot_DTLH", iSmType, iSmId, iRpcId, iSel ), Form("Clu Tot of Rpc #%03d in Sm %03d of type %d under Selector %02d; log(#DeltaT (ns)); TOT [a.u.]", iRpcId, iSmId, iSmType, iSel ), 100, 0., 10., 100, fdTOTMin, fdTOTMax); fhTRpcCluSizeDTLastHits[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_Size_DTLH", iSmType, iSmId, iRpcId, iSel ), Form("Clu size of Rpc #%03d in Sm %03d of type %d under Selector %02d; log(#DeltaT (ns)); size [strips]", iRpcId, iSmId, iSmType, iSel ), 100, 0., 10., 10, 0.5, 10.5); fhTRpcCluMemMulDTLastHits[iDetIndx][iSel] = new TH2F( Form("cl_SmT%01d_sm%03d_rpc%03d_Sel%02d_MemMul_DTLH", iSmType, iSmId, iRpcId, iSel ), Form("Clu Memorized Multiplicity of Rpc #%03d in Sm %03d of type %d under Selector %02d; log(#DeltaT (ns)); TOff [ns]", iRpcId, iSmId, iSmType, iSel ), 100, 0., 10., 10, 0, 10 ); } } } // MC reference fhHitsPerTracks = new TH1I( "Clus_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( "Clus_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( "Clus_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( "Clus_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( "Clus_TofTimePtVsHit", "Time resolution for TofHits containing Digis from a single MC Track; t(1st Mc Point) -tTofHit [ns]", 2000, 0.0, 50.0, 2000, 0.0, 50.0 ); } else { fhTimeResSingHits = new TH1I( "Clus_TofTimeResClust", "Time resolution for TofHits containing Digis from a single TofPoint; tMcPoint -tTofHit [ns]", 10000, -25.0, 25.0 ); fhTimeResSingHitsB = new TH2I( "Clus_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( "Clus_TofTimePtVsHit", "Time resolution for TofHits containing Digis from a single TofPoint; tMcPoint -tTofHit [ps]", 2000, 0.0, 50.0, 2000, 0.0, 50.0 ); } // else of if( kTRUE == fDigiBdfPar->ClustUseTrackId() ) fhClusterSize = new TH1I( "Clus_ClusterSize", "Cluster Size distribution; Cluster Size [Strips]", 100, 0.5, 100.5 ); fhClusterSizeType = new TH2I( "Clus_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( "Clus_TrackMul", "Number of MC tracks generating the cluster; MC Tracks multiplicity []", 100, 0.5, 100.5 ); fhClusterSizeMulti = new TH2I( "Clus_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( "Clus_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( "Clus_HiTrkMulPos", "Position of Clusters with >1 MC tracks generating the cluster; X [cm]; Y [cm]", 1500, -750, 750, 1000, -500, 500 ); fhAllTrkMulPos = new TH2D( "Clus_AllTrkMulPos", "Position of all clusters generating the cluster; X [cm]; Y [cm]", 1500, -750, 750, 1000, -500, 500 ); fhMultiTrkProbPos = new TH2D( "Clus_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() ) gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager! return kTRUE; } Bool_t CbmTofEventClusterizer::FillHistos() { fhClustBuildTime->Fill( fStop.GetSec() - fStart.GetSec() + (fStop.GetNanoSec() - fStart.GetNanoSec())/1e9 ); Int_t iNbTofHits = fTofHitsColl->GetEntries(); CbmTofHit *pHit; //gGeoManager->SetTopVolume( gGeoManager->FindVolumeFast("tof_v14a") ); gGeoManager->CdTop(); if(0GetNbDet(); for(Int_t iDetIndx=0; iDetIndx 0) iDetMul++; dCluMul += fviClusterMul[iSmType][iSm][iRpc]; // total hit multiplicity in event fhRpcCluMul[iDetIndx]->Fill(fviClusterMul[iSmType][iSm][iRpc]); // } } fhCluMulCorDutSel->Fill(fviClusterMul[fSelId][fSelSm][fSelRpc],fviClusterMul[fDutId][fDutSm][fDutRpc]); // do input distributions first for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; if (StartAnalysisTime == 0.) { StartAnalysisTime = pHit->GetTime(); LOG(info) << "StartAnalysisTime set to "<GetAddress() & DetMask); std::map::iterator it=fDetIdIndexMap.find(iDetId); if (it == fDetIdIndexMap.end()) continue; // continue for invalid detector index Int_t iDetIndx=it->second; //fDetIdIndexMap[iDetId]; Int_t iSmType = CbmTofAddress::GetSmType( iDetId ); Int_t iSm = CbmTofAddress::GetSmId( iDetId ); Int_t iRpc = CbmTofAddress::GetRpcId( iDetId ); Int_t iCh = CbmTofAddress::GetChannelId( pHit->GetAddress() ); Double_t dTimeAna=(pHit->GetTime() - StartAnalysisTime)/1.E9; //LOG(debug)<<"TimeAna "<GetTime()<<", "<Fill(dTimeAna,1./fhRpcCluRate[iDetIndx]->GetBinWidth(1)); // deal with spill structure Double_t dTimeAna10s = pHit->GetTime() - fdStartAna10s; if(iHitInd == 0) { fhEvCluMul->Fill(dTimeAna,dCluMul); if ( iDetMul > 5 && fviClusterMul[fiBeamRefType][fiBeamRefSm][fiBeamRefDet]>0) { // FIXME: hardwired constants if(dTLEvt == 0) dTLEvt=pHit->GetTime(); Double_t dDTLEvt=pHit->GetTime()-dTLEvt; dTLEvt=pHit->GetTime(); if( dDTLEvt> fdSpillBreak*1.E9 && dTimeAna10s > fdSpillDuration*1.E9 ) { //if( dDTLEvt> 5.E8 && dTimeAna10s > 10.) { fdStartAna10s=pHit->GetTime(); iNSpill ++; LOG(info)<<"Resetting 10s rate histo for spill " <GetNbDet(); iDet++) { if ( NULL == fhRpcCluRate10s[iDet]) continue; fhRpcCluRate10s[iDet]->Reset("ICES"); } dTimeAna10s=0.; } } } if(fdStartAna10s>0.) fhRpcCluRate10s[iDetIndx]->Fill(dTimeAna10s/1.E9,1./fhRpcCluRate10s[iDetIndx]->GetBinWidth(1)); if(fdMemoryTime>0. && fvLastHits[iSmType][iSm][iRpc][iCh].size()==0) LOG(fatal)< hit not stored in memory for TSRC %d%d%d%d", iSmType,iSm,iRpc,iCh) ; //CheckLHMemory(); if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1) { // check for outdated hits //std::list::iterator it0=fvLastHits[iSmType][iSm][iRpc][iCh].begin(); //std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); //CbmTofHit* pH0 = *it0; //CbmTofHit* pHL = *(--itL); CbmTofHit* pH0 = fvLastHits[iSmType][iSm][iRpc][iCh].front(); CbmTofHit* pHL = fvLastHits[iSmType][iSm][iRpc][iCh].back(); if(pH0->GetTime() > pHL->GetTime()) LOG(warning)<GetTime() - pH0->GetTime() ) ; // while( (*((std::list::iterator) fvLastHits[iSmType][iSm][iRpc][iCh].begin()))->GetTime()+fdMemoryTime < pHit->GetTime() while( fvLastHits[iSmType][iSm][iRpc][iCh].size() > 2. || fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()+fdMemoryTime < pHit->GetTime() ) { LOG(debug) << " pop from list size "<< fvLastHits[iSmType][iSm][iRpc][iCh].size() <GetTime()- fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()) //(*((std::list::iterator) fvLastHits[iSmType][iSm][iRpc][iCh].begin()))->GetTime()) ; if( fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != pHit->GetAddress() ) LOG(fatal)<GetAddress(), fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()) ; fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete(); fvLastHits[iSmType][iSm][iRpc][iCh].pop_front(); } } //fvLastHits[iSmType][iSm][iRpc][iCh].size()>1) // plot remaining time difference to previous hits if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1){ // check for previous hits in memory time interval CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchColl->At( iHitInd ); Double_t dTotSum=0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink+=2) { // loop over digis CbmLink L0 = digiMatch->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); Int_t iDigInd1=(digiMatch->GetLink(iLink+1)).GetIndex(); CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); CbmTofDigiExp *pDig1 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd1)); dTotSum += pDig0->GetTot() + pDig1->GetTot(); } std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); itL--; for(size_t iH=0; iHFill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime())); fhRpcDTLastHits_CluSize[iDetIndx]->Fill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime()), digiMatch->GetNofLinks()/2.); fhRpcDTLastHits_Tot[iDetIndx]->Fill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime()),dTotSum); } } } // iHitInd loop end // do reference first dTRef = dDoubleMax; fTRefHits=0; for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if (NULL==pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); if( fiBeamRefAddr == iDetId ){ if(fviClusterMul[fiBeamRefType][fiBeamRefSm][fiBeamRefDet]>fiBeamRefMulMax) break; // Check Tot CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchColl->At(iHitInd); Double_t TotSum=0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink+=2){ // loop over digis CbmLink L0 = digiMatch->GetLink(iLink); //vDigish.at(ivDigInd); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofCalDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); TotSum += pDig0->GetTot(); } } TotSum /= (0.5 * digiMatch->GetNofLinks()); if( TotSum > fhRpcCluTot[iIndexDut]->GetYaxis()->GetXmax()) continue; // ignore too large clusters fTRefHits=1; if(pHit->GetTime() < dTRef) { dTRef = pHit->GetTime(); pBeamRef=pHit; } iBeamRefMul++; }else{ //additional reference type multiplicity if(fiBeamRefType == CbmTofAddress::GetSmType( iDetId ) ) iBeamAddRefMul++; } } LOG(debug) <<"CbmTofEventClusterizer::FillHistos: BRefMul: " <-1) {iR0=fDutRpc; iRl=fDutRpc+1;} for(Int_t iRpc=iR0; iRpc0 && iDutMul < fiCluMulMax ) { LOG(debug1)<<"Found selector 0, NbHits "<At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); LOG(debug1)<GetTime(), dTTrig[iSel]); //if( fDutId == CbmTofAddress::GetSmType( iDetId )) if( fDutAddr == iDetId) { if(pHit->GetTime() < dTTrig[iSel]) { if(TMath::Abs(pBeamRef->GetTime()-pHit->GetTime()) < fdDelTofMax) { // < fhTRpcCluTOff[iIndexDut][iSel]->GetYaxis()->GetXmax()) { dTTrig[iSel] = pHit->GetTime(); pTrig[iSel] = pHit; BSel[iSel]=kTRUE; } } } } if(BSel[iSel]) LOG(debug)<GetAddress(),dTTrig[iSel]); } break; case 1 : // MRef & BRef iRl=fviClusterMul[fSelId][fSelSm].size(); if(fIdMode==0 && fSelRpc>-1) {iR0=fSelRpc; iRl=fSelRpc+1;} for(Int_t iRpc=iR0; iRpc0 && iRefMul < fiCluMulMax ) { LOG(debug1)<<"CbmTofEventClusterizer::FillHistos(): Found selector 1"; dTTrig[iSel]=dDoubleMax; for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); if( fSelAddr == iDetId ) { if(pHit->GetTime() < dTTrig[iSel]) { if(TMath::Abs(pBeamRef->GetTime()-pHit->GetTime()) < fdDelTofMax) { // < fhTRpcCluTOff[iIndexDut][iSel]->GetYaxis()->GetXmax()) { dTTrig[iSel] = pHit->GetTime(); pTrig[iSel] = pHit; BSel[iSel]=kTRUE; } } } } if(BSel[iSel]) LOG(debug)<GetAddress(),dTTrig[iSel]); } break; default : LOG(info)<<"CbmTofEventClusterizer::FillHistos: selection not implemented "<10){ dTTrig[iSel]=dTRef; Int_t iReqMul=(fTRefMode-fTRefMode%10)/10; if ( iDetMul < iReqMul ) BSel[iSel]=0; } } // iSel - loop end if (fSel2Id > -1 ) { // confirm selector by independent match for (Int_t iSel=0; iSel 0 && fviClusterMul[fSel2Id][fSel2Sm][fSel2Rpc] < fiCluMulMax) for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); if( fSel2Addr == iDetId ) { Double_t dzscal=1.; if(fEnableMatchPosScaling) dzscal=pHit->GetZ()/pTrig[iSel]->GetZ(); Double_t dSEl2dXdz=(pHit->GetX()-pTrig[iSel]->GetX())/(pHit->GetZ()-pTrig[iSel]->GetZ()); Double_t dSEl2dYdz=(pHit->GetY()-pTrig[iSel]->GetY())/(pHit->GetZ()-pTrig[iSel]->GetZ()); if ( iDetId == fiBeamRefAddr || ( TMath::Sqrt(TMath::Power(pHit->GetX()-dzscal*pTrig[iSel]->GetX(),2.) +TMath::Power(pHit->GetY()-dzscal*pTrig[iSel]->GetY(),2.) )At( iHitInd ); if (NULL==pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); if( fiBeamRefType == CbmTofAddress::GetSmType( iDetId )){ if(fiBeamRefSm == CbmTofAddress::GetSmId( iDetId )) { Double_t dDT2=0.; Double_t dNT=0.; for (Int_t iSel=0; iSelGetTime()-dTTrig[iSel],2); dNT++; } } if( dNT > 0) if( dDT2/dNT < dRefChi2 ) { fTRefHits=1; dTRef = pHit->GetTime(); dRefChi2 = dDT2/dNT; } } } } */ UInt_t uTriggerPattern=1; if(NULL != fTrbHeader) uTriggerPattern=fTrbHeader->GetTriggerPattern(); else { for (Int_t iSel=0; iSelGetAddress() & DetMask ) ) ); } } for (Int_t iSel=0; iSel0) { for(UInt_t uChannel = 0; uChannel < 16; uChannel++) { if( uTriggerPattern & (0x1 << uChannel) ) { fhSeldT[iSel]->Fill(dTTrig[iSel]-dTRef,uChannel); } } } } } } // 0At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & DetMask); std::map::iterator it=fDetIdIndexMap.find(iDetId); if (it == fDetIdIndexMap.end()) continue; // continue for invalid detector index Int_t iDetIndx=it->second; //fDetIdIndexMap[iDetId]; Int_t iSmType = CbmTofAddress::GetSmType( iDetId ); Int_t iSm = CbmTofAddress::GetSmId( iDetId ); Int_t iRpc = CbmTofAddress::GetRpcId( iDetId ); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType ); if(-1Fill(fviClusterMul[iSmType][iSm][iRpc],w); } } if(fviClusterMul[iSmType][iSm][iRpc] > fiCluMulMax) continue; // skip this event if(iBeamRefMul == 0) break; Int_t iChId = pHit->GetAddress(); fChannelInfo = fDigiPar->GetCell( iChId ); Int_t iCh = CbmTofAddress::GetChannelId(iChId); if(NULL == fChannelInfo){ LOG(error) << "Invalid Channel Pointer for ChId " << Form(" 0x%08x ",iChId)<<", Ch "<local trafo gGeoManager->FindNode(fChannelInfo->GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ()); LOG(debug1)<<"Hit info: " <GetX(),pHit->GetY(),pHit->GetTime(),fChannelInfo->GetX(),fChannelInfo->GetY(), iHitInd ) ; Double_t hitpos[3]; hitpos[0]=pHit->GetX(); hitpos[1]=pHit->GetY(); hitpos[2]=pHit->GetZ(); Double_t hitpos_local[3]; TGeoNode* cNode= gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos, hitpos_local); LOG(debug1)<< Form(" MasterToLocal for %d, %d%d%d, node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)", iDetIndx, iSmType, iSm, iRpc, cNode, hitpos[0], hitpos[1], hitpos[2], hitpos_local[0], hitpos_local[1], hitpos_local[2]) ; fhRpcCluPosition[iDetIndx]->Fill((Double_t)iCh,hitpos_local[1]); //pHit->GetY()-fChannelInfo->GetY()); fhSmCluPosition[iSmType]->Fill((Double_t)(iSm*iNbRpc+iRpc),hitpos_local[1]); for (Int_t iSel=0; iSelFill((Double_t)iCh,hitpos_local[1]); //pHit->GetY()-fChannelInfo->GetY()); fhTSmCluPosition[iSmType][iSel]->Fill((Double_t)(iSm*iNbRpc+iRpc),hitpos_local[1]); } if(TMath::Abs(hitpos_local[1])>fChannelInfo->GetSizey()*fPosYMaxScal) continue; Double_t dTimeAna=(pHit->GetTime() - StartAnalysisTime)/1.E9; if(dTRef !=0.) fhRpcCluTimeEvol[iDetIndx]->Fill(dTimeAna,pHit->GetTime()-dTRef); fhRpcCluPositionEvol[iDetIndx]->Fill(dTimeAna,hitpos_local[1]); //LOG(info) << "Fill TEvol at " << dTimeAna ; LOG(debug1)<<" TofDigiMatchColl entries:" <GetEntries() ; if(iHitInd>fTofDigiMatchColl->GetEntries()){ LOG(error)<<" Inconsistent DigiMatches for Hitind " <GetEntries() ; } CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchColl->At( iHitInd ); LOG(debug1)<<" got " <GetNofLinks()<< " matches for iCh "<Fill((Double_t)iCh,digiMatch->GetNofLinks()/2.); for (Int_t iSel=0; iSelFill((Double_t)iCh,digiMatch->GetNofLinks()/2.); if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1){ // check for previous hits in memory time interval std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); itL--; for(size_t iH=0; iHFill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime()), digiMatch->GetNofLinks()/2.); } } } Double_t TotSum=0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink++){ // loop over digis CbmLink L0 = digiMatch->GetLink(iLink); //vDigish.at(ivDigInd); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofCalDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); TotSum += pDig0->GetTot(); } } Double_t dMeanTimeSquared=0.; Double_t dNstrips=0.; Double_t dDelTof=0.; Double_t dTcor[iNSel]; Double_t dTTcor[iNSel]; Double_t dZsign[iNSel]; Double_t dzscal=1.; //Double_t dDist=0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink+=2) { // loop over digis CbmLink L0 = digiMatch->GetLink(iLink); //vDigish.at(ivDigInd); Int_t iDigInd0=L0.GetIndex(); Int_t iDigInd1=(digiMatch->GetLink(iLink+1)).GetIndex(); //vDigish.at(ivDigInd+1); //LOG(debug1)<<" " << iDigInd0<<", "<GetEntries() && iDigInd1 < fTofCalDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); CbmTofDigiExp *pDig1 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd1)); if((Int_t)pDig0->GetType()!=iSmType){ LOG(error)<GetChannel()<<", Side "<GetSide() <<", StripSide "<<(Double_t)iCh*2.+pDig0->GetSide() <<" Digi 1 "<GetChannel()<<", Side "<GetSide() <<", StripSide "<<(Double_t)iCh*2.+pDig1->GetSide() <<", Tot0 " << pDig0->GetTot() <<", Tot1 "<GetTot(); fhRpcCluTot[iDetIndx]->Fill(pDig0->GetChannel()*2.+pDig0->GetSide(),pDig0->GetTot()); fhRpcCluTot[iDetIndx]->Fill(pDig1->GetChannel()*2.+pDig1->GetSide(),pDig1->GetTot()); Int_t iCh0=pDig0->GetChannel(); Int_t iCh1=pDig1->GetChannel(); Int_t iS0=pDig0->GetSide(); Int_t iS1=pDig1->GetSide(); if(iCh0 != iCh1 || iS0==iS1){ LOG(fatal)<GetTime()) <GetTime()) ; continue; } if(0>iCh0 || fDigiBdfPar->GetNbChan( iSmType, iRpc )<=iCh0) { LOG(error)<GetNofLinks()>2 ) //&& digiMatch->GetNofLinks()<8 ) // FIXME: hardcoded limits on CluSize { dNstrips+=1.; dMeanTimeSquared += TMath::Power(0.5*(pDig0->GetTime()+pDig1->GetTime())-pHit->GetTime(),2); // fhRpcCluAvWalk[iDetIndx]->Fill(0.5*(pDig0->GetTot()+pDig1->GetTot()), // 0.5*(pDig0->GetTime()+pDig1->GetTime())-pHit->GetTime()); fhRpcCluAvLnWalk[iDetIndx]->Fill(TMath::Log10(0.5*(pDig0->GetTot()+pDig1->GetTot())), 0.5*(pDig0->GetTime()+pDig1->GetTime())-pHit->GetTime()); Double_t dTotWeigth=(pDig0->GetTot()+pDig1->GetTot())/TotSum; Double_t dCorWeigth=1.-dTotWeigth; fhRpcCluDelTOff[iDetIndx]->Fill(pDig0->GetChannel(), dCorWeigth*(0.5*(pDig0->GetTime()+pDig1->GetTime())-pHit->GetTime())); Double_t dDelPos=0.5*(pDig0->GetTime()-pDig1->GetTime())*fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc); if(0==pDig0->GetSide()) dDelPos *= -1.; fhRpcCluDelPos[iDetIndx]->Fill(pDig0->GetChannel(),dCorWeigth*(dDelPos-hitpos_local[1])); fhRpcCluWalk[iDetIndx][iCh0][iS0]->Fill(pDig0->GetTot(),pDig0->GetTime()-(pHit->GetTime() -(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))); fhRpcCluWalk[iDetIndx][iCh1][iS1]->Fill(pDig1->GetTot(),pDig1->GetTime()-(pHit->GetTime() -(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))); fhRpcCluAvWalk[iDetIndx]->Fill(pDig0->GetTot(),pDig0->GetTime()-(pHit->GetTime() -(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))); fhRpcCluAvWalk[iDetIndx]->Fill(pDig1->GetTot(),pDig1->GetTime()-(pHit->GetTime() -(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))); } // end of Clustersize > 1 condition LOG(debug1)<<" fhTRpcCluTot: Digi 0 "<GetChannel()<<", Side "<GetSide() <<", StripSide "<<(Double_t)iCh*2.+pDig0->GetSide() <<" Digi 1 "<GetChannel()<<", Side "<GetSide() <<", StripSide "<<(Double_t)iCh*2.+pDig1->GetSide(); for (Int_t iSel=0; iSelGetAddress() == pTrig[iSel]->GetAddress() ) continue; fhTRpcCluTot[iDetIndx][iSel]->Fill(pDig0->GetChannel()*2.+pDig0->GetSide(),pDig0->GetTot()); fhTRpcCluTot[iDetIndx][iSel]->Fill(pDig1->GetChannel()*2.+pDig1->GetSide(),pDig1->GetTot()); if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1){ // check for previous hits in memory time interval std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); itL--; for(size_t iH=0; iHFill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime()), pDig0->GetTot()); fhTRpcCluTotDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(pHit->GetTime()-(*itL)->GetTime()), pDig1->GetTot()); } } if(iLink==0) { // Fill histo only once (for 1. digi entry) if(fEnableMatchPosScaling) dzscal=pHit->GetZ()/pTrig[iSel]->GetZ(); fhTRpcCludXdY[iDetIndx][iSel]->Fill(pHit->GetX()-dzscal*pTrig[iSel]->GetX(), pHit->GetY()-dzscal*pTrig[iSel]->GetY()); dZsign[iSel] = 1.; if(pHit->GetZ() < pTrig[iSel]->GetZ()) dZsign[iSel]=-1.; } //// look for geometrical match with selector hit if( iSmType==fiBeamRefType // to get entries in diamond/BeamRef histos //if( iSmType == 5 // FIXME, to get entries in diamond histos || TMath::Sqrt(TMath::Power(pHit->GetX()-dzscal*pTrig[iSel]->GetX(),2.) +TMath::Power(pHit->GetY()-dzscal*pTrig[iSel]->GetY(),2.))GetX()-(pTrig[iSel]->GetX()+ddXdZ[iSel]*(pHit->GetZ()-(pTrig[iSel]->GetZ()))),2.) +TMath::Power(pHit->GetY()-(pTrig[iSel]->GetY()+ddYdZ[iSel]*(pHit->GetZ()-(pTrig[iSel]->GetZ()))),2.)) > 0.5*fdCaldXdYMax) continue; // refine position selection cut in cosmic measurement dTcor[iSel]=0.; // precaution if (dTRef !=0. && TMath::Abs(dTRef-dTTrig[iSel])GetX()-pTrig[iSel]->GetX(),2.) +TMath::Power(pHit->GetY()-pTrig[iSel]->GetY(),2.) +TMath::Power(pHit->GetZ()-pTrig[iSel]->GetZ(),2.)); */ // determine correction value //if(fiBeamRefAddr != iDetId) // do not do this for reference counter itself if(fTRefMode<11) // do not do this for trigger counter itself { Double_t dTentry=dTRef-dTTrig[iSel]+fdDelTofMax; Int_t iBx = dTentry/2./fdDelTofMax*nbClDelTofBinX; if(iBx<0) iBx=0; if(iBx>nbClDelTofBinX-1) iBx=nbClDelTofBinX-1; Double_t dBinWidth=2.*fdDelTofMax/nbClDelTofBinX; Double_t dDTentry=dTentry-((Double_t)iBx)*dBinWidth; Int_t iBx1=0; dDTentry < 0 ? iBx1=iBx-1 : iBx1=iBx+1; Double_t w0=1.-TMath::Abs(dDTentry)/dBinWidth; Double_t w1=1.-w0; if(iBx1<0) iBx1=0; if(iBx1>nbClDelTofBinX-1) iBx1=nbClDelTofBinX-1; dDelTof=fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx][iSel]*w0 + fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx1][iSel]*w1; //dDelTof *= dDist; // has to be consistent with fhTRpcCluDelTof filling LOG(debug1)< DelT %6.1f", iSmType, iSm, iRpc, iSel, dTRef-dTTrig[iSel], iBx, iBx1, dDTentry, dDelTof) ; } dTTcor[iSel]=dDelTof*dZsign[iSel]; dTcor[iSel]=pHit->GetTime()-dDelTof-dTTrig[iSel]; // Double_t dAvTot=0.5*(pDig0->GetTot()+pDig1->GetTot()); (VF) not used } // if(iLink==0) if ( dTcor[iSel]==0. ) continue; LOG(debug)<GetNofLinks(),iSel,iSmType,iSm,iRpc, iCh0, iCh1, iS0, iS1, dTTrig[iSel], dDelTof, fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->GetEntries(), fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->GetEntries()) ; if( fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->GetEntries() != fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->GetEntries() ) LOG(error) << Form(" Inconsistent walk histograms -> debugging necessary ... for %d, %d, %d, %d, %d, %d, %d ", fiNevtBuild, iDetIndx, iSel, iCh0, iCh1, iS0, iS1) ; LOG(debug1)<GetTot(),pDig0->GetTime() +((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTcor[iSel]-dTTrig[iSel], iS1,pDig1->GetTot(),pDig1->GetTime() +((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTcor[iSel]-dTTrig[iSel]) ; fhTRpcCluWalk[iDetIndx][iSel][iCh0][iS0]->Fill(pDig0->GetTot(), //(pDig0->GetTime()+((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]); // dTcor[iSel]+(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)); dTcor[iSel]); fhTRpcCluWalk[iDetIndx][iSel][iCh1][iS1]->Fill(pDig1->GetTot(), //(pDig1->GetTime()+((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]); //dTcor[iSel]+(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)); dTcor[iSel]); fhTRpcCluWalk2[iDetIndx][iSel]->Fill(pDig0->GetTot(),pDig1->GetTot(),dTcor[iSel]); fhTRpcCluAvWalk[iDetIndx][iSel]->Fill(pDig0->GetTot(), //(pDig0->GetTime()+((1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]); //dTcor[iSel]+(1.-2.*pDig0->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)); dTcor[iSel]); fhTRpcCluAvWalk[iDetIndx][iSel]->Fill(pDig1->GetTot(), //(pDig1->GetTime()+((1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc))-dTTrig[iSel])-dTTcor[iSel]); //dTcor[iSel]+(1.-2.*pDig1->GetSide())*hitpos_local[1]/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)); dTcor[iSel]); if(iLink==0){ // Fill histo only once (for 1. digi entry) //fhTRpcCluDelTof[iDetIndx][iSel]->Fill(dTRef-dTTrig[iSel],dTcor[iSel]/dDist); fhTRpcCluDelTof[iDetIndx][iSel]->Fill(dTRef-dTTrig[iSel],dTcor[iSel]); fhTSmCluTOff[iSmType][iSel]->Fill((Double_t)(iSm*iNbRpc+iRpc),dTcor[iSel]); fhTSmCluTRun[iSmType][iSel]->Fill(fdEvent,dTcor[iSel]); if( iDetId != (pTrig[iSel]->GetAddress() & DetMask) ){ // transform matched hit-pair back into detector frame hitpos[0]=pHit->GetX()-dzscal*pTrig[iSel]->GetX() + fChannelInfo->GetX(); hitpos[1]=pHit->GetY()-dzscal*pTrig[iSel]->GetY() + fChannelInfo->GetY(); hitpos[2]=pHit->GetZ(); gGeoManager->MasterToLocal(hitpos, hitpos_local); // transform into local frame fhRpcCluDelMatPos[iDetIndx]->Fill((Double_t)iCh,hitpos_local[1]); fhRpcCluDelMatTOff[iDetIndx]->Fill((Double_t)iCh,(pHit->GetTime()-dTTrig[iSel])-dTTcor[iSel]); } } // iLink==0 condition end } // position condition end } // Match condition end } // closing of selector loop } // digi index range check condition end else { LOG(error)<<"CbmTofEventClusterizer::FillHistos: invalid digi index "<GetEntries() // << " in event " << XXX ; } } // iLink digi loop end; if (1GetTime(),2); Double_t dVar=dMeanTimeSquared/(dNstrips-1); //if(dVar<0.) dVar=0.; Double_t dTrms=TMath::Sqrt(dVar); LOG(debug)<Fill((Double_t)iCh,dTrms); pHit->SetTimeError(dTrms); } LOG(debug1)<<" Fill Time of iDetIndx "<GetAddress(),hitpos_local[1]) <<" for |y| <" <GetYaxis()->GetXmax(); if(TMath::Abs(hitpos_local[1])< (fhRpcCluPosition[iDetIndx]->GetYaxis()->GetXmax())){ if(dTRef !=0. && fTRefHits==1){ fhRpcCluTOff[iDetIndx]->Fill((Double_t)iCh,pHit->GetTime()-dTRef); fhSmCluTOff[iSmType]->Fill((Double_t)(iSm*iNbRpc+iRpc),pHit->GetTime()-dTRef); for (Int_t iSel=0; iSelGetTime()-dTTrig[iSel],fvLastHits[iSmType][iSm][iRpc][iCh].size()); if( pHit->GetAddress() == pTrig[iSel]->GetAddress() ) continue; if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1) { // check for previous hits in memory time interval std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); itL--; for(size_t iH=0; iHGetTime()-(*itL)->GetTime()); } } // LOG(debug1); // fill Time Offset histograms without velocity spread (DelTof) correction if ( pBeamRef != NULL ) if(TMath::Abs(pBeamRef->GetTime()-pTrig[iSel]->GetTime()) < fdDelTofMax) { // if(TMath::Abs(pBeamRef->GetTime()-pTrig[iSel]->GetTime()) < fhTRpcCluTOff[iIndexDut][iSel]->GetYaxis()->GetXmax()) { fhTRpcCluTOff[iDetIndx][iSel]->Fill((Double_t)iCh,pHit->GetTime()-dTTrig[iSel]);// -dTTcor[iSel] only valid for matches if( digiMatch->GetNofLinks()>0 ) fhTRpcCluTofOff[iDetIndx][iSel]->Fill((Double_t)iCh, pHit->GetTime()-dTTrig[iSel]); // valid for beam experiments // pHit->GetTime()-pBeamRef->GetTime()); // shift cluster time to beamcounter time // pHit->GetTime()-pBeamRef->GetTime()-fdToDAv*pTrig[iSel]->GetR());// valid for beam experiments } if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>1){ // check for previous hits in memory time interval std::list::iterator itL=fvLastHits[iSmType][iSm][iRpc][iCh].end(); itL--; for(Int_t iH=0; iH<1; iH++){ // use only last hit // for(Int_t iH=0; iHGetTime()-(*itL)->GetTime(); if(dTsinceLast > fdMemoryTime) LOG(fatal)< %f", iSmType,iSm,iRpc,iCh,dTsinceLast,fdMemoryTime ) ; fhTRpcCluTOffDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(dTsinceLast), pHit->GetTime()-dTTrig[iSel]); fhTRpcCluMemMulDTLastHits[iDetIndx][iSel]->Fill(TMath::Log10(dTsinceLast), fvLastHits[iSmType][iSm][iRpc][iCh].size()-1); } } } } } } // iHitInd hit loop end for( Int_t iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++ ){ for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ ) { LOG(debug1)<<"CbmTofEventClusterizer::FillHistos: " <Fill( fviClusterSize[iSmType][iRpc][uCluster]); fhClusterSizeType->Fill( fviClusterSize[iSmType][iRpc][uCluster], 40*iSmType + iRpc ); //FIXME - hardwired constant if( kFALSE) // 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(kFALSE) // 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); } // if(0cd(); // fhClustBuildTime->Write(); for(Int_t iDetIndx=0; iDetIndx< fDigiBdfPar->GetNbDet(); iDetIndx++){ if(NULL == fhRpcCluMul[iDetIndx]) continue; // fhRpcCluMul[iDetIndx]->Write(); // fhRpcCluRate[iDetIndx]->Write(); // fhRpcCluPosition[iDetIndx]->Write(); // fhRpcCluDelPos[iDetIndx]->Write(); // fhRpcCluTOff[iDetIndx]->Write(); // fhRpcCluDelTOff[iDetIndx]->Write(); // fhRpcCluTrms[iDetIndx]->Write(); // fhRpcCluTot[iDetIndx]->Write(); // fhRpcCluAvWalk[iDetIndx]->Write(); // fhRpcCluAvLnWalk[iDetIndx]->Write(); // fhRpcDTLastHits[iDetIndx]->Write(); // fhRpcDTLastHits_Tot[iDetIndx]->Write(); // fhRpcDTLastHits_CluSize[iDetIndx]->Write(); LOG(debug)<<"Write triggered Histos for Det Ind "<GetDetUId(iDetIndx)); for (Int_t iSel=0; iSelWrite(); // fhTRpcCluPosition[iDetIndx][iSel]->Write(); // fhTRpcCluTOff[iDetIndx][iSel]->Write(); // fhTRpcCluTot[iDetIndx][iSel]->Write(); // fhTRpcCluAvWalk[iDetIndx][iSel]->Write(); } Int_t iUniqueId = fDigiBdfPar->GetDetUId( iDetIndx ); Int_t iSmAddr = iUniqueId & SelMask; Int_t iSmType = CbmTofAddress::GetSmType( iUniqueId ); Int_t iSm = CbmTofAddress::GetSmId( iUniqueId ); Int_t iRpc = CbmTofAddress::GetRpcId( iUniqueId ); Int_t iNent =0; if(fCalSel > -1){ if(NULL == fhTRpcCluAvWalk[iDetIndx][fCalSel]) continue; iNent = (Int_t) fhTRpcCluAvWalk[iDetIndx][fCalSel]->GetEntries(); } else{ if(NULL == fhRpcCluAvWalk[iDetIndx]) continue; iNent = (Int_t) fhRpcCluAvWalk[iDetIndx]->GetEntries(); } if(0==iNent){ LOG(info)<<"WriteHistos: No entries in Walk histos for " <<"Smtype"< Comment to remove warning because set but never used TH2 *htempTot = NULL; TProfile *htempTot_pfx = NULL; TH1 *htempTot_Mean = NULL; TH1 *htempTot_Off = NULL; if(-1ProfileX("_pfx",1,fhRpcCluPosition[iDetIndx]->GetNbinsY()); //htempPos = fhTRpcCluPosition[iDetIndx][fCalSel]; //htempPos_pfx = fhTRpcCluPosition[iDetIndx][fCalSel]->ProfileX("_pfx",1,fhTRpcCluPosition[iDetIndx][fCalSel]->GetNbinsY()); htempTOff = fhTRpcCluTOff[iDetIndx][fCalSel]; // -> Comment to remove warning because set but never used htempTOff_pfx = htempTOff->ProfileX("_pfx",1,fhTRpcCluTOff[iDetIndx][fCalSel]->GetNbinsY()); htempTOff_px = htempTOff->ProjectionX("_px",1,fhTRpcCluTOff[iDetIndx][fCalSel]->GetNbinsY()); for(Int_t iCh=0; iChGetNbinsX(); iCh++) { // use peak value to prevent out of update range TH1* htempTOff_py = htempTOff->ProjectionY("_py",iCh+1,iCh+1); Double_t Ymax=htempTOff_py->GetMaximum(); if (Ymax > 0.) { Int_t iBmax=htempTOff_py->GetMaximumBin(); Double_t dTOffmax=htempTOff_py->GetXaxis()->GetBinCenter(iBmax); if ( TMath::Abs(dTOffmax) > 0.3*htempTOff_py->GetXaxis()->GetXmax()) { LOG(debug) << "Use Maximum of TOff in ch " << iCh << " of histo " << htempTOff-> GetName() << ": " << dTOffmax << ", " << htempTOff_py->GetXaxis()->GetXmax() <<" instead of " << htempTOff_pfx->GetBinContent(iCh+1) ; htempTOff_pfx->SetBinContent(iCh+1,dTOffmax); htempTOff_pfx->SetBinEntries(iCh+1,1); } } } htempTot = fhTRpcCluTot[iDetIndx][fCalSel]; htempTot_pfx = fhTRpcCluTot[iDetIndx][fCalSel]->ProfileX("_pfx",1,fhTRpcCluTot[iDetIndx][fCalSel]->GetNbinsY()); hAvPos_pfx = fhTSmCluPosition[iSmType][fCalSel]->ProfileX("_pfx",1,fhTSmCluPosition[iSmType][fCalSel]->GetNbinsY()); hAvTOff_pfx = fhTSmCluTOff[iSmType][fCalSel]->ProfileX("_pfx",1,fhTSmCluTOff[iSmType][fCalSel]->GetNbinsY(),"s"); }else // all triggers { htempPos = fhRpcCluPosition[iDetIndx]; htempTot = fhRpcCluTot[iDetIndx]; htempTot_pfx = fhRpcCluTot[iDetIndx]->ProfileX("_pfx",1,fhRpcCluTot[iDetIndx]->GetNbinsY()); hAvPos_pfx = fhSmCluPosition[iSmType]->ProfileX("_pfx",1,fhSmCluPosition[iSmType]->GetNbinsY()); hAvTOff_pfx = fhSmCluTOff[iSmType]->ProfileX("_pfx",1,fhSmCluTOff[iSmType]->GetNbinsY()); switch(fCalSel) { case -1: // take corrections from untriggered distributions htempPos_pfx = fhRpcCluPosition[iDetIndx]->ProfileX("_pfx",1,fhRpcCluPosition[iDetIndx]->GetNbinsY()); // htempTOff = fhRpcCluTOff[iDetIndx]; // -> Comment to remove warning because set but never used htempTOff_pfx = fhRpcCluTOff[iDetIndx]->ProfileX("_pfx",1,fhRpcCluTOff[iDetIndx]->GetNbinsY(),"s"); htempTOff_px = fhRpcCluTOff[iDetIndx]->ProjectionX("_px",1,fhRpcCluTOff[iDetIndx]->GetNbinsY()); break; case -2: //take corrections from Cluster deviations htempPos_pfx = fhRpcCluDelPos[iDetIndx]->ProfileX("_pfx",1,fhRpcCluDelPos[iDetIndx]->GetNbinsY()); // htempTOff = fhRpcCluDelTOff[iDetIndx]; // -> Comment to remove warning because set but never used htempTOff_pfx = fhRpcCluDelTOff[iDetIndx]->ProfileX("_pfx",1,fhRpcCluDelTOff[iDetIndx]->GetNbinsY()); htempTOff_px = fhRpcCluDelTOff[iDetIndx]->ProjectionX("_px",1,fhRpcCluDelTOff[iDetIndx]->GetNbinsY()); break; case -3: // take corrections from deviations to matched trigger hit htempPos_pfx = fhRpcCluDelMatPos[iDetIndx]->ProfileX("_pfx",1,fhRpcCluDelMatPos[iDetIndx]->GetNbinsY()); // htempTOff = fhRpcCluDelMatTOff[iDetIndx]; // -> Comment to remove warning because set but never used htempTOff_pfx = fhRpcCluDelMatTOff[iDetIndx]->ProfileX("_pfx",1,fhRpcCluDelMatTOff[iDetIndx]->GetNbinsY()); htempTOff_px = fhRpcCluDelMatTOff[iDetIndx]->ProjectionX("_px",1,fhRpcCluDelMatTOff[iDetIndx]->GetNbinsY()); break; case -4: // shift all detectors without match requirement to beam counter times { Int_t iCalSel=0; htempPos_pfx = fhTRpcCluPosition[iDetIndx][iCalSel]->ProfileX("_pfx",1,fhTRpcCluPosition[iDetIndx][iCalSel]->GetNbinsY()); htempTOff_pfx = fhTRpcCluTOff[iDetIndx][iCalSel]->ProfileX("_pfx",1,fhTRpcCluTofOff[iDetIndx][iCalSel]->GetNbinsY(),"s"); htempTOff_px = fhTRpcCluTofOff[iDetIndx][iCalSel]->ProjectionX("_px",1,fhTRpcCluTofOff[iDetIndx][iCalSel]->GetNbinsY()); } break; case -5: // shift all detectors without match requirement to beam counter times { Int_t iCalSel=1; htempPos_pfx = fhTRpcCluPosition[iDetIndx][iCalSel]->ProfileX("_pfx",1,fhTRpcCluPosition[iDetIndx][iCalSel]->GetNbinsY()); htempTOff_pfx = fhTRpcCluTOff[iDetIndx][iCalSel]->ProfileX("_pfx",1,fhTRpcCluTofOff[iDetIndx][iCalSel]->GetNbinsY(),"s"); htempTOff_px = fhTRpcCluTofOff[iDetIndx][iCalSel]->ProjectionX("_px",1,fhTRpcCluTofOff[iDetIndx][iCalSel]->GetNbinsY()); } break; } } if(NULL == htempPos_pfx) { LOG(info)<<"WriteHistos: Projections not available, continue " ; continue; } htempTot_Mean=htempTot_pfx->ProjectionX("_Mean"); htempTot_Off =htempTot_pfx->ProjectionX("_Off"); htempPos_pfx->SetName( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Pos_pfx",iSmType,iSm,iRpc) ); htempTOff_pfx->SetName( Form("cl_CorSmT%01d_sm%03d_rpc%03d_TOff_pfx",iSmType,iSm,iRpc) ); htempTot_pfx->SetName( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_pfx",iSmType,iSm,iRpc) ); htempTot_Mean->SetName( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Mean",iSmType,iSm,iRpc) ); htempTot_Off->SetName( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Tot_Off",iSmType,iSm,iRpc) ); hAvPos_pfx->SetName( Form("cl_CorSmT%01d_Pos_pfx",iSmType) ); hAvTOff_pfx->SetName( Form("cl_CorSmT%01d_TOff_pfx",iSmType) ); switch(fCalMode%10){ case 0 : { // Initialize htempTot_Off->Reset(); // prepare TotOffset histo TH1 *hbins[200]; Int_t nbins=htempTot->GetNbinsX(); for (int i=0;iProjectionY(Form("bin%d",i+1),i+1,i+1); /* Double_t Ymax=hbins[i]->GetMaximum();*/ Int_t iBmax=hbins[i]->GetMaximumBin(); TAxis *xaxis = hbins[i]->GetXaxis(); Double_t Xmax=xaxis->GetBinCenter(iBmax); Double_t XOff=Xmax-fTotPreRange; XOff = (Double_t)(Int_t)XOff; if(XOff<0) XOff=0; htempTot_Off->SetBinContent(i+1,XOff); Double_t Xmean=htempTot_Mean->GetBinContent(i+1); if(XmeanSetBinContent(i+1,(Xmean-XOff)); if ( htempTot_Mean->GetBinContent(i+1) != (Xmean-XOff)) LOG(warning)<<"Tot numbers not stored properly for " <GetBinContent(i+1),Xmean-XOff) ; } htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); } break; case 1 : //save offsets, update walks { Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); LOG(debug)<<"WriteHistos: restore Offsets and Gains and save Walk for " <<"Smtype"<Reset(); //reset to restore means of original histos htempTOff_pfx->Reset(); htempTot_Mean->Reset(); htempTot_Off->Reset(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { Double_t YMean=fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*0.5 *(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]-fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); Double_t TMean=0.5*(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]+fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); htempPos_pfx->Fill(iCh,YMean); if( ((TProfile *)htempPos_pfx)->GetBinContent(iCh+1)!=YMean) { LOG(error)<<"WriteHistos: restore unsuccessful! ch "<GetBinContent(iCh)<<","<GetBinContent(iCh+1) <<","<GetBinContent(iCh+2) <<", expected "<Fill(iCh,TMean); for(Int_t iSide=0; iSide<2; iSide++){ htempTot_Mean->SetBinContent(iCh*2+1+iSide, fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] ); //nh +1 empirical(?) htempTot_Off->SetBinContent(iCh*2+1+iSide,fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]); } } LOG(debug1)<<" Offset, gain restoring done ... "; htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); for(Int_t iSel=0; iSelcd(); // TH1D *hCorDelTof =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); gDirectory->cd( curdir->GetPath() ); if (NULL!=hCorDelTof) { TH1D *hCorDelTofout=(TH1D*)hCorDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); hCorDelTofout->Write(); }else { LOG(debug)<<" No CorDelTof histo " <GetNbChan(iSmType,iRpc)<<" channels" ; for( Int_t iCh=0; iCh< fDigiBdfPar->GetNbChan( iSmType, iRpc ); iCh++){ TH2 *h2tmp0; TH2 *h2tmp1; if(!fEnableAvWalk){ if(-1GetEntries(); if(iCh==0) // condition to print message only once LOG(info)<Write(); // h2tmp1->Write(); if(-1ProfileX("_pfx",1,h2tmp0->GetNbinsY()); TProfile *htmp1 = h2tmp1->ProfileX("_pfx",1,h2tmp1->GetNbinsY()); TH1D *h1tmp0 = h2tmp0->ProjectionX("_px",1,h2tmp0->GetNbinsY()); TH1D *h1tmp1 = h2tmp1->ProjectionX("_px",1,h2tmp1->GetNbinsY()); TH1D *h1ytmp0 = h2tmp0->ProjectionY("_py",1,nbClWalkBinX); // preserve means TH1D *h1ytmp1 = h2tmp1->ProjectionY("_py",1,nbClWalkBinX); Double_t dWMean0 = h1ytmp0->GetMean(); Double_t dWMean1 = h1ytmp1->GetMean(); Double_t dWMean = 0.5*(dWMean0+dWMean1); Int_t iWalkUpd=2; // Walk update mode flag //if(5==iSmType || 8==iSmType || 2==iSmType) iWalkUpd=0; // keep both sides consistent for diamonds and pads if(5==iSmType || 8==iSmType) iWalkUpd=0; // keep both sides consistent for diamonds and pads (Cern2016) for(Int_t iWx=0; iWxGetBinContent(iWx+1)>WalkNHmin && h1tmp1->GetBinContent(iWx+1)>WalkNHmin ){ // preserve y - position (difference) on average Double_t dWcor=(((TProfile *)htmp0)->GetBinContent(iWx+1) + ((TProfile *)htmp1)->GetBinContent(iWx+1))*0.5; fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]+=dWcor-dWMean; fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]+=dWcor-dWMean; LOG(info) << Form("Walk for TSR %d%d%d%d Tot %d set to %f",iSmType,iSm,iRpc,iCh,iWx, fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) ; } break; case 1 : if(h1tmp0->GetBinContent(iWx+1)>WalkNHmin && h1tmp1->GetBinContent(iWx+1)>WalkNHmin) { Double_t dWcor0 = ((TProfile *)htmp0)->GetBinContent(iWx+1)-dWMean0; Double_t dWcor1 = ((TProfile *)htmp1)->GetBinContent(iWx+1)-dWMean1; Double_t dWcor = 0.5*(dWcor0 + dWcor1); fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]+=dWcor; //-dWMean0; fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]+=dWcor; //-dWMean1; if(iCh==0 && iSmType==9 && iSm==0 && h1tmp0->GetBinContent(iWx+1)>WalkNHmin) LOG(debug) <<"Update Walk Sm = "<GetBinContent(iWx+1)<<" cts: " <GetBinContent(iWx+1) << " - " << dWMean0 <<" -> "<GetBinContent(iWx+1) << " - " << dWMean1 <<" -> "<GetBinContent(iWx+1)>WalkNHmin && h1tmp1->GetBinContent(iWx+1)>WalkNHmin) { Double_t dWcor0 = ((TProfile *)htmp0)->GetBinContent(iWx+1)-dWMean0; Double_t dWcor1 = ((TProfile *)htmp1)->GetBinContent(iWx+1)-dWMean1; //Double_t dWcor = 0.5*(dWcor0 + dWcor1); fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]+=dWcor0; fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]+=dWcor1; } break; default: ; } } h1tmp0->Reset(); h1tmp1->Reset(); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); Int_t iTry=3; while(iTry-- > 0 && h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { h1tmp0->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); } if (iTry==0) LOG(error)<<"writing not successful for "<GetName()<<", attempts left: "<GetBinContent(iWx+1) <<", expected "< 0 && h1tmp1->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]) { h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); } if (iTry==0) LOG(error)<<"writing not successful for "<GetName()<<", attempts left: "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); h1tmp0->Smooth(iNWalkSmooth); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); h1tmp1->Smooth(iNWalkSmooth); h1tmp1->Write(); } } }else{ // preserve whatever is there for fCalSmAddr ! for( Int_t iCh = 0; iCh < fDigiBdfPar->GetNbChan( iSmType, iRpc ); iCh++ ) // restore old values { // TProfile *htmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used // TProfile *htmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used TH1D *h1tmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProjectionX("_px",1,nbClWalkBinY); TH1D *h1tmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProjectionX("_px",1,nbClWalkBinY); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); if( h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { LOG(error)<<"WriteHistos: restore unsuccessful! iWx "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp0->Write(); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp1->Write(); h1tmp1->Write(); } } } break; case 2 : //update time offsets from positions and times with Sm averages, save walks and DELTOF { Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); if((fCalSmAddr < 0) || (fCalSmAddr != iSmAddr)){ // select detectors for updating offsets LOG(debug)<<"WriteHistos: (case 2) update Offsets and keep Gains, Walk and DELTOF for " <<"Smtype"<GetBinContent(iSm*iNbRpc+iRpc+1); if(dVscal==0.) dVscal=1.; Double_t YMean=((TProfile *)hAvPos_pfx)->GetBinContent(iB+1); //nh +1 empirical(?) htempPos_py=htempPos->ProjectionY(Form("%s_py",htempPos->GetName()),1,iNbCh); if(htempPos_py->GetEntries() > fdYFitMin && fPosYMaxScal < 1.1) { LOG(debug1)<GetName(),(Int_t)htempPos_py->GetEntries()) ; CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, 0); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(warning)<GetSizey()); TF1 *ff=htempPos_py->GetFunction("YBox"); if(NULL != ff){ LOG(info) << "FRes YBox "<GetEntries()<<" entries in TSR "<GetChisquare()/ff->GetNDF() << Form(", striplen (%5.2f), %4.2f: %7.2f +/- %5.2f, pos res %5.2f +/- %5.2f at y_cen = %5.2f +/- %5.2f", fChannelInfo->GetSizey(),dVscal, 2.*ff->GetParameter(1),2.*ff->GetParError(1), ff->GetParameter(2),ff->GetParError(2), ff->GetParameter(3),ff->GetParError(3)) ; if( TMath::Abs(fChannelInfo->GetSizey()-2.*ff->GetParameter(1))/fChannelInfo->GetSizey()<0.2 && TMath::Abs(ff->GetParError(1)/ff->GetParameter(1))<0.2 ) // && ff->GetChisquare() < 500.) //FIXME - constants! { if(TMath::Abs(ff->GetParameter(3)-YMean)<0.5*fChannelInfo->GetSizey()){ YMean=ff->GetParameter(3); Double_t dV =dVscal*fChannelInfo->GetSizey()/(2.*ff->GetParameter(1)); fhSmCluSvel[iSmType]->Fill((Double_t)(iSm*iNbRpc+iRpc),dV); } } } } TH1* hAvTOff_py = fhSmCluTOff[iSmType]->ProjectionY("_py",iB+1,iB+1); Double_t Ymax=hAvTOff_py->GetMaximum(); Double_t dTOffmax=0.; if (Ymax > 0.) { Int_t iBmax=hAvTOff_py->GetMaximumBin(); dTOffmax=hAvTOff_py->GetXaxis()->GetBinCenter(iBmax); } Double_t TMean=((TProfile *)hAvTOff_pfx)->GetBinContent(iB+1); if(TMath::Abs(dTOffmax-TMean) > 2.*TMean) { TMean=dTOffmax; LOG(debug) << "Use peak position for TOff of TSR "<GetBinError(iB+1); Double_t dTYOff=YMean/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) ; if (fiBeamRefAddr == iUniqueId ) TMean=0.; // don't shift reference counter LOG(debug)< Correct TSR %d%d%d by TMean=%8.2f, TYOff=%8.2f, TWidth=%8.2f, ",iSmType,iSm,iRpc,TMean,dTYOff,TWidth); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) // update Offset and Gain { fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] += -dTYOff + TMean ; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] += +dTYOff + TMean ; LOG(debug)<<"FillCalHist:" <<" SmT "<< iSmType<<" Sm "< " << Form(" %f, %f, %f, %f ",fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1], fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1]) ; } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) } htempPos_pfx->Reset(); //reset to store new values htempTOff_pfx->Reset(); htempTot_Mean->Reset(); htempTot_Off->Reset(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) // store new values { Double_t YMean=fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*0.5 *(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]-fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); Double_t TMean=0.5*(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]+fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); htempPos_pfx->Fill(iCh,YMean); if( ((TProfile *)htempPos_pfx)->GetBinContent(iCh+1)!=YMean) { LOG(error)<<"WriteHistos: restore unsuccessful! ch "<GetBinContent(iCh)<<","<GetBinContent(iCh+1) <<","<GetBinContent(iCh+2) <<", expected "<Fill(iCh,TMean); for(Int_t iSide=0; iSide<2; iSide++){ htempTot_Mean->SetBinContent(iCh*2+1+iSide, fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] ); //nh +1 empirical(?) htempTot_Off->SetBinContent(iCh*2+1+iSide,fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]); } } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) LOG(debug1)<<" Updating done ... write to file "; htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); // store old DELTOF histos LOG(debug)<<" Copy old DelTof histos from "<< gDirectory->GetName()<<" to file "; for(Int_t iSel=0; iSelcd(); // TH1D *hCorDelTof =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); gDirectory->cd( curdir->GetPath() ); if (NULL!=hCorDelTof) { TH1D *hCorDelTofout=(TH1D*)hCorDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); hCorDelTofout->Write(); }else { LOG(debug)<<" No CorDelTof histo " <ProfileX("_pfx",1,nbClWalkBinY); (VF) not used // TProfile *htmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used TH1D *h1tmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProjectionX("_px",1,nbClWalkBinY); TH1D *h1tmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProjectionX("_px",1,nbClWalkBinY); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); if( h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { LOG(error)<<"WriteHistos: restore unsuccessful! iWx "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp0->Write(); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp1->Write(); h1tmp1->Write(); } } break; case 3 : //update offsets, gains, save walks and DELTOF { Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); if((fCalSmAddr < 0) || (fCalSmAddr != iSmAddr) ){ // select detectors for updating offsets LOG(info)<<"WriteHistos (calMode==3): update Offsets and Gains, keep Walk and DelTof for " <<"Smtype"<GetMean(2); } */ Double_t dVscal=1.; Double_t dVW=1.; if(0) // NULL != fhSmCluSvel[iSmType]) { dVscal=fhSmCluSvel[iSmType]->GetBinContent(iSm*iNbRpc+iRpc+1); if(dVscal==0.) dVscal=1.; dVW=fhSmCluSvel[iSmType]->GetBinEffectiveEntries(iSm*iNbRpc+iRpc+1); dVW *= 50. ; //(Double_t)iNbCh; if(dVW<0.1) dVW=0.1; } // determine average values htempPos_py=htempPos->ProjectionY(Form("%s_py",htempPos->GetName()),1,iNbCh); Double_t dYMeanAv=0.; Double_t dYMeanFit=0.; if(htempPos_py->GetEntries() > fdYFitMin && fPosYMaxScal < 1.1) { dYMeanAv=htempPos_py->GetMean(); LOG(debug1)<GetName(),(Int_t)htempPos_py->GetEntries()) ; CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, 0); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(warning)<GetSizey()); TF1 *ff=htempPos_py->GetFunction("YBox"); if(NULL != ff){ if( TMath::Abs(fChannelInfo->GetSizey()-2.*ff->GetParameter(1))/fChannelInfo->GetSizey()<0.2 && TMath::Abs(ff->GetParError(1)/ff->GetParameter(1))<0.2 ) { Double_t dV =dVscal*fChannelInfo->GetSizey()/(2.*ff->GetParameter(1)); LOG(info) << "FAvRes YBox "<GetEntries()<<" entries in TSR "<fCstatu <<", chi2 "<GetChisquare()/ff->GetNDF() << Form(", striplen (%5.2f): %7.2f+/-%5.2f, pos res %5.2f+/-%5.2f at y_cen = %5.2f+/-%5.2f", fChannelInfo->GetSizey(), 2.*ff->GetParameter(1),2.*ff->GetParError(1), ff->GetParameter(2),ff->GetParError(2), ff->GetParameter(3),ff->GetParError(3)) ; if(TMath::Abs(ff->GetParameter(3)-dYMeanAv)<0.5*fChannelInfo->GetSizey()){ dYMeanFit=ff->GetParameter(3); fhSmCluSvel[iSmType]->Fill((Double_t)(iSm*iNbRpc+iRpc),dV,dVW); for (Int_t iPar=0; iPar<4; iPar++) fhSmCluFpar[iSmType][iPar]->Fill((Double_t)(iSm*iNbRpc+iRpc),ff->GetParameter(2+iPar)); } }else { LOG(info) << "FAvBad YBox "<GetEntries()<<" entries in "<GetChisquare() << Form(", striplen (%5.2f), %4.2f: %7.2f +/- %5.2f, pos res %5.2f +/- %5.2f at y_cen = %5.2f +/- %5.2f", fChannelInfo->GetSizey(),dVscal, 2.*ff->GetParameter(1),2.*ff->GetParError(1), ff->GetParameter(2),ff->GetParError(2), ff->GetParameter(3),ff->GetParError(3)); } } else { LOG(info) << "FAvFailed for TSR "<GetBinContent(iCh+1); //set default YMean += dYShift; htempPos_py=htempPos->ProjectionY(Form("%s_py%02d",htempPos->GetName(),iCh),iCh+1,iCh+1); if(htempPos_py->GetEntries() > fdYFitMin && fPosYMaxScal < -1.1) { //disabled LOG(debug1)<GetName(),iCh,(Int_t)htempPos_py->GetEntries()) ; CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, iCh); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(warning)<GetBinContent(iSm*iNbRpc+iRpc+1); //LOG(info) << Form("Call yFit with %6.3f, %6.3f, %6.3f, %6.3f",fp[0],fp[1],fp[2],fp[3]) // ; Double_t *fpp = &fp[0]; fit_ybox(htempPos_py,0.5*fChannelInfo->GetSizey(),fpp); TF1 *ff=htempPos_py->GetFunction("YBox"); if(NULL != ff){ if( TMath::Abs(fChannelInfo->GetSizey()-2.*ff->GetParameter(1))/fChannelInfo->GetSizey()<0.1 && TMath::Abs(ff->GetParError(1)/ff->GetParameter(1))<0.05 ) //&& ff->GetChisquare() < 200.) //FIXME - constants! { if(TMath::Abs(ff->GetParameter(3)-YMean)<0.5*fChannelInfo->GetSizey()){ YMean=ff->GetParameter(3); Double_t dV =dVscal*fChannelInfo->GetSizey()/(2.*ff->GetParameter(1)); fhSmCluSvel[iSmType]->Fill((Double_t)(iSm*iNbRpc+iRpc),dV,dVW); LOG(info) << "FRes YBox "<GetEntries()<<" entries in "<GetChisquare() << Form(", striplen (%5.2f), %4.2f -> %4.2f, %4.1f: %7.2f+/-%5.2f, pos res %5.2f+/-%5.2f at y_cen = %5.2f+/-%5.2f", fChannelInfo->GetSizey(),dVscal,dV,dVW, 2.*ff->GetParameter(1),2.*ff->GetParError(1), ff->GetParameter(2),ff->GetParError(2), ff->GetParameter(3),ff->GetParError(3)); for (Int_t iPar=0; iPar<4; iPar++) fhSmCluFpar[iSmType][iPar]->Fill((Double_t)(iSm*iNbRpc+iRpc),ff->GetParameter(2+iPar)); } } else { //YMean=0.; // no new info available - did not help! LOG(info) << "FBad YBox "<GetEntries()<<" entries in "<GetChisquare() << Form(", striplen (%5.2f), %4.2f: %7.2f +/- %5.2f, pos res %5.2f +/- %5.2f at y_cen = %5.2f +/- %5.2f", fChannelInfo->GetSizey(),dVscal, 2.*ff->GetParameter(1),2.*ff->GetParError(1), ff->GetParameter(2),ff->GetParError(2), ff->GetParameter(3),ff->GetParError(3)); } } } Double_t TMean=((TProfile *)htempTOff_pfx)->GetBinContent(iCh+1); Double_t dTYOff=YMean/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) ; if (fiBeamRefAddr == iUniqueId ) { // don't shift reference counter // if (fiBeamRefType == iSmType && fiBeamRefSm == iSm && fiBeamRefDet == iRpc) { // don't shift reference counter on average if(iCh==0) { Double_t dW=0.; for( Int_t iRefCh = 0; iRefCh < iNbCh; iRefCh++ ) { if( 0 != ((TH1 *)htempTOff_px)->GetBinContent(iRefCh+1)) { dW+= ((TH1 *)htempTOff_px)->GetBinContent(iRefCh+1); TWMean+= ((TProfile *)htempTOff_pfx)->GetBinContent(iRefCh+1) * ((TH1 *)htempTOff_px)->GetBinContent(iRefCh+1); } } if(dW >0.) TWMean /= dW; else TWMean=0.; } if( htempTOff_px->GetBinContent(iCh+1) > 0. ) LOG(info) <GetBinContent(iCh+1), TMean, ((TProfile *)hAvTOff_pfx)->GetBinContent(iSm*iNbRpc+iRpc+1), TWMean, fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]+TMean-TWMean); // TMean-=((TProfile *)hAvTOff_pfx)->GetBinContent(iSm*iNbRpc+iRpc+1); TMean -= TWMean; } if(htempTOff_px->GetBinContent(iCh+1)>WalkNHmin){ fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] += -dTYOff + TMean; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] += +dTYOff + TMean; LOG(info)< new Off %8.3f,%8.3f ", fCalMode,fCalSel,fTRefMode,iSmType,iSm,iRpc,iCh,htempTOff_px->GetBinContent(iCh+1), dTYOff,TMean, fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]) ; } /* Double_t TotMean=((TProfile *)htempTot_pfx)->GetBinContent(iCh+1); //nh +1 empirical(!) if(0.001 < TotMean){ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0] *= fdTTotMean / TotMean; fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1] *= fdTTotMean / TotMean; } */ for(Int_t iSide=0; iSide<2; iSide++){ Int_t ib=iCh*2+1+iSide; TH1 * hbin=htempTot->ProjectionY(Form("bin%d",ib),ib,ib); if(100>hbin->GetEntries()) continue; //request min number of entries /* Double_t Ymax=hbin->GetMaximum();*/ Int_t iBmax=hbin->GetMaximumBin(); TAxis *xaxis = hbin->GetXaxis(); Double_t Xmax=xaxis->GetBinCenter(iBmax)/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]; Double_t XOff=Xmax-fTotPreRange; if(0){//TMath::Abs(XOff - fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide])>100){ LOG(warning)<<"XOff changed for " <GetBinContent(ib); Double_t TotMean=hbin->GetMean(); if(15==iSmType){ LOG(warning)<<"Gain for " <GetBinContent(ib), fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide],fdTTotMean / TotMean); } if(0.001 < TotMean){ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] *= fdTTotMean / TotMean; } } if(5==iSmType && fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] != fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]) {// diamond LOG(warning)<<"CbmTofEventClusterizer::FillCalHist:" <<" SmT "<< iSmType<<" Sm "< " <Reset(); //reset to store new values htempTOff_pfx->Reset(); htempTot_Mean->Reset(); htempTot_Off->Reset(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) // store new values { Double_t YMean=fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*0.5 *(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]-fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); Double_t TMean=0.5*(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]+fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); htempPos_pfx->Fill(iCh,YMean); if( ((TProfile *)htempPos_pfx)->GetBinContent(iCh+1)!=YMean) { LOG(error)<<"WriteHistos: restore unsuccessful! ch "<GetBinContent(iCh)<<","<GetBinContent(iCh+1) <<","<GetBinContent(iCh+2) <<", expected "<Fill(iCh,TMean); for(Int_t iSide=0; iSide<2; iSide++){ htempTot_Mean->SetBinContent(iCh*2+1+iSide, fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] ); htempTot_Off->SetBinContent(iCh*2+1+iSide,fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]); } // htempTot_pfx->Fill(iCh,fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1]); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) LOG(debug1)<<" Updating done ... write to file "; htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); // store old DELTOF histos LOG(debug)<<" Copy old DelTof histos from "<< gDirectory->GetName()<<" to file "; for(Int_t iSel=0; iSelcd(); // TH1D *hCorDelTof =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); gDirectory->cd( curdir->GetPath() ); if (NULL!=hCorDelTof) { TH1D *hCorDelTofout=(TH1D*)hCorDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); hCorDelTofout->Write(); }else { LOG(debug)<<" No CorDelTof histo " <ProfileX("_pfx",1,nbClWalkBinY); (VF) not used // TProfile *htmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used TH1D *h1tmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProjectionX("_px",1,nbClWalkBinY); TH1D *h1tmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProjectionX("_px",1,nbClWalkBinY); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); if( h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { LOG(error)<<"WriteHistos: restore unsuccessful! iWx "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp0->Write(); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp1->Write(); h1tmp1->Write(); } } break; case 4 : //update DelTof, save offsets, gains and walks { Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); LOG(debug)<<"WriteHistos: restore Offsets, Gains and Walk, save DelTof for " <<"Smtype"<Reset(); //reset to restore mean of original histos htempTOff_pfx->Reset(); htempTot_Mean->Reset(); htempTot_Off->Reset(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) { Double_t YMean=fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*0.5 *(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]-fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); Double_t TMean=0.5*(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]+fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); htempPos_pfx->Fill(iCh,YMean); if( ((TProfile *)htempPos_pfx)->GetBinContent(iCh+1)!=YMean) { LOG(error)<<"WriteHistos: restore unsuccessful! ch "<GetBinContent(iCh)<<","<GetBinContent(iCh+1) <<","<GetBinContent(iCh+2) <<", expected "<Fill(iCh,TMean); for(Int_t iSide=0; iSide<2; iSide++){ htempTot_Mean->SetBinContent(iCh*2+1+iSide, fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] ); //nh +1 empirical(?) htempTot_Off->SetBinContent(iCh*2+1+iSide,fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]); } } LOG(debug1)<<" Restoring of Offsets and Gains done ... "; htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); // restore walk histos for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) // restore old values { if(NULL == fhRpcCluWalk[iDetIndx][iCh][0]) break; // TProfile *htmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used // TProfile *htmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used TH1D *h1tmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProjectionX("_px",1,nbClWalkBinY); TH1D *h1tmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProjectionX("_px",1,nbClWalkBinY); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); if( h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { LOG(error)<<"WriteHistos: restore unsuccessful! iWx "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp0->Write(); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp1->Write(); h1tmp1->Write(); } // generate/update DelTof corrections if((fCalSmAddr < 0 && -fCalSmAddr != iSmAddr) || (fCalSmAddr == iSmAddr)) // select detectors for determination of DelTof correction { if( fiBeamRefAddr == iSmAddr ) continue; // no DelTof correction for Diamonds for(Int_t iSel=0; iSelGetEntries(); // h2tmp->Write(); TProfile *htmp = h2tmp->ProfileX("_pfx",1,h2tmp->GetNbinsY()); TH1D *h1tmp = h2tmp->ProjectionX("_px",1,h2tmp->GetNbinsY()); /* TH1D *h1ytmp = h2tmp->ProjectionY("_py",1,h2tmp->GetNbinsX());*/ Double_t dDelMean=0.;//h1ytmp->GetMean();// inspect normalization, interferes with mode 3 for diamonds, nh, 10.01.2015 (?) Double_t dNEntriesSum=0.; for(Int_t iBx=0; iBxGetBinContent(iBx+1); if(dNEntries > WalkNHmin) // modify, request sufficient # of entries fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx][iSel]+=((TProfile *)htmp)->GetBinContent(iBx+1); dDelMean += fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx][iSel]*dNEntries; dNEntriesSum += dNEntries; } dDelMean /= dNEntriesSum; LOG(debug)<SetBinContent(iBx+1,fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx][iSel]-dDelMean); //h1tmp->SetBinContent(iBx+1,fvCPDelTof[iSmType][iSm*iNbRpc+iRpc][iBx][iSel]); } h1tmp->SetName(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof", iSmType, iSm, iRpc,iSel)); h1tmp->Write(); } }else{ // copy existing histograms for(Int_t iSel=0; iSelcd(); // TH1D *hCorDelTof =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); gDirectory->cd( curdir->GetPath() ); if (NULL!=hCorDelTof) { TH1D *hCorDelTofout=(TH1D*)hCorDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); LOG(debug)<<" Save existing CorDelTof histo " <Write(); }else { LOG(debug)<<" No CorDelTof histo " <GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); if((fCalSmAddr < 0) || (fCalSmAddr != iSmAddr) ){ // select detectors for updating offsets LOG(info)<<"WriteHistos (calMode==5): update Offsets and Gains, keep Walk and DelTof for " <<"Smtype"<GetBinContent(iCh+1); //nh +1 empirical(?) Double_t TMean=0.; Double_t dTYOff=YMean/fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) ; if(htempTOff_px->GetBinContent(iCh+1)>WalkNHmin){ fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] += -dTYOff + TMean; fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1] += +dTYOff + TMean; } LOG(debug3)<GetBinContent(iCh+1), fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0], fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]) ; /* Double_t TotMean=((TProfile *)htempTot_pfx)->GetBinContent(iCh+1); //nh +1 empirical(!) if(0.001 < TotMean){ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][0] *= fdTTotMean / TotMean; fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1] *= fdTTotMean / TotMean; } */ for(Int_t iSide=0; iSide<2; iSide++){ Int_t ib=iCh*2+1+iSide; TH1 * hbin=htempTot->ProjectionY(Form("bin%d",ib),ib,ib); if(100>hbin->GetEntries()) continue; //request min number of entries /* Double_t Ymax=hbin->GetMaximum();*/ Int_t iBmax=hbin->GetMaximumBin(); TAxis *xaxis = hbin->GetXaxis(); Double_t Xmax=xaxis->GetBinCenter(iBmax)/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]; Double_t XOff=Xmax-fTotPreRange; if(0){//TMath::Abs(XOff - fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide])>100){ LOG(warning)<<"XOff changed for " <GetBinContent(ib); Double_t TotMean=hbin->GetMean(); if(15==iSmType){ LOG(warning)<<"Gain for " <GetBinContent(ib), fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide],fdTTotMean / TotMean) ; } if(0.001 < TotMean){ fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] *= fdTTotMean / TotMean; } } if(5==iSmType && fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0] != fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]) {// diamond LOG(warning)<<"CbmTofEventClusterizer::FillCalHist:" <<" SmT "<< iSmType<<" Sm "< " <Reset(); //reset to store new values htempTOff_pfx->Reset(); htempTot_Mean->Reset(); htempTot_Off->Reset(); for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) // store new values { Double_t YMean=fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)*0.5 *(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]-fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); Double_t TMean=0.5*(fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][1]+fvCPTOff[iSmType][iSm*iNbRpc+iRpc][iCh][0]); htempPos_pfx->Fill(iCh,YMean); if( ((TProfile *)htempPos_pfx)->GetBinContent(iCh+1)!=YMean) { LOG(error)<<"WriteHistos: restore unsuccessful! ch "<GetBinContent(iCh)<<","<GetBinContent(iCh+1) <<","<GetBinContent(iCh+2) <<", expected "<Fill(iCh,TMean); for(Int_t iSide=0; iSide<2; iSide++){ htempTot_Mean->SetBinContent(iCh*2+1+iSide, fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][iSide] ); htempTot_Off->SetBinContent(iCh*2+1+iSide,fvCPTotOff[iSmType][iSm*iNbRpc+iRpc][iCh][iSide]); } // htempTot_pfx->Fill(iCh,fdTTotMean/fvCPTotGain[iSmType][iSm*iNbRpc+iRpc][iCh][1]); } // for( Int_t iCh = 0; iCh < iNbCh; iCh++ ) LOG(debug1)<<" Updating done ... write to file "; htempPos_pfx->Write(); htempTOff_pfx->Write(); // htempTot_pfx->Write(); htempTot_Mean->Write(); htempTot_Off->Write(); // store old DELTOF histos LOG(debug)<<" Copy old DelTof histos from "<< gDirectory->GetName()<<" to file "; for(Int_t iSel=0; iSelcd(); // TH1D *hCorDelTof =(TH1D*) gDirectory->FindObjectAny( Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); gDirectory->cd( curdir->GetPath() ); if (NULL!=hCorDelTof) { TH1D *hCorDelTofout=(TH1D*)hCorDelTof->Clone(Form("cl_CorSmT%01d_sm%03d_rpc%03d_Sel%02d_DelTof",iSmType,iSm,iRpc,iSel)); hCorDelTofout->Write(); }else { LOG(debug)<<" No CorDelTof histo " <ProfileX("_pfx",1,nbClWalkBinY); (VF) not used // TProfile *htmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProfileX("_pfx",1,nbClWalkBinY); (VF) not used TH1D *h1tmp0 = fhRpcCluWalk[iDetIndx][iCh][0]->ProjectionX("_px",1,nbClWalkBinY); TH1D *h1tmp1 = fhRpcCluWalk[iDetIndx][iCh][1]->ProjectionX("_px",1,nbClWalkBinY); for(Int_t iWx=0; iWxSetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]); h1tmp1->SetBinContent(iWx+1,fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][1][iWx]); if( h1tmp0->GetBinContent(iWx+1)!=fvCPWalk[iSmType][iSm*iNbRpc+iRpc][iCh][0][iWx]) { LOG(error)<<"WriteHistos: restore unsuccessful! iWx "<GetBinContent(iWx+1) <<", expected "<SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S0_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp0->Write(); h1tmp0->Write(); h1tmp1->SetName(Form("Cor_SmT%01d_sm%03d_rpc%03d_Ch%03d_S1_Walk_px", iSmType, iSm, iRpc, iCh )); // htmp1->Write(); h1tmp1->Write(); } } break; default: LOG(debug)<<"WriteHistos: update mode " <Write(); // fhDigSpacDifClust->Write(); // fhDigTimeDifClust->Write(); // fhDigDistClust->Write(); // fhClustSizeDifX->Write(); // fhClustSizeDifY->Write(); // fhChDifDifX->Write(); // fhChDifDifY->Write(); // fhNbSameSide->Write(); // fhNbDigiPerChan->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() ) for (Int_t iS=0; iS< fDigiBdfPar->GetNbSmTypes(); iS++){ if (NULL == fhSmCluPosition[iS]) continue; // fhSmCluPosition[iS]->Write(); // fhSmCluTOff[iS]->Write(); fhSmCluSvel[iS]->Write(); for (Int_t iPar=0; iPar<4; iPar++) fhSmCluFpar[iS][iPar]->Write(); for (Int_t iSel=0; iSelWrite(); // fhTSmCluTOff[iS][iSel]->Write(); // fhTSmCluTRun[iS][iSel]->Write(); } } gDirectory->cd( oldir->GetPath() ); fHist->Close(); return kTRUE; } Bool_t CbmTofEventClusterizer::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 CbmTofEventClusterizer::BuildClusters() { Int_t iMess =0; //gGeoManager->SetTopVolume( gGeoManager->FindVolumeFast("tof_v14a") ); gGeoManager->CdTop(); if(NULL == fTofDigisColl) { LOG(info) <<" No RawDigis defined ! Check! " ; return kFALSE; } fiNevtBuild++; LOG(debug)<<"Build clusters from " <GetEntries()<<" digis in event "<GetEntries(); if (iNbTofDigi > 100000) { LOG(warning) << "Too many digis in event " << fiNevtBuild ; return kFALSE; } if( kFALSE ) { // Duplicate type "5" - digis Int_t iNbDigi=iNbTofDigi; for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { CbmTofDigiExp *pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd ); if( pDigi->GetType() == 5 ) { CbmTofDigiExp *pDigiN = new((*fTofDigisColl)[iNbDigi++]) CbmTofDigiExp( *pDigi ); pDigiN->SetAddress(pDigi->GetSm(),pDigi->GetRpc(),pDigi->GetChannel(), (0 == pDigi->GetSide()) ? 1 : 0,pDigi->GetType()); } } iNbTofDigi = fTofDigisColl->GetEntries(); // Update } if( kTRUE ) { for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { CbmTofDigiExp *pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd ); LOG(debug)<GetAddress()) <<" SmT " << pDigi->GetType() <<" Sm " << pDigi->GetSm() <<" Rpc "<< pDigi->GetRpc() <<" Ch " << pDigi->GetChannel() <<" S " << pDigi->GetSide() <<" : " << pDigi->ToString() // <<" Time "<GetTime() // <<" Tot " <GetTot() ; Int_t iDetIndx= fDigiBdfPar->GetDetInd( pDigi->GetAddress() ); if (fDigiBdfPar->GetNbDet()-1< %d,0 ",iDetIndx,fDigiBdfPar->GetNbDet()) ; break; } if (NULL == fhRpcDigiCor[iDetIndx] ) { if ( 100GetAddress()) ; continue; } size_t iDigiCh=pDigi->GetChannel()*2+pDigi->GetSide(); if (iDigiCh < fvTimeLastDigi[iDetIndx].size()) { if(fvTimeLastDigi[iDetIndx][iDigiCh] > 0) { if( fdStartAna10s > 0.) { Double_t dTimeAna10s = (pDigi->GetTime() - fdStartAna10s)/1.E9; if ( dTimeAna10s < fdSpillDuration) fhRpcDigiDTLD[iDetIndx]->Fill( iDigiCh, (pDigi->GetTime()-fvTimeLastDigi[iDetIndx][iDigiCh])/1.E9 ); } } fvTimeLastDigi[iDetIndx][iDigiCh] = pDigi->GetTime(); } Double_t dTDifMin=dDoubleMax; CbmTofDigiExp *pDigi2Min=NULL; // for (Int_t iDigI2 =iDigInd+1; iDigI2At( iDigI2 ); if( iDetIndx == fDigiBdfPar->GetDetInd( pDigi2->GetAddress() )){ if(0.==pDigi->GetSide() && 1.==pDigi2->GetSide()){ fhRpcDigiCor[iDetIndx]->Fill( pDigi->GetChannel(), pDigi2->GetChannel() ); } else { if (1.==pDigi->GetSide() && 0.==pDigi2->GetSide()){ fhRpcDigiCor[iDetIndx]->Fill( pDigi2->GetChannel(), pDigi->GetChannel() ); } } if( pDigi->GetSide() != pDigi2->GetSide() ) { if (pDigi->GetChannel() == pDigi2->GetChannel()) { Double_t dTDif=TMath::Abs(pDigi->GetTime()-pDigi2->GetTime()); if(dTDifGetChannel() - pDigi2->GetChannel()) == 1 ) { // opposite side missing, neighbouring channel has hit on opposite side // FIXME // check that same side digi of neighbouring channel is absent Int_t iDigI3 =0; for (; iDigI3At( iDigI3 ); if(pDigi3->GetSide() == pDigi->GetSide() && pDigi2->GetChannel() == pDigi3->GetChannel()) break; } if(iDigI3 == iNbTofDigi) // same side neighbour did not fire { Int_t iCorMode=0; // Missing hit correction mode switch(iCorMode){ case 0: // no action break; case 1: // shift found hit LOG(debug2) << Form("shift channel %d%d%d%d%d and ",(Int_t)pDigi->GetType(),(Int_t)pDigi->GetSm(),(Int_t)pDigi->GetRpc(),(Int_t)pDigi->GetChannel(),(Int_t)pDigi->GetSide()) << Form(" %d%d%d%d%d ",(Int_t)pDigi2->GetType(),(Int_t)pDigi2->GetSm(),(Int_t)pDigi2->GetRpc(),(Int_t)pDigi2->GetChannel(),(Int_t)pDigi2->GetSide()) ; //if(pDigi->GetTime() < pDigi2->GetTime()) if(pDigi->GetSide() == 0) pDigi2->SetAddress(pDigi->GetSm(),pDigi->GetRpc(),pDigi->GetChannel(),1-pDigi->GetSide(),pDigi->GetType()); else pDigi->SetAddress(pDigi2->GetSm(),pDigi2->GetRpc(),pDigi2->GetChannel(),1-pDigi2->GetSide(),pDigi2->GetType()); LOG(debug2) << Form("resultchannel %d%d%d%d%d and ",(Int_t)pDigi->GetType(),(Int_t)pDigi->GetSm(),(Int_t)pDigi->GetRpc(),(Int_t)pDigi->GetChannel(),(Int_t)pDigi->GetSide()) << Form(" %d%d%d%d%d ",(Int_t)pDigi2->GetType(),(Int_t)pDigi2->GetSm(),(Int_t)pDigi2->GetRpc(),(Int_t)pDigi2->GetChannel(),(Int_t)pDigi2->GetSide()) ; break; case 2: // insert missing hits CbmTofDigiExp *pDigiN = new((*fTofDigisColl)[iNbTofDigi++]) CbmTofDigiExp( *pDigi ); pDigiN->SetAddress(pDigi->GetSm(),pDigi->GetRpc(),pDigi2->GetChannel(),pDigi->GetSide(),pDigi->GetType()); pDigiN->SetTot(pDigi2->GetTot()); CbmTofDigiExp *pDigiN2 = new((*fTofDigisColl)[iNbTofDigi++]) CbmTofDigiExp( *pDigi2 ); pDigiN2->SetAddress(pDigi2->GetSm(),pDigi2->GetRpc(),pDigi->GetChannel(),pDigi2->GetSide(),pDigi2->GetType()); pDigiN2->SetTot(pDigi->GetTot()); break; } } } } } } if( pDigi2Min !=NULL ){ CbmTofDetectorInfo xDetInfo(kTof, pDigi->GetType(), pDigi->GetSm(), pDigi->GetRpc(), 0, pDigi->GetChannel()); Int_t iChId = fTofId->SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(warning)<GetSigVel(pDigi->GetType(),pDigi->GetSm(),pDigi->GetRpc()) * dTDifMin * 0.5 < fPosYMaxScal*fChannelInfo->GetSizey()) { //check consistency if(8==pDigi->GetType() || 5==pDigi->GetType()){ if(pDigi->GetTime() != pDigi2Min->GetTime()){ if(fiMsgCnt-- >0){ LOG(warning)<<" BuildClusters: Inconsistent duplicated digis in event " << fiNevtBuild <<", Ind: "<ToString() ; LOG(warning)<<" "<ToString() ; } pDigi2Min->SetTot(pDigi->GetTot()); pDigi2Min->SetTime(pDigi->GetTime()); } } // average ToTs! temporary fix, FIXME /* Double_t dAvTot=0.5*(pDigi->GetTot()+pDigi2Min->GetTot()); pDigi->SetTot(dAvTot); pDigi2Min->SetTot(dAvTot); LOG(debug)<<" BuildClusters: TDif "<ToString() ; LOG(debug)<<" "<ToString() ; */ } } } } // Calibrate RawDigis if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) { CbmTofDigiExp *pDigi; // CbmTofDigiExp *pCalDigi=NULL; (VF) not used CalibRawDigis(); // Then loop over the digis array and store the Digis in separate vectors for // each RPC modules iNbTofDigi = fTofCalDigisColl->GetEntries(); for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { pDigi = (CbmTofDigiExp*) fTofCalDigisColl->At( iDigInd ); LOG(debug1)<<"AC " // After Calibration <GetAddress())<<" TSRC " <GetType() <GetSm() <GetRpc() <GetChannel())<<" " <GetSide()<<" " <GetTime())<<" " <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(), pDigi->GetRpc()) > 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); } else { LOG(info)<<"Skip2 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++ ) // inspect digi array Int_t iNbDet=fDigiBdfPar->GetNbDet(); for(Int_t iDetIndx=0; iDetIndxGetNbChan(iSmType,iRpc); for(Int_t iStrip=0; iStripGetNbRpc(iSmType) + iRpc] [iStrip].size() ; //LOG(info)<<"Inspect TSRC "< 0 ) { fhRpcDigiMul[iDetIndx] ->Fill(iStrip,iDigiMul); if(iDigiMul==1) { fhRpcDigiStatus[iDetIndx] ->Fill(iStrip,0); if(iStrip > 0) if ( fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc(iSmType ) + iRpc] [iStrip-1].size() > 1) { fhRpcDigiStatus[iDetIndx] ->Fill(iStrip,1); if( TMath::Abs( fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc( iSmType ) + iRpc] [iStrip][0]->GetTime() -fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc( iSmType ) + iRpc] [iStrip-1][0]->GetTime()) < fMaxTimeDist ) fhRpcDigiStatus[iDetIndx] ->Fill(iStrip,3); } if(iStrip < iNbStrips-2) { if ( fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc( iSmType ) + iRpc] [iStrip+1].size() > 1) { fhRpcDigiStatus[iDetIndx] ->Fill(iStrip,2); if( TMath::Abs( fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc( iSmType ) + iRpc] [iStrip][0]->GetTime() -fStorDigiExp[iSmType] [iSm*fDigiBdfPar->GetNbRpc( iSmType ) + iRpc] [iStrip+1][0]->GetTime()) < fMaxTimeDist ) fhRpcDigiStatus[iDetIndx] ->Fill(iStrip,4); } } } } } } BuildHits(); } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) else { LOG(error)<<" Compressed Digis not implemented ... "; return kFALSE; // not implemented properly yet } return kTRUE; } Bool_t CbmTofEventClusterizer::MergeClusters() { // Merge clusters from neigbouring Rpc within a (Super)Module if(NULL == fTofHitsColl) { LOG(info) <<" No Hits defined ! Check! " ; return kFALSE; } // inspect hits for( Int_t iHitInd = 0; iHitInd < fTofHitsColl->GetEntries(); iHitInd++) { CbmTofHit *pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & SelMask); Int_t iSmType = CbmTofAddress::GetSmType( iDetId ); Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); if(iSmType != 5 && iSmType != 8 ) continue; // only merge diamonds and Pad LOG(debug)<<"MergeClusters: in SmT "<< iSmType << " for "<< iNbRpc << " Rpcs" ; if(iNbRpc>1) { // check for possible mergers Int_t iSm = CbmTofAddress::GetSmId( iDetId ); Int_t iRpc = CbmTofAddress::GetRpcId( iDetId ); Int_t iChId = pHit->GetAddress(); fChannelInfo = fDigiPar->GetCell( iChId ); Int_t iCh = CbmTofAddress::GetChannelId(iChId); LOG(debug)<<"MergeClusters: Check for mergers in " <GetEntries(); iHitInd2++) { CbmTofHit *pHit2 = (CbmTofHit*) fTofHitsColl->At( iHitInd2 ); if(NULL == pHit2) continue; Int_t iDetId2 = (pHit2->GetAddress() & SelMask); Int_t iSmType2 = CbmTofAddress::GetSmType( iDetId2 ); if(iSmType2 == iSmType) { Int_t iSm2 = CbmTofAddress::GetSmId( iDetId2 ); if(iSm2 == iSm || iSmType==5 ){ Int_t iRpc2 = CbmTofAddress::GetRpcId( iDetId2 ); if(TMath::Abs(iRpc-iRpc2)==1 || iSm2!=iSm){ // Found neighbour Int_t iChId2 = pHit2->GetAddress(); // CbmTofCell *fChannelInfo2 = fDigiPar->GetCell( iChId2 ); (VF) not used Int_t iCh2 = CbmTofAddress::GetChannelId(iChId2); Double_t xPos=pHit->GetX(); Double_t yPos=pHit->GetY(); Double_t tof=pHit->GetTime(); Double_t xPos2=pHit2->GetX(); Double_t yPos2=pHit2->GetY(); Double_t tof2=pHit2->GetTime(); LOG(debug)<<"MergeClusters: Found hit in neighbour " <At(iHitInd); Double_t dTot=0; for (Int_t iLink=0; iLinkGetNofLinks(); iLink+=2){ // loop over digis CbmLink L0 = digiMatch->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); Int_t iDigInd1=(digiMatch->GetLink(iLink+1)).GetIndex(); if (iDigInd0 < fTofCalDigisColl->GetEntries() && iDigInd1 < fTofCalDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); CbmTofDigiExp *pDig1 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd1)); dTot += pDig0->GetTot(); dTot += pDig1->GetTot(); } } CbmMatch* digiMatch2=(CbmMatch *)fTofDigiMatchColl->At(iHitInd2); Double_t dTot2=0; for (Int_t iLink=0; iLinkGetNofLinks(); iLink+=2){ // loop over digis CbmLink L0 = digiMatch2->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); Int_t iDigInd1=(digiMatch2->GetLink(iLink+1)).GetIndex(); if (iDigInd0 < fTofCalDigisColl->GetEntries() && iDigInd1 < fTofCalDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd0)); CbmTofDigiExp *pDig1 = (CbmTofDigiExp*) (fTofCalDigisColl->At(iDigInd1)); dTot2 += pDig0->GetTot(); dTot2 += pDig1->GetTot(); digiMatch->AddLink(CbmLink(pDig0->GetTot(),iDigInd0,fiOutputTreeEntry,fiFileIndex)); digiMatch->AddLink(CbmLink(pDig1->GetTot(),iDigInd1,fiOutputTreeEntry,fiFileIndex)); } } LOG(debug)<<"MergeClusters: Found merger in neighbour " <GetEntries()) <SetX(dxPosM); pHit->SetY(dyPosM); pHit->SetTime(dtofM); // remove merged hit at iHitInd2 and update digiMatch fTofHitsColl->RemoveAt( iHitInd2 ); fTofDigiMatchColl->RemoveAt( iHitInd2 ); fTofDigiMatchColl->Compress(); fTofHitsColl->Compress(); LOG(debug)<<"MergeClusters: Compress TClonesArrays to " <GetEntries()<<", " <GetEntries() ; /* for(Int_t i=iHitInd2; iGetEntries(); i++){ // update RefLinks CbmTofHit *pHiti = (CbmTofHit*) fTofHitsColl->At( i ); pHiti->SetRefId(i); } */ //check merged hit (cluster) //pHit->Print(); } } } } } } } return kTRUE; } static Double_t f1_xboxe(double *x, double *par){ double xx = x[0]; double wx = 1. - par[4]*TMath::Power(xx+par[5],2); double xboxe = par[0]*0.25*(1.+TMath::Erf(( xx+par[1]-par[3])/par[2])) *(1.+TMath::Erf((-xx+par[1]+par[3])/par[2])); return xboxe*wx ; } void CbmTofEventClusterizer::fit_ybox(const char *hname) { TH1 *h1; h1=(TH1 *)gROOT->FindObjectAny(hname); if(NULL != h1){ fit_ybox(h1,0.); } } void CbmTofEventClusterizer::fit_ybox(TH1 *h1, Double_t ysize) { Double_t* fpar=NULL; fit_ybox(h1, ysize, fpar); } void CbmTofEventClusterizer::fit_ybox(TH1 *h1, Double_t ysize, Double_t* fpar=NULL) { TAxis *xaxis = h1->GetXaxis(); Double_t Ymin= xaxis->GetXmin(); Double_t Ymax= xaxis->GetXmax(); TF1 *f1=new TF1("YBox",f1_xboxe,Ymin,Ymax,6); Double_t yini=(h1->GetMaximum()+h1->GetMinimum())*0.5; if (ysize==0.) ysize=Ymax*0.8; f1->SetParameters(yini,ysize*0.5,1.,0.,0.,0.); // f1->SetParLimits(1,ysize*0.8,ysize*1.2); f1->SetParLimits(2,0.2,3.); f1->SetParLimits(3,-4.,4.); if(fpar != NULL) { Double_t fp[4]; for(Int_t i=0; i<4; i++) fp[i]=*fpar++; for(Int_t i=0; i<4; i++) f1->SetParameter(2+i,fp[i]); LOG(debug) << "Ini Fpar for " << h1->GetName() << " with " << Form(" %6.3f %6.3f %6.3f %6.3f ",fp[0],fp[1],fp[2],fp[3]) ; } h1->Fit("YBox","Q"); double res[10]; double err[10]; res[9]=f1->GetChisquare(); for (int i=0; i<6; i++) { res[i]=f1->GetParameter(i); err[i]=f1->GetParError(i); //cout << " FPar "<< i << ": " << res[i] << ", " << err[i] << endl; } LOG(debug) << "YBox Fit of "<GetName()<<" ended with chi2 = "<(fDigiBdfPar->GetNbSmTypes())) LOG(fatal)<GetNbSmTypes()) ; for (Int_t iSmType=0; iSmTypeGetNbSmTypes(); iSmType++ ){ if(fvLastHits[iSmType].size() != static_cast(fDigiBdfPar->GetNbSm( iSmType ))) LOG(fatal)<GetNbSm( iSmType ), iSmType) ; for( Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm( iSmType); iSm++ ){ if(fvLastHits[iSmType][iSm].size() != static_cast(fDigiBdfPar->GetNbRpc( iSmType ))) LOG(fatal)<GetNbRpc( iSmType ),iSmType,iSm) ; for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ ){ if(fvLastHits[iSmType][iSm][iRpc].size() != static_cast(fDigiBdfPar->GetNbChan( iSmType, iRpc ))) LOG(fatal)<GetNbChan( iSmType, iRpc ),iSmType,iSm,iRpc) ; for (Int_t iCh=0; iCh< fDigiBdfPar->GetNbChan( iSmType, iRpc ); iCh++ ) if(fvLastHits[iSmType][iSm][iRpc][iCh].size()>0){ CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, iCh); Int_t iAddr = fTofId->SetDetectorInfo( xDetInfo ); if( fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr ) LOG(fatal)<GetAddress(), fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()) ; } } } } LOG(debug) << Form("LH check passed for event %8.0f",fdEvent) ; } void CbmTofEventClusterizer::CleanLHMemory() { if(fvLastHits.size() != static_cast(fDigiBdfPar->GetNbSmTypes())) LOG(fatal)<GetNbSmTypes()) ; for (Int_t iSmType=0; iSmTypeGetNbSmTypes(); iSmType++ ){ if(fvLastHits[iSmType].size() != static_cast(fDigiBdfPar->GetNbSm( iSmType ))) LOG(fatal)<GetNbSm( iSmType ), iSmType) ; for( Int_t iSm = 0; iSm < fDigiBdfPar->GetNbSm( iSmType); iSm++ ){ if(fvLastHits[iSmType][iSm].size() != static_cast(fDigiBdfPar->GetNbRpc( iSmType ))) LOG(fatal)<GetNbRpc( iSmType ),iSmType,iSm) ; for( Int_t iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc( iSmType); iRpc++ ){ if(fvLastHits[iSmType][iSm][iRpc].size() != static_cast(fDigiBdfPar->GetNbChan( iSmType, iRpc ))) LOG(fatal)<GetNbChan( iSmType, iRpc ),iSmType,iSm,iRpc) ; for (Int_t iCh=0; iCh< fDigiBdfPar->GetNbChan( iSmType, iRpc ); iCh++ ) while(fvLastHits[iSmType][iSm][iRpc][iCh].size()>0){ CbmTofDetectorInfo xDetInfo(kTof, iSmType, iSm, iRpc, 0, iCh); Int_t iAddr = fTofId->SetDetectorInfo( xDetInfo ); if( fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetAddress() != iAddr ) LOG(fatal)<GetAddress(), fvLastHits[iSmType][iSm][iRpc][iCh].front()->GetTime()) ; fvLastHits[iSmType][iSm][iRpc][iCh].front()->Delete(); fvLastHits[iSmType][iSm][iRpc][iCh].pop_front(); } } } } LOG(info) << Form("LH cleaning done after %8.0f events",fdEvent) ; } Bool_t CbmTofEventClusterizer::AddNextChan(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iLastChan, Double_t dLastPosX, Double_t dLastPosY, Double_t dLastTime, Double_t dLastTotS){ // Int_t iNbSm = fDigiBdfPar->GetNbSm( iSmType); (VF) not used Int_t iNbRpc = fDigiBdfPar->GetNbRpc( iSmType); Int_t iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); // Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); (VF) not used Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm,iRpc,0,0,iSmType); Int_t iDetIndx = fDetIdIndexMap[iDetId]; // Detector Index Int_t iCh=iLastChan+1; while( fvDeadStrips[iDetIndx] & ( 1 << iCh )) { LOG(debug) << "Skip channel "<= iNbCh) return kFALSE; } LOG(debug1) << Form("Inspect channel TSRC %d%d%d%d at time %f, pos %f, size ",iSmType,iSm,iRpc,iCh,dLastTime,dLastPosY) << fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ; if (iCh == iNbCh) return kFALSE; if (0 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) return kFALSE; if( 0 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) fhNbDigiPerChan->Fill( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ); if( 1 < fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ) { Bool_t AddedHit=kFALSE; for(size_t i1=0; i1GetSide() == (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][i2])->GetSide() ) { i2++; continue; } // endif same side // 2 Digis, both sides present CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][i1]; CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][i2]; Double_t dTime = 0.5 * ( xDigiA->GetTime() + xDigiB->GetTime() ) ; if(TMath::Abs(dTime-dLastTime)SetDetectorInfo( xDetInfo ); fChannelInfo = fDigiPar->GetCell( iChId ); gGeoManager->FindNode(fChannelInfo->GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ()); Double_t dTimeDif = xDigiA->GetTime() - xDigiB->GetTime(); Double_t dPosY=0.; if( 1 == xDigiA->GetSide() ) dPosY = fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDif * 0.5; else dPosY = -fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDif * 0.5; if(TMath::Abs(dPosY - dLastPosY) < fdMaxSpaceDist ) { // append digi pair to current cluster Double_t dNClHits=(Double_t)(vDigiIndRef.size()/2); Double_t dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex(); Double_t dTotS = xDigiA->GetTot() + xDigiB->GetTot(); Double_t dNewTotS = (dLastTotS + dTotS); dLastPosX=(dLastPosX*dLastTotS + dPosX*dTotS)/dNewTotS; dLastPosY=(dLastPosY*dLastTotS + dPosY*dTotS)/dNewTotS; dLastTime=(dLastTime*dLastTotS + dTime*dTotS)/dNewTotS; dLastTotS = dNewTotS; // attach selected digis from pool Int_t Ind1 = fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][i1]; Int_t Ind2 = fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][i2]; vDigiIndRef.push_back( Ind1 ); vDigiIndRef.push_back( Ind2 ); // remove selected digis from pool fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+i1); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+i1); std::vector::iterator it; it = find (fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin(), fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].end(), Ind2); if (it != fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].end()){ auto ipos = it - fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin(); LOG(debug1)<<"Found i2 "<Rndm()-0.5)*fChannelInfo->GetSizex()*0.5; hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey()*0.5; } */ Double_t hitpos[3]; /*TGeoNode* cNode = */gGeoManager->GetCurrentNode(); /*TGeoHMatrix* cMatrix = */gGeoManager->GetCurrentMatrix(); gGeoManager->LocalToMaster(hitpos_local, hitpos); TVector3 hitPos(hitpos[0],hitpos[1],hitpos[2]); TVector3 hitPosErr(0.5,0.5,0.5); // FIXME including positioning uncertainty Int_t iChm=floor(dLastPosX/fChannelInfo->GetSizex())+iNbCh/2; if(iChm<0) iChm=0; if(iChm >iNbCh-1) iChm=iNbCh-1; iDetId = CbmTofAddress::GetUniqueAddress(iSm,iRpc,iChm,0,iSmType); Int_t iNbChanInHit=vDigiIndRef.size()/2; fviClusterMul[iSmType][iSm][iRpc]++; LOG(debug1)<<"Save A-Hit " << Form("%2d %2d 0x%08x %3d t %f, y %f ", fiNbHits,iNbChanInHit,iDetId,iLastChan,dLastTime,dLastPosY) <<", DigiSize: "< Tot // output hit new((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit); if(fdMemoryTime > 0.) { // memorize hit LH_store(iSmType,iSm,iRpc,iChm,pHit); }else{ pHit->Delete(); } CbmMatch* digiMatch = new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(); for (size_t i=0; iAt(vDigiIndRef.at(i))))->GetTot(); digiMatch->AddLink(CbmLink(dTot,vDigiIndRef.at(i),fiOutputTreeEntry,fiFileIndex)); } fiNbHits++; vDigiIndRef.clear(); return kTRUE; } void CbmTofEventClusterizer::LH_store(Int_t iSmType, Int_t iSm, Int_t iRpc, Int_t iChm, CbmTofHit *pHit){ if( fvLastHits[iSmType][iSm][iRpc][iChm].size() == 0 ) fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit); else { Double_t dLastTime=pHit->GetTime(); if(dLastTime >= fvLastHits[iSmType][iSm][iRpc][iChm].back()->GetTime()){ fvLastHits[iSmType][iSm][iRpc][iChm].push_back(pHit); LOG(debug)<GetAddress(), dLastTime, dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime() ) ; }else{ if(dLastTime >= fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime()){// hit has to be inserted in the proper place std::list::iterator it; for (it=fvLastHits[iSmType][iSm][iRpc][iChm].begin(); it != fvLastHits[iSmType][iSm][iRpc][iChm].end(); ++it) if( (*it)->GetTime()>dLastTime ) break; fvLastHits[iSmType][iSm][iRpc][iChm].insert(--it,pHit); Double_t deltaTime=dLastTime - (*it)->GetTime(); LOG(debug)<GetAddress(), deltaTime) ; }else{ // this hit is first Double_t deltaTime=dLastTime - fvLastHits[iSmType][iSm][iRpc][iChm].front()->GetTime(); LOG(debug)<GetAddress(), deltaTime) ; if(deltaTime==0.){ // remove hit, otherwise double entry? pHit->Delete(); }else{ fvLastHits[iSmType][iSm][iRpc][iChm].push_front(pHit); } } } } } Bool_t CbmTofEventClusterizer::BuildHits() { // 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; //vPtsRef.clear(); vDigiIndRef.clear(); // CbmTofCell *fTrafoCell=NULL; (VF) not used // Int_t iTrafoCell=-1; (VF) not used 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 dTimeDif = 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 iNbCh = fDigiBdfPar->GetNbChan( iSmType, iRpc ); Int_t iChType = fDigiBdfPar->GetChanType( iSmType, iRpc ); LOG(debug2)<<"RPC - Loop " << Form(" %3d %3d %3d %3d ",iSmType,iSm,iRpc,iChType) ; fviClusterMul[iSmType][iSm][iRpc]=0; Int_t iChId = 0; Int_t iDetId = CbmTofAddress::GetUniqueAddress(iSm,iRpc,0,0,iSmType);; Int_t iDetIndx = fDetIdIndexMap[iDetId]; // Detector Index if( 0 == iChType ){ // 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; // -> Comment to remove warning because set but never used dLastPosY = 0.0; dLastTime = 0.0; LOG(debug2)<<"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)<<"VDigisize " << Form(" T %3d Sm %3d R %3d Ch %3d Size %3lu ", iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; if (0 == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc].size()) continue; if( fvDeadStrips[iDetIndx] & ( 1 << iCh )) continue; // skip over dead channels 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() ) { while( (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetSide() == (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetSide() ) { // Not one Digi of each end! fiNbSameSide++; if(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()>2) { LOG(debug) << "SameSide Digis! on TSRC " << iSmType<GetTime()) << ", "<GetTime()) <<", DeltaT " <<(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1])->GetTime() - (fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0])->GetTime() <<", array size: " << fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ; if ( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][2]->GetSide() == fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]->GetSide() ) { LOG(debug) << "3 consecutive SameSide Digis! on TSRC " << iSmType<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() <<", array size: " << fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ; 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()); }else { if( fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][2]->GetTime() -fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]->GetTime() > fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][2]->GetTime() -fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1]->GetTime()) { 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()); } else { LOG(debug) << Form("Ev %8.0f, digis not properly time ordered, TSRCS %d%d%d%d%d ", fdEvent,iSmType,iSm,iRpc,iCh,(Int_t)fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]->GetSide()) ; fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+1); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].erase(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+1); } } }else{ LOG(debug2)<<"SameSide Erase fStor entries(d) "< fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) break; continue; } // same condition side end LOG(debug2) << "digis processing for " << Form(" SmT %3d Sm %3d Rpc %3d Ch %3d # %3lu ",iSmType,iSm,iRpc,iCh, fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) ; if(2 > fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) { LOG(debug)<SetDetectorInfo( xDetInfo ); Int_t iUCellId=CbmTofAddress::GetUniqueAddress(iSm,iRpc,iCh,0,iSmType); LOG(debug1)<< Form(" TSRC %d%d%d%d size %3lu ", iSmType,iSm,iRpc,iCh,fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()) << Form(" ChId: 0x%08x 0x%08x ",iChId,iUCellId) ; fChannelInfo = fDigiPar->GetCell( iChId ); if(NULL == fChannelInfo){ LOG(error)<<"CbmTofEventClusterizer::BuildClusters: no geometry info! " << Form(" %3d %3d %3d %3d 0x%08x 0x%08x ",iSmType, iSm, iRpc, iCh, iChId,iUCellId) ; break; } TGeoNode *fNode= // prepare local->global trafo gGeoManager->FindNode(fChannelInfo->GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ()); LOG(debug2)<GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ(),fNode) ; // fNode->Print(); CbmTofDigiExp * xDigiA = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][0]; CbmTofDigiExp * xDigiB = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][1]; LOG(debug2) << " "<ToString(); LOG(debug2) << " "<ToString(); dTimeDif = ( xDigiA->GetTime() - xDigiB->GetTime() ) ; if(5==iSmType && dTimeDif !=0.) { // FIXME -> Overflow treatment in calib/tdc/TMbsCalibTdcTof.cxx LOG(debug)<<"CbmTofEventClusterizer::BuildClusters: Diamond hit in " << iSm <<" with inconsistent digits " << xDigiA->GetTime() << ", " << xDigiB->GetTime() << " -> "<ToString(); LOG(debug) << " "<ToString(); } if( 1 == xDigiA->GetSide() ) // 0 is the top side, 1 is the bottom side dPosY = fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDif * 0.5; else // 0 is the bottom side, 1 is the top side dPosY = -fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDif * 0.5; if(TMath::Abs(dPosY) > fChannelInfo->GetSizey() && fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size()>2) { LOG(debug)<<"Hit candidate outside correlation window, check for better possible digis, " <<" mul "<< fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].size() ; CbmTofDigiExp * xDigiC = fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh][2]; Double_t dPosYN=0.; Double_t dTimeDifN=0; if( xDigiC->GetSide()==xDigiA->GetSide() ) dTimeDifN= xDigiC->GetTime() - xDigiB->GetTime(); else dTimeDifN= xDigiA->GetTime() - xDigiC->GetTime(); if( 1 == xDigiA->GetSide() ) dPosYN = fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDifN * 0.5; else dPosYN = -fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc) * dTimeDifN * 0.5; if(TMath::Abs(dPosYN)GetSide() <<", yPosNext "<GetSide()==xDigiA->GetSide() ) { xDigiA=xDigiC; 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()); }else{ xDigiB=xDigiC; fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].erase(++(fStorDigiExp[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+1)); fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].erase(++(fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh].begin()+1)); } } } if(xDigiA->GetSide() == xDigiB->GetSide()){ LOG(fatal)<<"Wrong combinations of digis " << fStorDigiInd[iSmType][iSm*iNbRpc+iRpc][iCh][0]<<"," << fStorDigiInd[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(); // use local coordinates, (0,0,0) is in the center of counter ? dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex(); dPosZ=0.; LOG(debug1) <<"NbChanInHit " << Form(" %3d %3d %3d %3d %3d 0x%p %1.0f Time %f PosX %f PosY %f Svel %f ", iNbChanInHit,iSmType,iRpc,iCh,iLastChan,xDigiA,xDigiA->GetSide(), dTime,dPosX,dPosY,fDigiBdfPar->GetSigVel(iSmType,iSm,iRpc)) // << Form( ", Offs %f, %f ",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) < fdMaxTimeDist && iLastChan == iCh - 1 && TMath::Abs( dPosY - dLastPosY) < fdMaxSpaceDist ) { // Add to cluster/hit dWeightedTime += dTime*dTotS; dWeightedPosX += dPosX*dTotS; dWeightedPosY += dPosY*dTotS; dWeightedPosZ += dPosZ*dTotS; dWeightsSum += dTotS; iNbChanInHit += 1; 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<<", " <Rndm()-0.5)*fChannelInfo->GetSizex(); hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey(); } */ Double_t hitpos[3]; TGeoNode* cNode = gGeoManager->GetCurrentNode(); /*TGeoHMatrix* cMatrix =*/ gGeoManager->GetCurrentMatrix(); //cNode->Print(); //cMatrix->Print(); gGeoManager->LocalToMaster(hitpos_local, hitpos); LOG(debug1)<< Form(" LocalToMaster for node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)", cNode, 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 TVector3 hitPosErr(0.5,0.5,0.5); // including positioning uncertainty /* TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation 0.5, // Use generic value 1.); */ //fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSigVel(iSmType,iRpc), // 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 // calc mean ch from dPosX=((Double_t)(-iNbCh/2 + iCh)+0.5)*fChannelInfo->GetSizex(); Int_t iChm=floor(dWeightedPosX/fChannelInfo->GetSizex())+iNbCh/2; if(iChm<0) iChm=0; if(iChm >iNbCh-1) iChm=iNbCh-1; 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] ); //} TString sRef=""; for (UInt_t i=0; i0){ CbmTofHit *pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits-1); if(iDetId == pHitL->GetAddress() && dWeightedTime==pHitL->GetTime()){ LOG(debug)<<"Store Hit twice? " <<" fiNbHits "<At(vDigiIndRef.at(i)); LOG(debug)<<" Digi "<ToString(); } CbmMatch* digiMatchL=(CbmMatch *)fTofDigiMatchColl->At(fiNbHits-1); for (Int_t i=0; iGetNofLinks();i++){ CbmLink L0 = digiMatchL->GetLink(i); Int_t iDigIndL=L0.GetIndex(); CbmTofDigiExp *pDigiC = (CbmTofDigiExp*) fTofCalDigisColl->At(iDigIndL); LOG(debug)<<" DigiL "<ToString(); } } } CbmTofHit *pHit = new CbmTofHit( iDetId, hitPos, hitPosErr, //local detector coordinates fiNbHits, // this number is used as reference!! dWeightedTime, vDigiIndRef.size(), // number of linked digis = 2*CluSize //vPtsRef.size(), // flag = number of TofPoints generating the cluster Int_t(dWeightsSum*10.)) ; //channel -> Tot //0) ; //channel // output hit new((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit); // memorize hit if(fdMemoryTime > 0.) { LH_store(iSmType,iSm,iRpc,iChm,pHit); }else{ pHit->Delete(); } /* new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(); CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits); */ CbmMatch* digiMatch = new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(); for (size_t i=0; iAt(vDigiIndRef.at(i))))->GetTot(); digiMatch->AddLink(CbmLink(dTot,vDigiIndRef.at(i),fiOutputTreeEntry,fiFileIndex)); } 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; dWeightsSum = dTotS; iNbChanInHit = 1; // Save pointer on CbmTofPoint // vPtsRef.push_back( (CbmTofPoint*)(xDigiA->GetLinks()) ); // Save next digi address LOG(debug2)<<" 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(debug)<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; dLastPosY = dPosY; dLastTime = dTime; if ( AddNextChan(iSmType,iSm,iRpc,iLastChan,dLastPosX,dLastPosY,dLastTime,dWeightsSum) ){ iNbChanInHit=0; // cluster already stored } } // while( 1 < 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)<<"finished V-RPC" << Form(" %3d %3d %3d %d %f %fx",iSmType,iSm,iRpc,fTofHitsColl->GetEntries(),dLastPosX,dLastPosY) ; } // else of if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) } // if( 0 == iChType) else { LOG(error)<<"=> Cluster building " <<"from digis to hits not implemented for pads, Sm type " <GetChanOrient( iSmType, iRpc ) ){ LOG(debug1)<<"H-Hit " ; } // if( 1 == fDigiBdfPar->GetChanOrient( iSmType, iRpc ) ) else { LOG(debug2)<<"V-Hit " ; // Save Hit dWeightedTime /= dWeightsSum; dWeightedPosX /= dWeightsSum; dWeightedPosY /= dWeightsSum; dWeightedPosZ /= dWeightsSum; //TVector3 hitPos(dWeightedPosX, dWeightedPosY, dWeightedPosZ); Double_t hitpos_local[3]={3*0.}; hitpos_local[0] = dWeightedPosX; hitpos_local[1] = dWeightedPosY; hitpos_local[2] = dWeightedPosZ; /* if( 5 == iSmType || 8 == iSmType) { // for PAD counters hitpos_local[0] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizex(); hitpos_local[1] = (gRandom->Rndm()-0.5)*fChannelInfo->GetSizey(); } */ Double_t hitpos[3]; TGeoNode* cNode= gGeoManager->GetCurrentNode(); /*TGeoHMatrix* cMatrix =*/ gGeoManager->GetCurrentMatrix(); //cNode->Print(); //cMatrix->Print(); gGeoManager->LocalToMaster(hitpos_local, hitpos); LOG(debug1)<< Form(" LocalToMaster for V-node %p: (%6.1f,%6.1f,%6.1f) ->(%6.1f,%6.1f,%6.1f)", cNode, hitpos_local[0], hitpos_local[1], hitpos_local[2], hitpos[0], hitpos[1], hitpos[2]) ; TVector3 hitPos(hitpos[0],hitpos[1],hitpos[2]); // Event 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(0.5,0.5,0.5); // including positioning uncertainty /* TVector3 hitPosErr( fChannelInfo->GetSizex()/TMath::Sqrt(12.0), // Single strips approximation 0.5, // Use generic value 1.); */ // fDigiBdfPar->GetFeeTimeRes() * fDigiBdfPar->GetSigVel(iSmType,iRpc), // 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 = iChId; Int_t iChm=floor(dWeightedPosX/fChannelInfo->GetSizex())+iNbCh/2; if(iChm<0) iChm=0; if(iChm >iNbCh-1) iChm=iNbCh-1; 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] ); TString cstr = "Save V-Hit "; cstr += Form(" %3d %3d 0x%08x %3d 0x%08x", // %3d %3d fiNbHits,iNbChanInHit,iDetId,iLastChan,iRefId); //vPtsRef.size(),vPtsRef[0]) // dWeightedTime,dWeightedPosY) cstr += Form(", DigiSize: %lu ",vDigiIndRef.size()); cstr += ", DigiInds: "; fviClusterMul[iSmType][iSm][iRpc]++; for (UInt_t i=0; i0){ CbmTofHit *pHitL = (CbmTofHit*) fTofHitsColl->At(fiNbHits-1); if(iDetId == pHitL->GetAddress() && dWeightedTime==pHitL->GetTime()) LOG(debug)<<"Store Hit twice? " <<" fiNbHits "< Tot // 0) ; //channel // vDigiIndRef); // output hit new((*fTofHitsColl)[fiNbHits]) CbmTofHit(*pHit); // memorize hit if(fdMemoryTime > 0.) { LH_store(iSmType,iSm,iRpc,iChm,pHit); }else{ pHit->Delete(); } /* new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(); CbmMatch* digiMatch = (CbmMatch *)fTofDigiMatchColl->At(fiNbHits); */ CbmMatch* digiMatch = new((*fTofDigiMatchColl)[fiNbHits]) CbmMatch(); for (size_t i=0; iAt(vDigiIndRef.at(i))))->GetTot(); digiMatch->AddLink(CbmLink(dTot,vDigiIndRef.at(i),fiOutputTreeEntry,fiFileIndex)); } 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(debug2)<<" Fini-A "<mChannelDeadTime; Int_t iNbTofDigi = fTofDigisColl->GetEntriesFast(); for( Int_t iDigInd = 0; iDigInd < iNbTofDigi; iDigInd++ ) { pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd ); if(fbSwapChannelSides && 5 != pDigi->GetType() && 8 != pDigi->GetType()) { pDigi->SetAddress(pDigi->GetSm(), pDigi->GetRpc(), pDigi->GetChannel(), (0 == pDigi->GetSide()) ? 1 : 0, pDigi->GetType()); } Int_t iAddr=pDigi->GetAddress(); LOG(debug1)<<"BC " // Before Calibration <GetAddress())<<" TSRC " <GetType() <GetSm() <GetRpc() <GetChannel())<<" " <GetSide()<<" " <GetTime())<<" " <GetTot() ; if(pDigi->GetType()==5 || pDigi->GetType() == 8) // for Pad counters generate fake digi to mockup a strip if(pDigi->GetSide()==1) continue; // skip one side to avoid double entries Bool_t bValid=kTRUE; std::map::iterator it; it = mChannelDeadTime.find(iAddr); if( it != mChannelDeadTime.end() ) { LOG(debug1)<<"CCT found valid ChannelDeadtime entry "<GetTime()-mChannelDeadTime[iAddr] ; if( (bValid = (pDigi->GetTime() > mChannelDeadTime[iAddr] + fdChannelDeadtime)) ) pCalDigi = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigiExp( *pDigi ); }else { pCalDigi = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigiExp( *pDigi ); } mChannelDeadTime[iAddr]=pDigi->GetTime(); if ( ! bValid ) continue; LOG(debug1)<<"DC " // After deadtime check. before Calibration <GetAddress())<<" TSRC " <GetType() <GetSm() <GetRpc() <GetChannel())<<" " <GetSide()<<" " <GetTime())<<" " <GetTot() ; if(fbPs2Ns) { pCalDigi->SetTime(pCalDigi->GetTime()/1000.); // for backward compatibility pCalDigi->SetTot(pCalDigi->GetTot()/1000.); // for backward compatibility } 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(), pDigi->GetRpc()) > pDigi->GetChannel() ) { LOG(debug2) <<" CluCal-Init: "<ToString(); // apply calibration vectors pCalDigi->SetTime(pCalDigi->GetTime()- // calibrate Digi Time fvCPTOff[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()]); LOG(debug2)<<" CluCal-TOff: "<ToString(); Double_t dTot = pCalDigi->GetTot()- // subtract Offset fvCPTotOff[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()]; if (dTot < 0.001) dTot=0.001; pCalDigi->SetTot(dTot * // calibrate Digi ToT fvCPTotGain[pDigi->GetType()] [pDigi->GetSm()*fDigiBdfPar->GetNbRpc( pDigi->GetType()) + pDigi->GetRpc()] [pDigi->GetChannel()] [pDigi->GetSide()]); // walk correction Double_t dTotBinSize = (fdTOTMax-fdTOTMin)/ nbClWalkBinX; Int_t iWx = (Int_t)((pCalDigi->GetTot()-fdTOTMin)/dTotBinSize); if (0>iWx) iWx=0; if (iWx>=nbClWalkBinX) iWx=nbClWalkBinX-1; Double_t dDTot = (pCalDigi->GetTot()-fdTOTMin)/dTotBinSize-(Double_t)iWx-0.5; Double_t dWT = fvCPWalk[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx]; if(dDTot > 0) { // linear interpolation to next bin if(iWx < nbClWalkBinX -1) { // linear interpolation to next bin dWT += dDTot * (fvCPWalk[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx+1] -fvCPWalk[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx]); //memory leak??? } }else // dDTot < 0, linear interpolation to next bin { if(0 < iWx) { // linear interpolation to next bin dWT -= dDTot * (fvCPWalk[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx-1] -fvCPWalk[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx]); //memory leak??? } } pCalDigi->SetTime(pCalDigi->GetTime() - dWT); // calibrate Digi Time LOG(debug2)<<" CluCal-Walk: "<ToString(); if(1) {//pDigi->GetType()==7 && pDigi->GetSm()==0){ LOG(debug)<<"BuildClusters: CalDigi " <GetType() <GetSm() <GetRpc() <GetChannel())) <GetSide() <GetTime()) <<", Tot "<GetTot() <<", TotGain "<< fvCPTotGain[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] <<", TotOff "<< fvCPTotOff[pCalDigi->GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] <<", Walk "<GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx] ; LOG(debug1)<<" dDTot "<GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx-1] <<", "<GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx] <<", "<GetType()] [pCalDigi->GetSm()*fDigiBdfPar->GetNbRpc( pCalDigi->GetType()) + pCalDigi->GetRpc()] [pCalDigi->GetChannel()] [pCalDigi->GetSide()] [iWx+1] <<" -> dWT = "<< dWT ; } } else { LOG(info)<<"Skip1 Digi " <<" Type "<GetType()<<" "<< fDigiBdfPar->GetNbSmTypes() <<" Sm " <GetSm()<<" " << fDigiBdfPar->GetNbSm(pDigi->GetType()) <<" Rpc " <GetRpc()<<" "<< fDigiBdfPar->GetNbRpc(pDigi->GetType()) <<" Ch " <GetChannel()<<" "<GetNbChan(pDigi->GetType(),0) ; } if(pCalDigi->GetType()==5 || pCalDigi->GetType() == 8) { // for Pad counters generate fake digi to mockup a strip CbmTofDigiExp *pCalDigi2 = new((*fTofCalDigisColl)[++iDigIndCal]) CbmTofDigiExp( *pCalDigi ); if(pCalDigi->GetSide()==0) pCalDigi2->SetAddress(pCalDigi->GetSm(),pCalDigi->GetRpc(),pCalDigi->GetChannel(),1,pCalDigi->GetType()); else pCalDigi2->SetAddress(pCalDigi->GetSm(),pCalDigi->GetRpc(),pCalDigi->GetChannel(),0,pCalDigi->GetType());; } } // for( Int_t iDigInd = 0; iDigInd < nTofDigi; iDigInd++ ) iNbTofDigi = fTofCalDigisColl->GetEntries(); // update because of added duplicted digis if(fTofCalDigisColl->IsSortable()) LOG(debug)<<"CbmTofEventClusterizer::BuildClusters: Sort "<GetEntries()<<" calibrated digis " ; if(iNbTofDigi>1){ fTofCalDigisColl->Sort(iNbTofDigi); // Time order again, in case modified by the calibration if(!fTofCalDigisColl->IsSorted()){ LOG(warning)<<"CbmTofEventClusterizer::BuildClusters: Sorting not successful " ; } } return kTRUE; } void CbmTofEventClusterizer::SetDeadStrips (Int_t iDet, Int_t ival){ if( fvDeadStrips.size() < static_cast(iDet+1) ) fvDeadStrips.resize(iDet+1); fvDeadStrips[iDet]=ival; }