// ------------------------------------------------------------------ // ----- CbmTofAnaTestbeam ----- // ----- Created 12/04/2014 by nh ----- // ------------------------------------------------------------------ #include "CbmTofAnaTestbeam.h" // TOF Classes and includes #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 "CbmTofTracklet.h" #include "CbmTofTestBeamClusterizer.h" #include "TMbsMappingTofPar.h" // CBMroot classes and includes // FAIR classes and includes #include "FairRunAna.h" #include "FairRuntimeDb.h" #include "FairRootManager.h" #include "FairLogger.h" // ROOT Classes and includes #include "Riostream.h" #include "TClonesArray.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" #include "TProfile.h" #include "TString.h" #include "TFile.h" #include "TMath.h" #include "TRandom.h" #include "TROOT.h" #include "TDirectory.h" #include "TGeoManager.h" const Int_t DetMask = 4194303; const Double_t DTDMAX=6000.; // diamond inspection range in ps Double_t dTDia; Double_t dDTD4Min=1.E8; //___________________________________________________________________ // // CbmTofAnaTestbeam // // Task for analysis of Testbeam data // // ------------------------------------------------------------------ CbmTofAnaTestbeam::CbmTofAnaTestbeam() : FairTask("HadronAnalysis"), fEvents(0), fGeoHandler(new CbmTofGeoHandler()), fTofId(NULL), fChannelInfo(NULL), fChannelInfoRef(NULL), fChannelInfoDut(NULL), fChannelInfoSel2(NULL), fMbsMappingPar(NULL), iNbSmTot(0), fvTypeSmOffs(), iNbRpcTot(0), fvSmRpcOffs(), iNbChTot(0), fvRpcChOffs(), fDigiPar(NULL), fDigiBdfPar(NULL), fTofDigisColl(NULL), fTofHitsColl(NULL), fTofDigiMatchColl(NULL), fTofTrackColl(NULL), fdDXMean(0.), fdDYMean(0.), fdDTMean(0.), fdDXWidth(0.), fdDYWidth(0.), fdDTWidth(0.), fhDT2(NULL), fhXX2(NULL), fhYY2(NULL), fhNMatch04(NULL), fhXX04(NULL), fhYY04(NULL), fhTT04(NULL), fhChi04(NULL), fhChiSel24(NULL), fhDXSel24(NULL), fhDYSel24(NULL), fhDTSel24(NULL), fhDXDY04(NULL), fhDXDT04(NULL), fhDYDT04(NULL), fhTofSel24(NULL), fhNMatch04sel(NULL), fhChi04best(NULL), fhDigiMul0best(NULL), fhDigiMul4best(NULL), fhDXDY04best(NULL), fhDXDT04best(NULL), fhDYDT04best(NULL), fhChiDT04best(NULL), fhDT24DT04best(NULL), fhDTD4DT04best(NULL), fhX0DT04best(NULL), fhY0DT04best(NULL), fhNMatchD4sel(NULL), fhChi04D4best(NULL), fhTofD4best(NULL), fhVelD4best(NULL), fhDigiMul0D4best(NULL), fhDigiMul4D4best(NULL), fhCluMul04D4best(NULL), fhDXDY04D4best(NULL), fhDXDT04D4best(NULL), fhDYDT04D4best(NULL), fhDistDT04D4best(NULL), fhTexpDT04D4best(NULL), fhCluSize0DT04D4best(NULL), fhCluSize4DT04D4best(NULL), fhTot0DT04D4best(NULL), fhTot4DT04D4best(NULL), fhChiDT04D4best(NULL), fhDT24DT04D4best(NULL), fhDTD4DT04D4best(NULL), fhX0DT04D4best(NULL), fhY0DT04D4best(NULL), fhDTMul4D4best(NULL), fhDTX4D4best(NULL), fhDTY4D4best(NULL), fhDXX4D4best(NULL), fhDXY4D4best(NULL), fhDYX4D4best(NULL), fhDYY4D4best(NULL), fhDTMul0D4best(NULL), fhDTX0D4best(NULL), fhDTY0D4best(NULL), fhDXX0D4best(NULL), fhDXY0D4best(NULL), fhDYX0D4best(NULL), fhDYY0D4best(NULL), fhChi04D4sbest(NULL), fhTofD4sbest(NULL), fhVelD4sbest(NULL), fhDigiMul0D4sbest(NULL), fhDigiMul4D4sbest(NULL), fhCluMul04D4sbest(NULL), fhDXDY04D4sbest(NULL), fhDXDT04D4sbest(NULL), fhDYDT04D4sbest(NULL), fhDistDT04D4sbest(NULL), fhTexpDT04D4sbest(NULL), fhCluSize0DT04D4sbest(NULL), fhCluSize4DT04D4sbest(NULL), fhTot0DT04D4sbest(NULL), fhTot4DT04D4sbest(NULL), fhChiDT04D4sbest(NULL), fhDT24DT04D4sbest(NULL), fhDTD4DT04D4sbest(NULL), fhX0DT04D4sbest(NULL), fhY0DT04D4sbest(NULL), fhDTMul4D4sbest(NULL), fhDTX4D4sbest(NULL), fhDTY4D4sbest(NULL), fhDXX4D4sbest(NULL), fhDXY4D4sbest(NULL), fhDYX4D4sbest(NULL), fhDYY4D4sbest(NULL), fhDTMul0D4sbest(NULL), fhDTX0D4sbest(NULL), fhDTY0D4sbest(NULL), fhDXX0D4sbest(NULL), fhDXY0D4sbest(NULL), fhDYX0D4sbest(NULL), fhDYY0D4sbest(NULL), fhNMatch24(NULL), fhNMatch24sel(NULL), fhDT24sel(NULL), fhChi24(NULL), fhXY24(NULL), fhDXDY24(NULL), fhDXDT24(NULL), fhDYDT24(NULL), fhXY0D4best(NULL), fhXY4D4best(NULL), fhXX04D4best(NULL), fhYY04D4best(NULL), fhXYSel2D4best(NULL), fhXY0D4sel(NULL), fhXY4D4sel(NULL), fhXYSel2D4sel(NULL), fhDTD4sel(NULL), fhTofD4sel(NULL), fhBRefMul(NULL), fhDTD4(NULL), fhXYPos(), fhDT04DX0_1(NULL), fhDT04DY0_1(NULL), fhDT04DT0_1(NULL), fhDT04DX4_1(NULL), fhDT04DY4_1(NULL), fhDT04DT4_1(NULL), fhDT04DX0_2(NULL), fhDT04DY0_2(NULL), fhDT04DT0_2(NULL), fhDT04DX4_2(NULL), fhDT04DY4_2(NULL), fhDT04DT4_2(NULL), fhDutPullX(NULL), fhDutPullY(NULL), fhDutPullZ(NULL), fhDutPullT(NULL), fhDutPullTB(NULL), fhDutChiFound(NULL), fhDutChiMissed(NULL), fhDutChiMatch(NULL), fhDutXY_Found(NULL), fhDutXY_Missed(NULL), fhDutXYDT(NULL), fStart(), fStop(), fCalParFileName(""), fCalParFile(NULL), fhDTD4DT04D4Off(NULL), fhDTX4D4Off(NULL), fhDTY4D4Off(NULL), fhDTTexpD4Off(NULL), fdMulDMax(0.), fdDTDia(0.), fdDTD4MAX(0.), fdMul0Max(0.), fdMul4Max(0.), fdCh4Sel(0.), fdDCh4Sel(0.), fdPosY4Sel(0.), fdPosY4SelOff(0.), fdChS2Sel(0.), fdDChS2Sel(0.), fdPosYS2Sel(0.), fdPosYS2SelOff(0.), fdSel2TOff(0.), fdHitDistMin(0.), fdTOffD4(0.), fdTShift(0.), fdChi2Lim(0.), fiCorMode(0), fiDut(0), fiMrpcRef(0), fiMrpcSel2(0), fiMrpcSel3(0), fiPlaSelect(0), fiBeamRefSmType(0), fiBeamRefSmId(0), fiDutNch(0), fSIGLIM(3.), fSIGT(100.), fSIGX(1.), fSIGY(1.), fFindTracks(NULL) { cout << "CbmTofTests: Task started " << endl; } // ------------------------------------------------------------------ // ------------------------------------------------------------------ CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose) : FairTask(name, verbose), fEvents(0), fGeoHandler(new CbmTofGeoHandler()), fTofId(NULL), fChannelInfo(NULL), fChannelInfoRef(NULL), fChannelInfoDut(NULL), fChannelInfoSel2(NULL), fMbsMappingPar(NULL), iNbSmTot(0), fvTypeSmOffs(), iNbRpcTot(0), fvSmRpcOffs(), iNbChTot(0), fvRpcChOffs(), fDigiPar(NULL), fDigiBdfPar(NULL), fTofDigisColl(NULL), fTofHitsColl(NULL), fTofDigiMatchColl(NULL), fTofTrackColl(NULL), fdDXMean(0.), fdDYMean(0.), fdDTMean(0.), fdDXWidth(0.), fdDYWidth(0.), fdDTWidth(0.), fhDT2(NULL), fhXX2(NULL), fhYY2(NULL), fhNMatch04(NULL), fhXX04(NULL), fhYY04(NULL), fhTT04(NULL), fhChi04(NULL), fhChiSel24(NULL), fhDXSel24(NULL), fhDYSel24(NULL), fhDTSel24(NULL), fhDXDY04(NULL), fhDXDT04(NULL), fhDYDT04(NULL), fhTofSel24(NULL), fhNMatch04sel(NULL), fhChi04best(NULL), fhDigiMul0best(NULL), fhDigiMul4best(NULL), fhDXDY04best(NULL), fhDXDT04best(NULL), fhDYDT04best(NULL), fhChiDT04best(NULL), fhDT24DT04best(NULL), fhDTD4DT04best(NULL), fhX0DT04best(NULL), fhY0DT04best(NULL), fhNMatchD4sel(NULL), fhChi04D4best(NULL), fhTofD4best(NULL), fhVelD4best(NULL), fhDigiMul0D4best(NULL), fhDigiMul4D4best(NULL), fhCluMul04D4best(NULL), fhDXDY04D4best(NULL), fhDXDT04D4best(NULL), fhDYDT04D4best(NULL), fhDistDT04D4best(NULL), fhTexpDT04D4best(NULL), fhCluSize0DT04D4best(NULL), fhCluSize4DT04D4best(NULL), fhTot0DT04D4best(NULL), fhTot4DT04D4best(NULL), fhChiDT04D4best(NULL), fhDT24DT04D4best(NULL), fhDTD4DT04D4best(NULL), fhX0DT04D4best(NULL), fhY0DT04D4best(NULL), fhDTMul4D4best(NULL), fhDTX4D4best(NULL), fhDTY4D4best(NULL), fhDXX4D4best(NULL), fhDXY4D4best(NULL), fhDYX4D4best(NULL), fhDYY4D4best(NULL), fhDTMul0D4best(NULL), fhDTX0D4best(NULL), fhDTY0D4best(NULL), fhDXX0D4best(NULL), fhDXY0D4best(NULL), fhDYX0D4best(NULL), fhDYY0D4best(NULL), fhChi04D4sbest(NULL), fhTofD4sbest(NULL), fhVelD4sbest(NULL), fhDigiMul0D4sbest(NULL), fhDigiMul4D4sbest(NULL), fhCluMul04D4sbest(NULL), fhDXDY04D4sbest(NULL), fhDXDT04D4sbest(NULL), fhDYDT04D4sbest(NULL), fhDistDT04D4sbest(NULL), fhTexpDT04D4sbest(NULL), fhCluSize0DT04D4sbest(NULL), fhCluSize4DT04D4sbest(NULL), fhTot0DT04D4sbest(NULL), fhTot4DT04D4sbest(NULL), fhChiDT04D4sbest(NULL), fhDT24DT04D4sbest(NULL), fhDTD4DT04D4sbest(NULL), fhX0DT04D4sbest(NULL), fhY0DT04D4sbest(NULL), fhDTMul4D4sbest(NULL), fhDTX4D4sbest(NULL), fhDTY4D4sbest(NULL), fhDXX4D4sbest(NULL), fhDXY4D4sbest(NULL), fhDYX4D4sbest(NULL), fhDYY4D4sbest(NULL), fhDTMul0D4sbest(NULL), fhDTX0D4sbest(NULL), fhDTY0D4sbest(NULL), fhDXX0D4sbest(NULL), fhDXY0D4sbest(NULL), fhDYX0D4sbest(NULL), fhDYY0D4sbest(NULL), fhNMatch24(NULL), fhNMatch24sel(NULL), fhDT24sel(NULL), fhChi24(NULL), fhXY24(NULL), fhDXDY24(NULL), fhDXDT24(NULL), fhDYDT24(NULL), fhXY0D4best(NULL), fhXY4D4best(NULL), fhXX04D4best(NULL), fhYY04D4best(NULL), fhXYSel2D4best(NULL), fhXY0D4sel(NULL), fhXY4D4sel(NULL), fhXYSel2D4sel(NULL), fhDTD4sel(NULL), fhTofD4sel(NULL), fhBRefMul(NULL), fhDTD4(NULL), fhXYPos(), fhDT04DX0_1(NULL), fhDT04DY0_1(NULL), fhDT04DT0_1(NULL), fhDT04DX4_1(NULL), fhDT04DY4_1(NULL), fhDT04DT4_1(NULL), fhDT04DX0_2(NULL), fhDT04DY0_2(NULL), fhDT04DT0_2(NULL), fhDT04DX4_2(NULL), fhDT04DY4_2(NULL), fhDT04DT4_2(NULL), fhDutPullX(NULL), fhDutPullY(NULL), fhDutPullZ(NULL), fhDutPullT(NULL), fhDutPullTB(NULL), fhDutChiFound(NULL), fhDutChiMissed(NULL), fhDutChiMatch(NULL), fhDutXY_Found(NULL), fhDutXY_Missed(NULL), fhDutXYDT(NULL), fStart(), fStop(), fCalParFileName(""), fCalParFile(NULL), fhDTD4DT04D4Off(NULL), fhDTX4D4Off(NULL), fhDTY4D4Off(NULL), fhDTTexpD4Off(NULL), fdMulDMax(0.), fdDTDia(0.), fdDTD4MAX(0.), fdMul0Max(0.), fdMul4Max(0.), fdCh4Sel(0.), fdDCh4Sel(0.), fdPosY4Sel(0.), fdPosY4SelOff(0.), fdChS2Sel(0.), fdDChS2Sel(0.), fdPosYS2Sel(0.), fdPosYS2SelOff(0.), fdSel2TOff(0.), fdHitDistMin(0.), fdTOffD4(0.), fdTShift(0.), fdChi2Lim(0.), fiCorMode(0), fiDut(0), fiMrpcRef(0), fiMrpcSel2(0), fiMrpcSel3(0), fiPlaSelect(0), fiBeamRefSmType(0), fiBeamRefSmId(0), fiDutNch(0), fSIGLIM(3.), fSIGT(100.), fSIGX(1.), fSIGY(1.), fFindTracks(NULL) { } // ------------------------------------------------------------------ // ------------------------------------------------------------------ CbmTofAnaTestbeam::~CbmTofAnaTestbeam() { // Destructor } // ------------------------------------------------------------------ /************************************************************************************/ // FairTasks inherited functions InitStatus CbmTofAnaTestbeam::Init() { if( kFALSE == RegisterInputs() ) return kFATAL; // fTofId = new ( CbmTofDetectorId )CbmTofDetectorId_v14a(); if( kFALSE == InitParameters() ) return kFATAL; if( kFALSE == LoadGeometry() ) return kFATAL; if( kFALSE == LoadCalParameter() ) return kFATAL; if( kFALSE == CreateHistos() ) return kFATAL; fFindTracks = CbmTofFindTracks::Instance(); if (NULL == fFindTracks) LOG(WARNING) << Form("CbmTofAnaTestbeam::Init : no FindTracks instance found") << FairLogger::endl; return kSUCCESS; } void CbmTofAnaTestbeam::SetParContainers() { LOG(INFO)<<" CbmTofAnaTestbeam => Get the digi parameters for tof"<GetRuntimeDb(); fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); fDigiBdfPar = (CbmTofDigiBdfPar*) (rtdb->getContainer("CbmTofDigiBdfPar")); } Bool_t CbmTofAnaTestbeam::LoadGeometry() { LOG(INFO)<<"CbmTofAnaTestbeam::LoadGeometry starting for " <GetNrOfModules() << " geometrically known detector cells, looking for type " << fiDut <GetNrOfModules(); ++icell) { Int_t cellId = fDigiPar->GetCellId(icell); // cellId is assigned in CbmTofCreateDigiPar fChannelInfo = fDigiPar->GetCell(cellId); /* Int_t smodule = fGeoHandler->GetSMType(cellId); // FIXME - wrong inline functions!!! Int_t smtype = fGeoHandler->GetSModule(cellId); // FIXME Int_t module = fGeoHandler->GetCounter(cellId); */ Int_t smtype = CbmTofAddress::GetSmType( cellId ); Int_t smodule = CbmTofAddress::GetSmId( cellId ); Int_t module = CbmTofAddress::GetRpcId( cellId ); LOG(DEBUG) <GetSizex() << ", y-size "<< fChannelInfo->GetSizey() << FairLogger::endl; if( smtype == fiDut ){ fiDutNch++; if(fChannelInfo->GetX() > dDutXmax) dDutXmax=fChannelInfo->GetX(); if(fChannelInfo->GetX() < dDutXmin) dDutXmin=fChannelInfo->GetX(); if(fiDutNch == 1) fChannelInfoDut = fChannelInfo; } } if (fiDutNch > 0){ Double_t dDutDx = (dDutXmax - dDutXmin)/(fiDutNch-1); LOG(INFO)<<"CbmTofAnaTestbeam::LoadGeometry Dut = " << fiDut << " with " << fiDutNch << " channels in x- direction from " << dDutXmin << " to " << dDutXmax <<", dx = "<< dDutDx < New event"<GetObject("CbmTofDigiExp"); if( NULL == fTofDigisColl) fTofDigisColl = (TClonesArray *) fManager->GetObject("CbmTofDigi"); if( NULL == fTofDigisColl) { LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!!"<GetObject("TofHit"); if( NULL == fTofHitsColl) { LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofHit TClonesArray!!!"<GetObject("TofDigiMatch"); if( NULL == fTofDigiMatchColl) { LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the Match TClonesArray!!!"<GetObject("TofTracks"); if( NULL == fTofTrackColl) { LOG(ERROR)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofTracklet TClonesArray!!!"<GetRuntimeDb(); Int_t iGeoVersion = fGeoHandler->Init(isSimulation); if( k14a > iGeoVersion ) { } fDigiPar = (CbmTofDigiPar*) (rtdb->getContainer("CbmTofDigiPar")); if( 0 == fDigiPar ) { LOG(ERROR)<<"CbmTofAnaTestbeam::InitParameters => Could not obtain the CbmTofDigiPar "<getContainer("CbmTofDigiBdfPar")); if( 0 == fDigiBdfPar ) { LOG(ERROR)<<"CbmTofAnaTestbeam::InitParameters => Could not obtain the CbmTofDigiBdfPar "<getContainer("TMbsMappingTofPar")); if( 0 == fMbsMappingPar ) { LOG(ERROR)<<"CbmTofAnaTestBeam::InitParameters => Could not obtain the TMbsMappingTofPar "<initContainers( ana->GetRunId() ); LOG(INFO)<<"CbmTofAnaTestbeam::InitParameter: currently " << fDigiPar->GetNrOfModules() << " digi cells " <FindObjectAny( Form("hDTD4DT04D4best_pfx_px")); if (NULL == fhtmp) { LOG(INFO)<<" Histo " << Form("hDTD4DT04D4best_pfx_px") << " not found. " <FindObjectAny( Form("hDTX4D4best_pfx_px")); if (NULL == fhtmpx) { LOG(INFO)<<" Histo " << Form("hDTX4D4best_pfx_px") << " not found. " <FindObjectAny( Form("hDTY4D4best_pfx_px")); if (NULL == fhtmpy) { LOG(INFO)<<" Histo " << Form("hDTY4D4best_pfx_px") << " not found. " <FindObjectAny( Form("hTexpDT04D4best_pfx_px")); if (NULL == fhtmpt) { LOG(INFO)<<" Histo " << Form("hTexpDT04D4best_pfx_px") << " not found. " <cd(); if(NULL != fhtmp) fhDTD4DT04D4Off=(TH1D *)fhtmp->Clone(); if(NULL != fhtmpx) fhDTX4D4Off=(TH1D *)fhtmpx->Clone(); if(NULL != fhtmpy) fhDTY4D4Off=(TH1D *)fhtmpy->Clone(); if(NULL != fhtmpt) fhDTTexpD4Off=(TH1D *)fhtmpt->Clone(); fCalParFile->Close(); // fhDTD4DT04D4Off->Draw(); if (fdDCh4Sel==0.) fdDCh4Sel=1000.; // open default window if (fdPosY4Sel==0.) fdPosY4Sel=10.; // open default window return kTRUE; } /************************************************************************************/ // ------------------------------------------------------------------ Bool_t CbmTofAnaTestbeam::CreateHistos() { // Create histogramms Double_t XDMAX=100.; Double_t YDMAX=100.; Double_t TDMAX=100000.; Double_t DTMAX=1000.; Double_t DXMAX=10.; Double_t DYMAX=10.; /* Double_t XMAX=100.;*/ /* Double_t YMAX=100.;*/ 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 ! // define histos here fhBRefMul = new TH1F( Form("hBRefMul"),Form("Multiplicity in Beam Reference counter ; Mul ()"), 50, 0., 50.); fhDTD4 = new TH1F( Form("hDTD4"),Form("reference time ; #Delta tD4 (ps)"), 100, -100000., 100000.); Int_t iNbDet=fMbsMappingPar->GetNbMappedDet(); fhXYPos.resize( iNbDet ); for(Int_t iDet=0; iDetcd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager! return kTRUE; } // ------------------------------------------------------------------ Bool_t CbmTofAnaTestbeam::FillHistos() { // Constants, TODO => put as parameter !!! /* Int_t kTOF=6;*/ // Declare variables outside the loop CbmTofHit *pHit; CbmTofHit *pHit1; CbmTofHit *pHit2; CbmTofHit *pHit3; CbmTofHit *pHit4; CbmTofHit *pHitRef =NULL; // May be used uninitialized later, to check! CbmTofHit *pHitSel2=NULL; CbmTofHit *pDia; CbmTofCell *fChannelInfo1; CbmTofCell *fChannelInfo2; CbmTofCell *fChannelInfo3; CbmTofCell *fChannelInfo4; /* Int_t iNbTofDigis;*/ Int_t iNbTofHits, iNbTofTracks; // iNbTofDigis = fTofDigisColl->GetEntriesFast(); iNbTofHits = fTofHitsColl->GetEntriesFast(); /* LOG(INFO)<UseExpandedDigi() ) { CbmTofDigiExp *pDigi; for( Int_t iDigInd = 0; iDigInd < iNbTofDigis; iDigInd++ ) { pDigi = (CbmTofDigiExp*) fTofDigisColl->At( iDigInd ); Int_t iSmType = pDigi->GetType(); Int_t iSm = pDigi->GetSm(); Int_t iRpc = pDigi->GetRpc(); Int_t iCh = pDigi->GetChannel(); } // for( Int_t iDigInd = 0; iDigInd < iNbTofDigis; iDigInd++ ) } // if( kTRUE == fDigiBdfPar->UseExpandedDigi() ) */ // Hits info Int_t iNbMatchedHits = 0; Int_t iNbMaxMatch=100; Double_t Zref=300.; /* Double_t Chi2MatchMin=1.E8;*/ Double_t Chi2List[iNbMaxMatch]; CbmTofHit *pChi2Hit1[iNbMaxMatch]; CbmTofHit *pChi2Hit2[iNbMaxMatch]; Chi2List[0]=1.E8; pChi2Hit1[0]=NULL; pChi2Hit2[0]=NULL; pDia=NULL; Int_t iNSel=1; Bool_t BSel[iNSel]; for(Int_t iSel=0;iSel0) dM4Max=fdMul4Max; Double_t dM0Max=100; // modify if (fdMul0Max>0) dM0Max=fdMul0Max; Double_t dMDMax=1; if (fdMulDMax>0) dMDMax=fdMulDMax; Double_t hitpos1[3], hitpos2[3], hitpos3[3], hitpos4[3]; Double_t hitpos1_local[3], hitpos2_local[3], hitpos3_local[3], hitpos4_local[3]; // find diamond reference for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & DetMask); Int_t iChId = pHit->GetAddress(); fChannelInfo = fDigiPar->GetCell( iChId ); Int_t iSmType=CbmTofAddress::GetSmType( iDetId ); if(NULL == fChannelInfo){ LOG(DEBUG) << "CbmTofAnaTestbeam::FillHistos: NULL Channel Pointer for ChId " << Form(" 0x%08x ",iChId) <GetTime()+fdTShift; pHit->SetTime(dTime); // shift beam reference times if(pHit->GetTime() != dTime) LOG(ERROR)<<" Hit time not updated "<GetTime()<<" -> "<GetTime() < dTDia) { dTDia = pHit->GetTime(); pDia = pHit; } } } } // reaction loop end; fhBRefMul->Fill(dMulD); LOG(DEBUG)<At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & DetMask); Int_t iChId = pHit->GetAddress(); fChannelInfo = fDigiPar->GetCell( iChId ); Int_t iSmType=CbmTofAddress::GetSmType( iDetId ); LOG(DEBUG)<GetX(), pHit->GetY(), pHit->GetZ(), pHit->GetTime()) <(iSmType)Fill(pHit->GetX(),pHit->GetY()); if(fiDut == iSmType) { dMul0++; Double_t xPos1=Zref/pHit->GetZ()*pHit->GetX(); Double_t yPos1=Zref/pHit->GetZ()*pHit->GetY(); Double_t tof1 =pHit->GetTime(); for( Int_t iHitInd2 = 0; iHitInd2 < iNbTofHits; iHitInd2++) { if(iHitInd2 != iHitInd){ pHit2 = (CbmTofHit*) fTofHitsColl->At( iHitInd2 ); if (NULL== pHit2) continue; Int_t iDetId2 = (pHit2->GetAddress() & DetMask); Int_t iChId2 = pHit2->GetAddress(); fChannelInfo2 = fDigiPar->GetCell( iChId2 ); if (fiMrpcRef == CbmTofAddress::GetSmType( iDetId2 )){ // MrpcRef Double_t xPos2=Zref/pHit2->GetZ()*pHit2->GetX(); Double_t yPos2=Zref/pHit2->GetZ()*pHit2->GetY(); Double_t tof2 =pHit2->GetTime(); Double_t dTcor=0.; if(fhDTD4DT04D4Off != NULL) dTcor=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(dTDia-tof2)); Double_t Chi2Match =TMath::Power((xPos1-xPos2-fdDXMean)/fdDXWidth,2.) +TMath::Power((yPos1-yPos2-fdDYMean)/fdDYWidth,2.) +TMath::Power((tof1-tof2-dTcor-fdDTMean)/fdDTWidth,2.); if (Chi2Match > 1.E8) continue; Chi2Match /= 3; LOG(DEBUG1)<<" Chi2 "< %f ", Chi2Match) < %f ", Chi2Match) <Fill(xPos1,xPos2); fhYY04->Fill(yPos1,yPos2); fhTT04->Fill(tof1,tof2); fhDXDY04->Fill(xPos1-xPos2,yPos1-yPos2); fhDXDT04->Fill(xPos1-xPos2,tof1-tof2); fhDYDT04->Fill(yPos1-yPos2,tof1-tof2); fhChi04->Fill(Chi2Match); for(Int_t iM=0; iM=iM; iMM--){ Chi2List[iMM]= Chi2List[iMM-1]; pChi2Hit1[iMM]= pChi2Hit1[iMM-1]; pChi2Hit2[iMM]= pChi2Hit2[iMM-1]; } Chi2List[iM]=Chi2Match; pChi2Hit1[iM]=pHit; pChi2Hit2[iM]=pHit2; Chi2List[iNbMatchedHits]=1.E8; if(iM>0){ if(Chi2Match == Chi2List[iM-1]){ LOG(DEBUG)<ToString() <ToString()<ToString() <Fill(Zref/pHit->GetZ()*pHit->GetX(),Zref/pHit2->GetZ()*pHit2->GetX()); fhYY02[iSm]->Fill(Zref/pHit->GetZ()*pHit->GetY(),Zref/pHit2->GetZ()*pHit2->GetY()); } } } //iHit2 != iHit condition end } // iHit2 loop end } // fiDut condition end if(fiMrpcRef == iSmType) { dMul4++; } if(fiBeamRefSmType == iSmType) { if( pHit != pDia ) continue; if( fdDTDia>0. ) { Double_t dDDia=0.; for( Int_t iHitInd1 = 0; iHitInd1 < iNbTofHits; iHitInd1++) if( iHitInd1!=iHitInd) { pHit1 = (CbmTofHit*) fTofHitsColl->At( iHitInd1 ); if(pHit1==NULL) continue; Int_t iDetId1 = (pHit1->GetAddress() & DetMask); Int_t iChId1 = pHit1->GetAddress(); fChannelInfo1 = fDigiPar->GetCell( iChId1 ); pHit1 = (CbmTofHit*) fTofHitsColl->At( iHitInd1 ); if( fiBeamRefSmType == CbmTofAddress::GetSmType( iDetId1 ) && fiBeamRefSmId != CbmTofAddress::GetSmId ( iDetId1 ) && TMath::Abs( pHit1->GetTime()-dTDia)GetTime()-dTDia; } } LOG(DEBUG)<At( iHitInd2 ); if(pHit2==NULL) continue; Int_t iDetId2 = (pHit2->GetAddress() & DetMask); Int_t iChId2 = pHit2->GetAddress(); fChannelInfo2 = fDigiPar->GetCell( iChId2 ); if(NULL == fChannelInfo2){ LOG(DEBUG) << "CbmTofAnaTestbeam::FillHistos: Invalid Channel Pointer for ChId2 " << Form(" 0x%08x ",iChId2) <GetTime()) <GetTime()-dTDia + fdTShift; fhDTD4->Fill(dDTD4); LOG(DEBUG1)<local trafo gGeoManager->FindNode(fChannelInfo2->GetX(),fChannelInfo2->GetY(),fChannelInfo2->GetZ()); Double_t hitpos[3], hitpos_local[3]; hitpos[0]=pHit2->GetX(); hitpos[1]=pHit2->GetY(); hitpos[2]=pHit2->GetZ(); /*TGeoNode* cNode=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos, hitpos_local); if( TMath::Abs(hitpos_local[1]-fdPosY4SelOff)GetSizey() &&TMath::Abs(dDTD4)GetZ()*pHit2->GetX(); Double_t yPos2=Zref/pHit2->GetZ()*pHit2->GetY(); Double_t tof2 =pHit2->GetTime(); Double_t dTcor=0.; // if(fhDTD4DT04D4Off != NULL) // dTcor=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(dTDia-tof2)); Double_t Chi2Max = fdChi2Lim; pHitSel2 = NULL; for( Int_t iHitInd3 = 0; iHitInd3 < iNbTofHits; iHitInd3++) if(iHitInd3 != iHitInd && iHitInd3 != iHitInd2) { pHit3 = (CbmTofHit*) fTofHitsColl->At( iHitInd3 ); if(pHit3==NULL) continue; Int_t iDetId3 = (pHit3->GetAddress() & DetMask); Int_t iChId3 = pHit3->GetAddress(); fChannelInfo3 = fDigiPar->GetCell( iChId3 ); if(NULL == fChannelInfo3){ LOG(DEBUG) << "CbmTofAnaTestbeam::FillHistos: Invalid Channel Pointer for ChId3 " << Form(" 0x%08x ",iChId3) <local trafo gGeoManager->FindNode(fChannelInfo3->GetX(),fChannelInfo3->GetY(),fChannelInfo3->GetZ()); hitpos3[0]=pHit3->GetX(); hitpos3[1]=pHit3->GetY(); hitpos3[2]=pHit3->GetZ(); /*TGeoNode* cNode3=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos3, hitpos3_local); if( TMath::Abs(hitpos3_local[1]-fdPosYS2SelOff)GetSizey() ){ Double_t xPos3=Zref/pHit3->GetZ()*pHit3->GetX(); Double_t yPos3=Zref/pHit3->GetZ()*pHit3->GetY(); Double_t tof3 =pHit3->GetTime(); Double_t Chi2Match =TMath::Power((xPos3-xPos2)/fdDXWidth,2.) +TMath::Power((yPos3-yPos2)/fdDYWidth,2.) +TMath::Power((tof3-tof2-dTcor-fdSel2TOff)/fdDTWidth,2.); Chi2Match /= 3; if (Chi2Match < Chi2Max) { fhChiSel24->Fill(Chi2Match); fhDXSel24->Fill(xPos3-xPos2); fhDYSel24->Fill(yPos3-yPos2); fhDTSel24->Fill(tof3-tof2-dTcor-fdSel2TOff); fhTofSel24->Fill(tof3-tof2-fdSel2TOff); Chi2Max = Chi2Match; BSel[0] = kTRUE; pHitSel2= pHit3; fChannelInfoSel2=fChannelInfo3; } } } } } // loop over third hit } } } } } //LOG(INFO)<<" TDia="<dM4Max || dMulD>dMDMax || dMul0>dM0Max) { BSel[0]=kFALSE; LOG(DEBUG) << Form(" Muls %4.0f, %4.0f, %4.0f, Matches %d",dMulD, dMul0, dMul4, iNbMatchedHits) << FairLogger::endl; } // normalisation distributions fhNMatch04->Fill(iNbMatchedHits); LOG(DEBUG)<Fill(iNbMatchedHits); // use as normalisation fhTofD4sel->Fill(pHitRef->GetTime()-pDia->GetTime()); // general normalisation fhDTD4sel->Fill(dDTD4Min); // general normalisation if(fChannelInfoDut != NULL){ // Project into Dut reference frame /*TGeoNode *fNodeDut=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfoDut->GetX(),fChannelInfoDut->GetY(),fChannelInfoDut->GetZ()); hitpos1[0]=fChannelInfoDut->GetZ()/pHitRef->GetZ()*pHitRef->GetX(); hitpos1[1]=fChannelInfoDut->GetZ()/pHitRef->GetZ()*pHitRef->GetY(); hitpos1[2]=fChannelInfoDut->GetZ(); /*TGeoNode* cNodeDut=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos1, hitpos1_local); //hitpos1_local[0] -= fiDutNch/2 * fChannelInfoDut->GetSizex(); fhXY0D4sel->Fill(hitpos1_local[0],hitpos1_local[1]); } // Monitor selected Reference Hit position /*TGeoNode *fNodeRef=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfoRef->GetX(),fChannelInfoRef->GetY(),fChannelInfoRef->GetZ()); hitpos2[0]=pHitRef->GetX(); hitpos2[1]=pHitRef->GetY(); hitpos2[2]=pHitRef->GetZ(); /*TGeoNode* cNodeRef=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos2, hitpos2_local); fhXY4D4sel->Fill(hitpos2_local[0],hitpos2_local[1]); if(NULL != pHitSel2){ /*Int_t iDetId3 = (pHitSel2->GetAddress() & DetMask);*/ Int_t iChId3 = pHitSel2->GetAddress(); fChannelInfo3 = fDigiPar->GetCell( iChId3 ); /*TGeoNode *fNode3=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo3->GetX(),fChannelInfo3->GetY(),fChannelInfo3->GetZ()); hitpos3[0]=pHitSel2->GetX(); hitpos3[1]=pHitSel2->GetY(); hitpos3[2]=pHitSel2->GetZ(); /*TGeoNode* cNode3=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos3, hitpos3_local); fhXYSel2D4sel->Fill(hitpos3_local[0],hitpos3_local[1]); } if(iNbMatchedHits>0){ // best match LOG(DEBUG)<GetAddress(),pChi2Hit2[0]->GetAddress(), Chi2List[0]) <Fill(hitpos3_local[0],hitpos3_local[1]); } pHit1=pChi2Hit1[0]; pHit2=pChi2Hit2[0]; Int_t iM0=0; if (pHit2 != pHitRef ) { LOG(DEBUG)<<" selector hit does not match reference hit for best match, chi2best " << Chi2List[0] << FairLogger::endl; for (iM0=1; iM0GetAddress() & DetMask); Int_t iChId1 = pHit1->GetAddress(); fChannelInfo1 = fDigiPar->GetCell( iChId1 ); Int_t iDetId2 = (pHit2->GetAddress() & DetMask); Int_t iChId2 = pHit2->GetAddress(); fChannelInfo2 = fDigiPar->GetCell( iChId2 ); Double_t xPos1=Zref/pHit1->GetZ()*pHit1->GetX(); Double_t yPos1=Zref/pHit1->GetZ()*pHit1->GetY(); Double_t tof1=pHit1->GetTime(); Double_t xPos2=Zref/pHit2->GetZ()*pHit2->GetX(); Double_t yPos2=Zref/pHit2->GetZ()*pHit2->GetY(); Double_t tof2=pHit2->GetTime(); Double_t dDist=TMath::Sqrt( TMath::Power(pHit1->GetX()-pHit2->GetX(),2) +TMath::Power(pHit1->GetY()-pHit2->GetY(),2) +TMath::Power(pHit1->GetZ()-pHit2->GetZ(),2) ); fhChi04D4best->Fill(Chi2List[iM0]); CbmMatch* digiMatch1=(CbmMatch *)fTofDigiMatchColl->At(pHit1->GetRefId()); fhDigiMul0D4best->Fill(digiMatch1->GetNofLinks()/2.); CbmMatch* digiMatch2=(CbmMatch *)fTofDigiMatchColl->At(pHit2->GetRefId()); fhDigiMul4D4best->Fill(digiMatch2->GetNofLinks()/2.); fhCluMul04D4best->Fill(dMul0,dMul4); // check for dependence in counter reference frame /*TGeoNode *fNode=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo2->GetX(),fChannelInfo2->GetY(),fChannelInfo2->GetZ()); hitpos2[0]=pHit2->GetX(); hitpos2[1]=pHit2->GetY(); hitpos2[2]=pHit2->GetZ(); /* TGeoNode* cNode=*/ gGeoManager->GetCurrentNode(); // -> Comment to remove warning because set but never used // if(0) cNode->Print(); gGeoManager->MasterToLocal(hitpos2, hitpos2_local); Double_t dTofD4 = fdTOffD4 + dDTD4Min; Double_t dInvVel = dTofD4/pHitRef->GetR(); // in ps/cm Double_t dDTexp = dDist*dInvVel; Double_t dTcor=0.; if(fhDTD4DT04D4Off != NULL) dTcor=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(-dDTD4Min)); if(fhDTX4D4Off != NULL) dTcor+=(Double_t)fhDTX4D4Off->GetBinContent(fhDTX4D4Off->FindBin(hitpos2_local[0])); if(fhDTY4D4Off != NULL) dTcor+=(Double_t)fhDTY4D4Off->GetBinContent(fhDTY4D4Off->FindBin(hitpos2_local[1])); if(fhDTTexpD4Off != NULL) dTcor+=(Double_t)fhDTTexpD4Off->GetBinContent(fhDTTexpD4Off->FindBin(dDTexp)); // LOG(INFO) << "dTcor for "<<-dDTD4<<" from "<Fill(dTofD4); if(dInvVel>0.) fhVelD4best->Fill(1000./dInvVel); fhChiDT04D4best->Fill(Chi2List[iM0],tof1-tof2-dTcor); fhDTD4DT04D4best->Fill(-dDTD4Min,tof1-tof2-dTcor); fhDTMul4D4best->Fill(dMul4,tof1-tof2-dTcor); fhXY4D4best->Fill(hitpos2_local[0],hitpos2_local[1]); fhDTX4D4best->Fill(hitpos2_local[0],tof1-tof2-dTcor); fhDTY4D4best->Fill(hitpos2_local[1],tof1-tof2-dTcor); fhDXX4D4best->Fill(hitpos2_local[0],xPos1-xPos2); fhDXY4D4best->Fill(hitpos2_local[1],xPos1-xPos2); fhDYX4D4best->Fill(hitpos2_local[0],yPos1-yPos2); fhDYY4D4best->Fill(hitpos2_local[1],yPos1-yPos2); fhCluSize4DT04D4best->Fill(digiMatch2->GetNofLinks()/2.,tof1-tof2-dTcor); Double_t dTot = 0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink++){ // loop over digis CbmLink L0 = digiMatch2->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0)); dTot += pDig0->GetTot(); LOG(DEBUG)< sum %f",iDetId2,iDigInd0,pDig0->GetTot(),dTot) <GetNofLinks(); // average time over threshold fhTot4DT04D4best->Fill(dTot,tof1-tof2-dTcor); CbmMatch* digiMatch0=(CbmMatch *)fTofDigiMatchColl->At(pHit1->GetRefId()); fhCluSize0DT04D4best->Fill(digiMatch0->GetNofLinks()/2.,tof1-tof2-dTcor); dTot = 0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink++){ // loop over digis CbmLink L0 = digiMatch0->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0)); dTot += pDig0->GetTot(); LOG(DEBUG1)< sum %f",iDetId1,iDigInd0,pDig0->GetTot(),dTot) <GetNofLinks(); // average time over threshold fhTot0DT04D4best->Fill(dTot,tof1-tof2-dTcor); fhDTMul0D4best->Fill(dMul0,tof1-tof2-dTcor); // check for dependence in counter reference frame /*TGeoNode *fNode1=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo1->GetX(),fChannelInfo1->GetY(),fChannelInfo1->GetZ()); hitpos1[0]=pHit1->GetX(); hitpos1[1]=pHit1->GetY(); hitpos1[2]=pHit1->GetZ(); /*TGeoNode* cNode1=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos1, hitpos1_local); fhXY0D4best->Fill(hitpos1_local[0],hitpos1_local[1]); fhXX04D4best->Fill(hitpos1_local[0],hitpos2_local[0]); fhYY04D4best->Fill(hitpos1_local[1],hitpos2_local[1]); fhDTX0D4best->Fill(hitpos1_local[0],tof1-tof2-dTcor); fhDTY0D4best->Fill(hitpos1_local[1],tof1-tof2-dTcor); fhDXX0D4best->Fill(hitpos1_local[0],xPos1-xPos2); fhDXY0D4best->Fill(hitpos1_local[1],xPos1-xPos2); fhDYX0D4best->Fill(hitpos1_local[0],yPos1-yPos2); fhDYY0D4best->Fill(hitpos1_local[1],yPos1-yPos2); fhDXDY04D4best->Fill(xPos1-xPos2,yPos1-yPos2); fhDXDT04D4best->Fill(xPos1-xPos2,tof1-tof2-dTcor); fhDYDT04D4best->Fill(yPos1-yPos2,tof1-tof2-dTcor); fhDistDT04D4best->Fill(dDist,tof1-tof2-dTcor); fhTexpDT04D4best->Fill(dDTexp,tof1-tof2-dTcor); fhX0DT04D4best->Fill(hitpos1_local[0],tof1-tof2-dTcor); fhY0DT04D4best->Fill(hitpos1_local[1],tof1-tof2-dTcor); if(iNbMatchedHits>1){ LOG(DEBUG)<1: %d with first chi2s = %12.1f, %12.1f, %12.1f, %12.1f",iNbMatchedHits, Chi2List[0],Chi2List[1],Chi2List[2],Chi2List[3]) <GetAddress(),pHit1->GetAddress(), pChi2Hit2[iM]->GetAddress(),pHit2->GetAddress() ) < 1.E3) break; // FIXME hardwired limit ! pHit3=pChi2Hit1[iM]; pHit4=pChi2Hit2[iM]; /*Int_t iDetId3 = (pHit1->GetAddress() & DetMask);*/ Int_t iChId3 = pHit1->GetAddress(); fChannelInfo3 = fDigiPar->GetCell( iChId3 ); /*Int_t iDetId4 = (pHit4->GetAddress() & DetMask);*/ Int_t iChId4 = pHit4->GetAddress(); fChannelInfo4 = fDigiPar->GetCell( iChId4 ); // check for dependence in counter reference frame /*TGeoNode *fNode4= */ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo4->GetX(),fChannelInfo4->GetY(),fChannelInfo4->GetZ()); hitpos4[0]=pChi2Hit2[iM]->GetX(); hitpos4[1]=pChi2Hit2[iM]->GetY(); hitpos4[2]=pChi2Hit2[iM]->GetZ(); /* cNode=*/ gGeoManager->GetCurrentNode(); // -> Comment to remove warning because set but never used gGeoManager->MasterToLocal(hitpos4, hitpos4_local); if(TMath::Abs(hitpos4_local[1])>fdPosY4Sel*fChannelInfo4->GetSizey()) continue; fhChi04D4sbest->Fill(Chi2List[iM]); Double_t xPos3=Zref/pChi2Hit1[iM]->GetZ()*pChi2Hit1[iM]->GetX(); Double_t yPos3=Zref/pChi2Hit1[iM]->GetZ()*pChi2Hit1[iM]->GetY(); Double_t tof3=pChi2Hit1[iM]->GetTime(); Double_t xPos4=Zref/pChi2Hit2[iM]->GetZ()*pChi2Hit2[iM]->GetX(); Double_t yPos4=Zref/pChi2Hit2[iM]->GetZ()*pChi2Hit2[iM]->GetY(); Double_t tof4=pChi2Hit2[iM]->GetTime(); Double_t dDist34=TMath::Sqrt( TMath::Power(pHit3->GetX()-pHit4->GetX(),2) +TMath::Power(pHit3->GetY()-pHit4->GetY(),2) +TMath::Power(pHit3->GetZ()-pHit4->GetZ(),2) ); CbmMatch* digiMatch3=(CbmMatch *)fTofDigiMatchColl->At(pHit3->GetRefId()); fhDigiMul0D4sbest->Fill(digiMatch3->GetNofLinks()/2.); CbmMatch* digiMatch4=(CbmMatch *)fTofDigiMatchColl->At(pHit4->GetRefId()); fhDigiMul4D4sbest->Fill(digiMatch4->GetNofLinks()/2.); fhCluMul04D4sbest->Fill(dMul0,dMul4); Double_t dTofD44 = fdTOffD4 + pHit4->GetTime()-dTDia; Double_t dInvVel4 = dTofD44/pHitRef->GetR(); // in ps/cm Double_t dDTexp4 = dDist34*dInvVel; Double_t dTcor4=0.; if(fhDTD4DT04D4Off != NULL) dTcor4=(Double_t)fhDTD4DT04D4Off->GetBinContent(fhDTD4DT04D4Off->FindBin(dTDia-pHit4->GetTime())); if(fhDTX4D4Off != NULL) dTcor4+=(Double_t)fhDTX4D4Off->GetBinContent(fhDTX4D4Off->FindBin(hitpos4_local[0])); if(fhDTY4D4Off != NULL) dTcor4+=(Double_t)fhDTY4D4Off->GetBinContent(fhDTY4D4Off->FindBin(hitpos4_local[1])); if(fhDTTexpD4Off != NULL) dTcor4+=(Double_t)fhDTTexpD4Off->GetBinContent(fhDTTexpD4Off->FindBin(dDTexp4)); fhTofD4sbest->Fill(dTofD44); if(dInvVel>0.) fhVelD4sbest->Fill(1000./dInvVel4); fhChiDT04D4sbest->Fill(Chi2List[iM],tof3-tof4-dTcor4); fhDTD4DT04D4sbest->Fill(dTDia-pHit4->GetTime(),tof3-tof4-dTcor4); fhDTMul4D4sbest->Fill(dMul4,tof3-tof4-dTcor4); fhDTX4D4sbest->Fill(hitpos4_local[0],tof3-tof4-dTcor4); fhDTY4D4sbest->Fill(hitpos4_local[1],tof3-tof4-dTcor4); fhDXX4D4sbest->Fill(hitpos4_local[0],xPos3-xPos4); fhDXY4D4sbest->Fill(hitpos4_local[1],xPos3-xPos4); fhDYX4D4sbest->Fill(hitpos4_local[0],yPos3-yPos4); fhDYY4D4sbest->Fill(hitpos4_local[1],yPos3-yPos4); fhCluSize4DT04D4sbest->Fill(digiMatch4->GetNofLinks()/2.,tof3-tof4-dTcor4); Double_t dTot4 = 0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink++){ // loop over digis CbmLink L0 = digiMatch4->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0)); dTot4 += pDig0->GetTot(); LOG(DEBUG)< sum %f",iDetId2,iDigInd0,pDig0->GetTot(),dTot4) <GetNofLinks(); // average time over threshold fhTot4DT04D4sbest->Fill(dTot4,tof3-tof4-dTcor4); fhCluSize0DT04D4sbest->Fill(digiMatch3->GetNofLinks()/2.,tof3-tof4-dTcor4); Double_t dTot3 = 0.; for (Int_t iLink=0; iLinkGetNofLinks(); iLink++){ // loop over digis CbmLink L0 = digiMatch3->GetLink(iLink); Int_t iDigInd0=L0.GetIndex(); if (iDigInd0 < fTofDigisColl->GetEntries()){ CbmTofDigiExp *pDig0 = (CbmTofDigiExp*) (fTofDigisColl->At(iDigInd0)); dTot3 += pDig0->GetTot(); LOG(DEBUG)< sum %f",iDetId1,iDigInd0,pDig0->GetTot(),dTot3) <GetNofLinks(); // average time over threshold fhTot0DT04D4sbest->Fill(dTot3,tof3-tof4-dTcor4); fhDTMul0D4sbest->Fill(dMul0,tof3-tof4-dTcor4); // check for dependence in counter reference frame /*TGeoNode *fNode3=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo3->GetX(),fChannelInfo3->GetY(),fChannelInfo3->GetZ()); hitpos3[0]=pChi2Hit1[iM]->GetX(); hitpos3[1]=pChi2Hit1[iM]->GetY(); hitpos3[2]=pChi2Hit1[iM]->GetZ(); /*TGeoNode* cNode3=*/ gGeoManager->GetCurrentNode(); gGeoManager->MasterToLocal(hitpos3, hitpos3_local); fhDTX0D4sbest->Fill(hitpos3_local[0],tof3-tof4-dTcor4); fhDTY0D4sbest->Fill(hitpos3_local[1],tof4-tof4-dTcor4); fhDXX0D4sbest->Fill(hitpos3_local[0],xPos3-xPos4); fhDXY0D4sbest->Fill(hitpos3_local[1],xPos3-xPos4); fhDYX0D4sbest->Fill(hitpos3_local[0],yPos3-yPos4); fhDYY0D4sbest->Fill(hitpos3_local[1],yPos3-yPos4); fhDXDY04D4sbest->Fill(xPos3-xPos4,yPos3-yPos4); fhDXDT04D4sbest->Fill(xPos3-xPos4,tof3-tof4-dTcor4); fhDYDT04D4sbest->Fill(yPos3-yPos4,tof3-tof4-dTcor4); fhDistDT04D4sbest->Fill(dDist34,tof3-tof4-dTcor4); fhTexpDT04D4sbest->Fill(dDTexp4,tof3-tof4-dTcor4); fhX0DT04D4sbest->Fill(hitpos3_local[0],tof3-tof4-dTcor4); fhY0DT04D4sbest->Fill(hitpos3_local[1],tof3-tof4-dTcor4); fhDT04DX0_2->Fill(hitpos1_local[0]-hitpos3_local[0],tof3-tof4-dTcor4); fhDT04DY0_2->Fill(hitpos1_local[1]-hitpos3_local[1],tof3-tof4-dTcor4); fhDT04DT0_2->Fill(tof1-tof3,tof3-tof4-dTcor4); fhDT04DX4_2->Fill(hitpos2_local[0]-hitpos4_local[0],tof3-tof4-dTcor4); fhDT04DY4_2->Fill(hitpos2_local[1]-hitpos4_local[1],tof3-tof4-dTcor4); fhDT04DT4_2->Fill(tof2-tof4,tof3-tof4-dTcor4); fhDT04DX0_1->Fill(hitpos1_local[0]-hitpos3_local[0],tof1-tof2-dTcor); fhDT04DY0_1->Fill(hitpos1_local[1]-hitpos3_local[1],tof1-tof2-dTcor); fhDT04DT0_1->Fill(tof1-tof3,tof1-tof2-dTcor); fhDT04DX4_1->Fill(hitpos2_local[0]-hitpos4_local[0],tof1-tof2-dTcor); fhDT04DY4_1->Fill(hitpos2_local[1]-hitpos4_local[1],tof1-tof2-dTcor); fhDT04DT4_1->Fill(tof2-tof4,tof1-tof2-dTcor); break; } } } } // end of if(iNbMatchedHits>0) } // Tracklet based analysis if(NULL!=fTofTrackColl){ iNbTofTracks = fTofTrackColl->GetEntries(); Int_t NStations = fFindTracks->GetNStations(); LOG(DEBUG)<0){ // Tracklet Analysis // prepare Dut Hit List Int_t iChIdDut = CbmTofAddress::GetUniqueAddress(0,0,0,0,fiDut); fChannelInfo = fDigiPar->GetCell( iChIdDut ); /*TGeoNode *fNode=*/ // prepare global->local trafo gGeoManager->FindNode(fChannelInfo->GetX(),fChannelInfo->GetY(),fChannelInfo->GetZ()); /*TGeoNode* cNode=*/ gGeoManager->GetCurrentNode(); Double_t dDutzPos=fChannelInfo->GetZ(); Double_t hitpos[3], hitpos_local[3]; vector vDutHit; for( Int_t iHitInd = 0; iHitInd < iNbTofHits; iHitInd++) { pHit = (CbmTofHit*) fTofHitsColl->At( iHitInd ); if(NULL == pHit) continue; Int_t iDetId = (pHit->GetAddress() & DetMask); Int_t iSmType=CbmTofAddress::GetSmType( iDetId ); if(iSmType == fiDut) { vDutHit.push_back(pHit); } } LOG(DEBUG)< > vTrkMap; //contains the tracks for a given hit std::vector > vHitMap; //contains the hits for a given track vTrkMap.resize(vDutHit.size()); vHitMap.resize(iNbTofTracks); for (Int_t iTrk=0; iTrkAt(iTrk); if(NULL == pTrk) continue; if (pTrk->GetNofHits() < NStations) continue; // Calculate positions and time in Dut plane Double_t dXex=pTrk->GetFitX(dDutzPos); Double_t dYex=pTrk->GetFitY(dDutzPos); Double_t dR=TMath::Sqrt(dXex*dXex + dYex*dYex + dDutzPos*dDutzPos); Double_t dTex=pTrk->GetFitT(dR); for (UInt_t i=0; iGetTime())/fSIGT,2) +TMath::Power(TMath::Abs(dXex-vDutHit[i]->GetX())/fSIGX,2) +TMath::Power(TMath::Abs(dYex-vDutHit[i]->GetY())/fSIGY,2))/3; LOG(DEBUG1)<GetTime()) << FairLogger::endl; if(dChi < fSIGLIM) { // acceptable match if(vHitMap[iTrk].size()>0) { Int_t iCnt=0; for ( std::map::iterator it=vHitMap[iTrk].begin(); it!=vHitMap[iTrk].end(); it++){ iCnt++; // LOG(DEBUG)< %6.2f ?",iTrk,iCnt, it->second,it->first,dChi) // << FairLogger::endl; if(it->first > dChi) { vHitMap[iTrk].insert(--it,std::pair(dChi,i)); //LOG(DEBUG)<second,it->first)<< FairLogger::endl; break; } } } else{ vHitMap[iTrk].insert(std::pair(dChi,i)); //LOG(DEBUG)<0) { for ( std::map::iterator it=vTrkMap[i].begin(); it!=vTrkMap[i].end(); it++){ if(it->first > dChi) { vTrkMap[i].insert(--it,std::pair(dChi,iTrk)); break; } } } else{ vTrkMap[i].insert(std::pair(dChi,iTrk)); } } // end of Chi condition if(vTrkMap[i].size()>0) LOG(DEBUG1)<second,vTrkMap[i].begin()->first) <0) LOG(DEBUG)<second,vHitMap[iTrk].begin()->first) < 0) for(Int_t iHit=0; static_cast(iHit)0){ Int_t iTrk=vTrkMap[iHit].begin()->second; // hit was assigned best to track iTrk if(vHitMap[iTrk].begin()->second == iHit) { // unique/consistent assignment LOG(DEBUG)< HitMap[%d]: uni %d, %6.4f ",iHit,iTrk,vHitMap[iTrk].begin()->second,vHitMap[iTrk].begin()->first) <::iterator it=vTrkMap[iHit].begin()++; it != vTrkMap[iHit].end(); it++){ Int_t iTrk1=it->second; if(iTrk != iTrk1) for ( std::map::iterator it1=vHitMap[iTrk1].begin()++; it1 != vHitMap[iTrk1].end(); it1++){ if(it1->second == iHit) { vHitMap[iTrk1].erase(it1); LOG(DEBUG1)<::iterator it=vHitMap[iTrk].begin()++; it != vHitMap[iTrk].end(); it++){ Int_t iHit1=it->second; if(iHit != iHit1) for ( std::map::iterator it1=vTrkMap[iHit1].begin()++; it1 != vTrkMap[iHit1].end(); it1++){ if(it1->second == iTrk) { vTrkMap[iHit1].erase(it1); LOG(DEBUG1)< HitMap[%d]: mis %d, %6.4f < %6.4f ",iHit,iTrk, vHitMap[iTrk].begin()->second,vHitMap[iTrk].begin()->first,vTrkMap[iHit].begin()->first) < TrkMap.size: %d ",iHit,(int)vTrkMap[iHit].size()) <At(iTrk); if(NULL == pTrk) continue; if (pTrk->GetNofHits() < NStations) continue; // Calculate positions and time in Dut plane /* Double_t dXex=pTrk->GetFitX(dDutzPos);*/ /* Double_t dYex=pTrk->GetFitY(dDutzPos);*/ /* Double_t dTex=pTrk->GetFitT(dDutzPos);*/ hitpos[0]=pTrk->GetFitX(dDutzPos); hitpos[1]=pTrk->GetFitY(dDutzPos); hitpos[2]=dDutzPos; gGeoManager->MasterToLocal(hitpos, hitpos_local); if(vHitMap[iTrk].size()>0){ // matched hit found LOG(DEBUG)<GetNofHits(),(int)vHitMap[iTrk].size(),vHitMap[iTrk].begin()->second,vHitMap[iTrk].begin()->first) <second]; Double_t dDX = pHit->GetX() - pTrk->GetFitX(pHit->GetZ()); // - tPar->GetX() - tPar->GetTx()*dDZ; Double_t dDY = pHit->GetY() - pTrk->GetFitY(pHit->GetZ()); // - tPar->GetTy()*dDZ; Double_t dDT = pHit->GetTime() - pTrk->GetFitT(pHit->GetR()); //pTrk->GetTdif(fStationType[iSt]); Double_t dDTB= pTrk->GetTdif(fiDut, pHit); // ignore pHit in calc of reference fhDutPullX->Fill(dDX); fhDutPullY->Fill(dDY); fhDutPullT->Fill(dDT); fhDutPullTB->Fill(dDTB); fhDutChiFound->Fill(pTrk->GetChiSq()); fhDutChiMatch->Fill(vHitMap[iTrk].begin()->first); fhDutXY_Found->Fill(hitpos_local[0],hitpos_local[1]); fhDutXYDT->Fill(hitpos_local[0],hitpos_local[1],dDTB); }else{ // no match for this track fhDutChiMissed->Fill(pTrk->GetChiSq()); fhDutXY_Missed->Fill(hitpos_local[0],hitpos_local[1]); } } } // #tracklets > 0 end } // TclonesArray existing end return kTRUE; } //FillHistos // ------------------------------------------------------------------ Bool_t CbmTofAnaTestbeam::WriteHistos() { LOG(INFO)<<"CbmTofAnaTestbeam::WriteHistos: ./tofAnaTestBeam.hst.root, mode = " << fiCorMode << FairLogger::endl; // Write histogramms to the file TDirectory * oldir = gDirectory; TFile *fHist = new TFile("./tofAnaTestBeam.hst.root","RECREATE"); fHist->cd(); switch(fiCorMode){ case 1 : { TProfile *htmp=fhDTD4DT04D4best->ProfileX(); TH1D *htmp1D=htmp->ProjectionX(); // htmp1D->Draw(); if(fhDTD4DT04D4Off != NULL){ // fhDTD4DT04D4Off->Draw("same"); // fhDTD4DT04D4Off->Write(); //LOG(INFO)<<"Update hDTD4DT04D4best"<GetNbinsX(); for (Int_t ix=0; ixGetBinContent(ix) + fhDTD4DT04D4Off->GetBinContent(ix); // Double_t dVal=fhDTD4DT04D4Off->GetBinContent(ix); LOG(DEBUG1)<<"Update hDTD4DT04D4best "<GetBinContent(ix)<<" + " << fhDTD4DT04D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl; htmp1D->SetBinContent(ix,dVal); } }else { LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTD4DT04D4Off not found " << FairLogger::endl; } // fhDTD4DT04D4best->Write(); htmp1D->Write(); if(fhDTX4D4Off != NULL) fhDTX4D4Off->Write(); if(fhDTY4D4Off != NULL) fhDTY4D4Off->Write(); if(fhDTTexpD4Off != NULL) fhDTTexpD4Off->Write(); } break; case 2 :{ TProfile *htmpx=fhDTX4D4best->ProfileX(); TH1D *htmpx1D=htmpx->ProjectionX(); if(fhDTX4D4Off != NULL){ Double_t nx=htmpx1D->GetNbinsX(); for (Int_t ix=0; ixGetBinContent(ix) + fhDTX4D4Off->GetBinContent(ix); LOG(DEBUG1)<<"Update hDTX4D4best "<GetBinContent(ix)<<" + " << fhDTX4D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl; htmpx1D->SetBinContent(ix,dVal); } }else { LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTX4D4Off not found " << FairLogger::endl; } htmpx1D->Write(); if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write(); if(fhDTY4D4Off != NULL) fhDTY4D4Off->Write(); if(fhDTTexpD4Off != NULL) fhDTTexpD4Off->Write(); } break; case 3 :{ TProfile *htmpx=fhDTY4D4best->ProfileX(); TH1D *htmpx1D=htmpx->ProjectionX(); if(fhDTY4D4Off != NULL){ Double_t nx=htmpx1D->GetNbinsX(); for (Int_t ix=0; ixGetBinContent(ix) + fhDTY4D4Off->GetBinContent(ix); LOG(DEBUG1)<<"Update hDTY4D4best "<GetBinContent(ix)<<" + " << fhDTY4D4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl; htmpx1D->SetBinContent(ix,dVal); } }else { LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTY4D4Off not found " << FairLogger::endl; } htmpx1D->Write(); if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write(); if(fhDTX4D4Off != NULL) fhDTX4D4Off->Write(); if(fhDTTexpD4Off != NULL) fhDTTexpD4Off->Write(); } break; case 4 :{ TProfile *htmpx=fhTexpDT04D4best->ProfileX(); TH1D *htmpx1D=htmpx->ProjectionX(); if(fhDTTexpD4Off != NULL){ Double_t nx=htmpx1D->GetNbinsX(); for (Int_t ix=0; ixGetBinContent(ix) + fhDTTexpD4Off->GetBinContent(ix); LOG(DEBUG1)<<"Update hDTTexpD4best "<GetBinContent(ix)<<" + " << fhDTTexpD4Off->GetBinContent(ix) << " -> " << dVal << FairLogger::endl; htmpx1D->SetBinContent(ix,dVal); } }else { LOG(WARNING)<<"CbmTofAnaTestbeam::WriteHistos: fhDTTexpD4Off not found " << FairLogger::endl; } htmpx1D->Write(); if(fhDTD4DT04D4Off != NULL) fhDTD4DT04D4Off->Write(); if(fhDTX4D4Off != NULL) fhDTX4D4Off->Write(); if(fhDTY4D4Off != NULL) fhDTY4D4Off->Write(); } break; default: ; } // fHist->Write(); if(0) { fhXX2->Write(); fhYY2->Write(); for (Int_t iDet=0; iDet<2; iDet++) { fhXX02[iDet]->Write(); fhYY02[iDet]->Write(); } fhXX04->Write(); fhYY04->Write(); } gDirectory->cd( oldir->GetPath() ); fHist->Close(); return kTRUE; } Bool_t CbmTofAnaTestbeam::DeleteHistos() { // Test class performance // Mapping return kTRUE; } ClassImp(CbmTofAnaTestbeam);