diff --git a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
index 420d21f4d755f3be42b1f56b07c8a197bd8640db..b83a210425b65873022631c3b5f4cffe81be97eb 100644
--- a/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
+++ b/analysis/PWGHAD/hadron/CbmHadronAnalysis.cxx
@@ -40,12 +40,19 @@ using namespace std;
#include "CbmHadron.h"
#include "CbmHadronAnalysis.h"
#include "CbmMatch.h"
+#include "CbmEvent.h"
#include "CbmVertex.h"
#include "CbmKFVertex.h"
#include "CbmStsKFTrackFitter.h"
+#include "CbmDigiManager.h"
-TClonesArray* fTofTracks; // CbmTofTrack array
+CbmDigiManager* fDigiMan; // TOF Input Digis
+TClonesArray* fEventsColl; // CBMEvents (time based)
+TClonesArray* fTofTracks; // CbmTofTrack array
+TClonesArray* fTofHitsColl;
+TClonesArray* fStsHitsColl;
+static Int_t iNbTs=0;
#define M2PI 0.019479835
#define M2KA 0.24371698
@@ -55,6 +62,7 @@ TClonesArray* fTofTracks; // CbmTofTrack array
Int_t nMCTracks=0, nTofPoints=0, nTofHits=0;
Int_t nTofTracks=0;
+ Int_t nStsHits=0;
Int_t nGlobTracks=0;
Int_t NMASS=3;
Float_t refMass[3]={0.139, 0.494, 0.938};
@@ -1256,22 +1264,38 @@ InitStatus CbmHadronAnalysis::Init()
// Task initialisation
FairRootManager* rootMgr = FairRootManager::Instance();
if(NULL == rootMgr) {
- cout << "-E- CbmHadronAnalysis::Init : "
- << "ROOT manager is not instantiated." << endl;
- return kFATAL;
+ cout << "-E- CbmHadronAnalysis::Init : ROOT manager is not instantiated." << endl;
+ return kFATAL;
}
+
+ fEventsColl = dynamic_cast<TClonesArray*>(rootMgr->GetObject("Event"));
+ if(NULL == fEventsColl)
+ fEventsColl = dynamic_cast<TClonesArray*>(rootMgr->GetObject("CbmEvent"));
+
+ if(NULL == fEventsColl)
+ LOG(info) << "CbmEvent not found in input file, assume eventwise input" ;
+ else
+ LOG(info) << "CbmEvent found in input file, timebased analysis" ;
+
fMCEventHeader = (FairMCEventHeader*) rootMgr->GetObject("MCEventHeader.");
if(NULL == fMCEventHeader) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no MC Header Info" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no MC Header Info" << endl;
}
fTofPoints = (TClonesArray*) rootMgr->GetObject("TofPoint");
if(NULL == fTofPoints) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no TOF point array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no TOF point array!" << endl;
}
-
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) {
+ LOG(error) << GetName() << ": No input digis!";
+ return kFATAL;
+ }else
+ LOG(info) << "DigiManager has Tof Digis";
+
+ /*
fTofDigis = (TClonesArray *) rootMgr->GetObject("TofDigi");
if( NULL == fTofDigis)
@@ -1279,30 +1303,29 @@ InitStatus CbmHadronAnalysis::Init()
if( NULL == fTofDigis)
fTofDigis = (TClonesArray *) rootMgr->GetObject("CbmTofDigi");
-
if( NULL == fTofDigis)
{
LOG(error)<<"CbmHadronAnalysis::RegisterInputs => Could not get the TofDigi TClonesArray!!!";
} // if( NULL == fTofDigis)
- else {
- fTofDigiMatchColl = (TClonesArray*) rootMgr->GetObject("TofDigiMatch");
- if(NULL == fTofDigiMatchColl) {
- cout << "-I- CbmHadronAnalysis::Init : "
- << "no TOF digiMatch array!" << endl;
- }
-
- fTofDigiMatchPointsColl = (TClonesArray*) rootMgr->GetObject("TofDigiMatchPoints");
- if(NULL == fTofDigiMatchPointsColl) {
- cout << "-I- CbmHadronAnalysis::Init : "
- << "no TOF digiMatchPoints array!" << endl;
+*/
+
+ fTofDigiMatchColl = (TClonesArray*) rootMgr->GetObject("TofCalDigiMatch");
+ if(NULL == fTofDigiMatchColl) {
+ cout << "-I- CbmHadronAnalysis::Init : no TOF CalDigiMatch array!" << endl;
+ fTofDigiMatchColl = (TClonesArray*) rootMgr->GetObject("TofDigiMatch");
+ if(NULL == fTofDigiMatchColl) {
+ cout << "-I- CbmHadronAnalysis::Init : no TOF digiMatch array!" << endl;
+ }
+ }
- }
+ fTofDigiMatchPointsColl = (TClonesArray*) rootMgr->GetObject("TofDigiMatchPoints");
+ if(NULL == fTofDigiMatchPointsColl) {
+ cout << "-I- CbmHadronAnalysis::Init : no TOF digiMatchPoints array!" << endl;
}
- fTofHits = (TClonesArray*) rootMgr->GetObject("TofHit");
- if(NULL == fTofHits) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no TOF Hit array!" << endl;
+ fTofHitsColl = (TClonesArray*) rootMgr->GetObject("TofHit");
+ if(NULL == fTofHitsColl) {
+ cout << "-W- CbmHadronAnalysis::Init : no TOF Hit array!" << endl;
}
fTofTracks = (TClonesArray*) rootMgr->GetObject("TofTrack");
@@ -1313,47 +1336,39 @@ InitStatus CbmHadronAnalysis::Init()
fTrdPoints = (TClonesArray*) rootMgr->GetObject("TrdPoint");
if(NULL == fTofPoints) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no TRD point array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no TRD point array!" << endl;
}
fTrdHits = (TClonesArray*) rootMgr->GetObject("TrdHit");
- if(NULL == fTofHits) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no TRD Hit array!" << endl;
+ if(NULL == fTrdHits) {
+ cout << "-W- CbmHadronAnalysis::Init : no TRD Hit array!" << endl;
}
fGlobalTracks = (TClonesArray*) rootMgr->GetObject("GlobalTrack");
if(NULL == fGlobalTracks) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no Global Track array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no Global Track array!" << endl;
}
fHadrons = (TClonesArray*) rootMgr->GetObject("Hadron");
if(NULL == fHadrons) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no Hadron array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no Hadron array!" << endl;
}
fStsTracks = (TClonesArray*) rootMgr->GetObject("StsTrack");
if(NULL == fStsTracks) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no STS Track array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no STS Track array!" << endl;
}
- fStsHits = (TClonesArray*) rootMgr->GetObject("StsHit");
- if(NULL == fStsHits) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no STS Hit array!" << endl;
+ fStsHitsColl = (TClonesArray*) rootMgr->GetObject("StsHit");
+ if(NULL == fStsHitsColl) {
+ cout << "-W- CbmHadronAnalysis::Init : no STS Hit array!" << endl;
}
fStsClusters = (TClonesArray*) rootMgr->GetObject("StsCluster");
if(NULL == fStsClusters) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no STS Cluster array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no STS Cluster array!" << endl;
}
fStsDigis = (TClonesArray*) rootMgr->GetObject("StsDigi");
if(NULL == fStsDigis) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no STS Digi array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : no STS Digi array!" << endl;
}
fStsDigiMatchColl = (TClonesArray*) rootMgr->GetObject("StsDigiMatch");
if(NULL == fStsDigiMatchColl) {
@@ -1372,7 +1387,14 @@ InitStatus CbmHadronAnalysis::Init()
LOG(warn) << "No primary vertex";
}
-
+ if(NULL == fEventsColl) {
+ fTofHits=fTofHitsColl;
+ fStsHits=fStsHitsColl;
+ } else { // for time based analysis generate temporary Event TClonesArrays
+ fTofHits = new TClonesArray("CbmTofHit",100);
+ fStsHits = new TClonesArray("CbmStsHit",100);
+ }
+
// --- MC data manager
CbmMCDataManager* mcManager =
(CbmMCDataManager*)rootMgr->GetObject("MCDataManager");
@@ -1384,15 +1406,13 @@ InitStatus CbmHadronAnalysis::Init()
}
fMCTracksColl = (TClonesArray*) rootMgr->GetObject("MCTrack");
- if(NULL == fMCTracks) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no MC track array!" << endl;
+ if(NULL == fMCTracksColl) {
+ cout << "-W- CbmHadronAnalysis::Init : " << "no MC track array!" << endl;
}
fStsPointsColl = (TClonesArray*) rootMgr->GetObject("StsPoint");
if(NULL == fStsPointsColl) {
- cout << "-W- CbmHadronAnalysis::Init : "
- << "no STS Point array!" << endl;
+ cout << "-W- CbmHadronAnalysis::Init : " << "no STS Point array!" << endl;
}
fTrackFitter.Init();
@@ -1401,6 +1421,7 @@ InitStatus CbmHadronAnalysis::Init()
if(kSUCCESS != status) {
return status;
}
+
InitStatus statf = ReadFlowFile();
if(kSUCCESS != statf) {
return statf;
@@ -1422,11 +1443,57 @@ InitStatus CbmHadronAnalysis::Init()
return kSUCCESS;
}
// ------------------------------------------------------------------
-
-
+// ------------------------------------------------------------------
+void CbmHadronAnalysis::Exec(Option_t* option)
+{
+ // Task execution on TS or Events
+ if(NULL != fEventsColl) { // working for STS and TOF only
+ iNbTs++;
+ LOG(debug) << "process TS "<<iNbTs<<" with " << fEventsColl->GetEntriesFast() << " events";
+ if( NULL == fTofHits || NULL == fTofHitsColl ) {
+ assert("Invalid pointer to Tof TClonesArray");
+ }
+ if( NULL == fStsHits || NULL == fStsHitsColl ) {
+ assert("Invalid pointer to Sts TClonesArray");
+ }
+ LOG(debug) << "TS contains " << fTofHitsColl->GetEntriesFast() << "TOF and "
+ << fStsHitsColl->GetEntriesFast() << " Sts hits";
+
+ for(Int_t iEvent = 0; iEvent < fEventsColl->GetEntriesFast(); iEvent++)
+ {
+ CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
+ // copy TOF hits
+ fTofHits->Clear();
+ Int_t iNbHits=0;
+ LOG(debug) << "Fill Tof array with mul " << tEvent->GetNofData(ECbmDataType::kTofHit);
+ for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++) {
+ Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
+ CbmTofHit* pHit = (CbmTofHit *)fTofHitsColl->At(iHitIndex);
+ new((*fTofHits)[iNbHits]) CbmTofHit(*pHit); // fill temporary working TClonesArray
+ iNbHits++;
+ }
+ // copy STS hits
+ iNbHits=0;
+ fStsHits->Clear();
+ LOG(debug) << "Fill Sts array with mul " << tEvent->GetNofData(ECbmDataType::kStsHit);
+ for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kStsHit); iHit++) {
+ Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kStsHit, iHit));
+ const CbmStsHit* pHit = (CbmStsHit *)fStsHitsColl->At(iHitIndex);
+ new((*fStsHits)[iNbHits]) CbmStsHit(*pHit); // fill temporary working TClonesArray
+ iNbHits++;
+ }
+ LOG(debug) << Form("analyze TS %d, Ev %d with %d STS, %d TOF hits",iNbTs, iEvent,
+ fStsHits->GetEntriesFast(), fTofHits->GetEntriesFast());
+ ExecEvent(option);
+ }
+ }
+ else { // event based analysis
+ ExecEvent(option);
+ }
+}
// ------------------------------------------------------------------
-void CbmHadronAnalysis::Exec(Option_t*)
+void CbmHadronAnalysis::ExecEvent(Option_t*)
{
// Task execution
@@ -1464,12 +1531,12 @@ void CbmHadronAnalysis::Exec(Option_t*)
}
Float_t yrp_mid=GetMidY(); // midrapidity -> update from simulation!
-
- nMCTracks = fMCTracksColl->GetEntriesFast();
- nTofPoints = fTofPoints->GetEntriesFast();
- nTofHits = fTofHits->GetEntriesFast();
+ if(fMCTracksColl != NULL) nMCTracks = fMCTracksColl->GetEntriesFast();
+ if(fTofPoints != NULL) nTofPoints = fTofPoints->GetEntriesFast();
+ if(fTofHits != NULL) nTofHits = fTofHits->GetEntriesFast();
if(fTofTracks != NULL) nTofTracks = fTofTracks->GetEntriesFast();
if(fGlobalTracks != NULL) nGlobTracks = fGlobalTracks->GetEntriesFast();
+ if(fStsHits != NULL) nStsHits = fStsHits->GetEntriesFast();
if(verbose>0){ //nh-debug
@@ -1478,28 +1545,42 @@ void CbmHadronAnalysis::Exec(Option_t*)
<< ", TofHit " << nTofHits
<< ", TofTrk " << nTofTracks
<< ", GlbTrk " << nGlobTracks
- << ", StsHit " << fStsHits->GetEntriesFast();
- LOG(debug) << "-D- b = "<<fMCEventHeader->GetB()<< ", phi = " << fMCEventHeader->GetRotZ();
+ << ", StsHit " << nStsHits;
+ if(fMCEventHeader != NULL)
+ LOG(debug) << "-D- b = "<<fMCEventHeader->GetB()<< ", phi = " << fMCEventHeader->GetRotZ();
}
if(FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug)) {
for (Int_t j =0; j<nTofHits; j++) {
- TofHit = (CbmTofHit*) fTofHits->At(j);
- if(NULL == TofHit) continue;
- LOG(DEBUG)<<Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ",
- j, TofHit->GetAddress(),TofHit->GetX(),TofHit->GetY(),TofHit->GetZ(),TofHit->GetTime());
-
+ TofHit = (CbmTofHit*) fTofHits->At(j);
+ if(NULL == TofHit) continue;
+ LOG(DEBUG)<<Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ",
+ j, TofHit->GetAddress(),TofHit->GetX(),TofHit->GetY(),TofHit->GetZ(),TofHit->GetTime());
}
}
- if(FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug)) {
+
+ Double_t dT0=0.;
+ for (Int_t j =0; j<nTofHits; j++) {
+ TofHit = (CbmTofHit*) fTofHits->At(j);
+ if(NULL == TofHit) continue;
+ if(TofHit->GetZ()==0.) dT0=TofHit->GetTime();
+ }
+ if( dT0 != 0.) {
+ for (Int_t j =0; j<nTofHits; j++) {
+ TofHit = (CbmTofHit*) fTofHits->At(j);
+ if(NULL == TofHit) continue;
+ TofHit->SetTime( TofHit->GetTime() - dT0);
+ }
+ }
+
+ if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug)) {
for (Int_t j =0; j<nTofHits; j++) {
- TofHit = (CbmTofHit*) fTofHits->At(j);
- if(NULL == TofHit) continue;
- LOG(DEBUG)<<Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ",
- j, TofHit->GetAddress(),TofHit->GetX(),TofHit->GetY(),TofHit->GetZ(),TofHit->GetTime());
-
+ TofHit = (CbmTofHit*) fTofHits->At(j);
+ if(NULL == TofHit) continue;
+ LOG(DEBUG)<<Form("TofHit %d, addr 0x%08x, x %6.1f, y %6.1f, z %6.1f, t %6.1f ",
+ j, TofHit->GetAddress(),TofHit->GetX(),TofHit->GetY(),TofHit->GetZ(),TofHit->GetTime());
}
}
-
+
if (bRecSec && nTofHits > 1) ReconstructSecondaries(); // independent method
// some local arrays
@@ -1512,7 +1593,7 @@ void CbmHadronAnalysis::Exec(Option_t*)
Int_t IndTofTrack_TofHit[nTofHits][MAXNHT]; // index of TofTracks assigned to specific TofHit
// generator level
- fa_mul_b_gen->Fill(fMCEventHeader->GetB(),nMCTracks);
+ if( NULL != fMCEventHeader ) fa_mul_b_gen->Fill(fMCEventHeader->GetB(),nMCTracks);
Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0;
for (Int_t k=0; k< nMCTracks; k++) { // inspect MCTracks
@@ -1782,7 +1863,7 @@ void CbmHadronAnalysis::Exec(Option_t*)
for (Int_t k=0; k< nMCTracks; k++) TrackP[k]=0; // reset track detected flag
Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0;
- fa_mul_b_poi->Fill(fMCEventHeader->GetB(),nTofPoints);
+ if(NULL != fMCEventHeader) fa_mul_b_poi->Fill(fMCEventHeader->GetB(),nTofPoints);
for (Int_t l =0; l<nTofPoints; l++) {
TofPoint = (CbmTofPoint*) fTofPoints->At(l);
@@ -2030,13 +2111,14 @@ void CbmHadronAnalysis::Exec(Option_t*)
while(phirp<-180.) {phirp+=360.;}
while(phirp>180.) {phirp-=360.;}
} // RP flattening end
- delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
- while(delrp<-180.) {delrp+=360.;}
- while(delrp> 180.) {delrp-=360.;}
- fa_phirp_poi->Fill(phirp); // 1D histo
- fa_delrp_b_poi->Fill(fMCEventHeader->GetB(),delrp);
- fa_cdelrp_b_poi->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
-
+ if (NULL != fMCEventHeader) {
+ delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
+ while(delrp<-180.) {delrp+=360.;}
+ while(delrp> 180.) {delrp-=360.;}
+ fa_phirp_poi->Fill(phirp); // 1D histo
+ fa_delrp_b_poi->Fill(fMCEventHeader->GetB(),delrp);
+ fa_cdelrp_b_poi->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
+ }
// TofHit level
for (Int_t k=0; k< nMCTracks; k++) TrackP[k]=0; // reset track detected flag
@@ -2044,7 +2126,7 @@ void CbmHadronAnalysis::Exec(Option_t*)
Int_t NT0=0; Float_t t0m_hit=0.;
Int_t NT0F=0; Float_t t0mf_hit=0.;
Int_t NT0NF=0; Float_t t0mnf_hit=0.;
- fa_mul_b_hit->Fill(fMCEventHeader->GetB(),nTofHits);
+ if (NULL != fMCEventHeader) fa_mul_b_hit->Fill(fMCEventHeader->GetB(),nTofHits);
Float_t T0MIN=0.;
Int_t NFHITS=10;
@@ -2139,7 +2221,7 @@ void CbmHadronAnalysis::Exec(Option_t*)
if(NULL == TofHit) continue;
// Int_t l = TofHit->GetRefId(); // pointer to Digi
Int_t l = j; // One CbmMatch per hit and in same order!!!!
- if(fTofDigiMatchColl != NULL) {
+ if(fTofDigiMatchColl != NULL && fTofDigiMatchPointsColl != NULL ) {
CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchColl->At(l);
// take first digi's point link
CbmLink L0 = digiMatch->GetLink(0);
@@ -2165,9 +2247,9 @@ void CbmHadronAnalysis::Exec(Option_t*)
}
*/
}else{
- lp=-1;
- LOG(WARNING)<<"No Link to MCTofPoint found for hit "<<j;
- continue;
+ lp=-1;
+ LOG(WARNING)<<"No Link to MCTofPoint found for hit "<<j;
+ continue;
}
TofPoint = (CbmTofPoint*) fTofPoints->At(lp);
Int_t k = TofPoint->GetTrackID();
@@ -2455,14 +2537,16 @@ void CbmHadronAnalysis::Exec(Option_t*)
fa_phirps_hit->Fill(phirp1*RADDEG); // 1D histo
fa_phirps_hit->Fill(phirp2*RADDEG); // 1D histo
delrp=(phirp1-phirp2);
- if(0){ //nh-debug
- cout << "<D-hit> Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
+ if (NULL != fMCEventHeader) {
+ if(0){ //nh-debug
+ cout << "<D-hit> Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
+ }
+ fa_cdrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
+ delrp=delrp*RADDEG;
+ if(delrp<-180.) delrp+=360.;
+ if(delrp>180.) delrp-=360.;
+ fa_drp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
}
- fa_cdrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
- delrp=delrp*RADDEG;
- if(delrp<-180.) delrp+=360.;
- if(delrp>180.) delrp-=360.;
- fa_drp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
}
phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2); // full reaction plane
while(phirp<-180.) {phirp+=360.;}
@@ -2483,32 +2567,34 @@ void CbmHadronAnalysis::Exec(Option_t*)
while(phirp<-180.) {phirp+=360.;}
while(phirp>180.) {phirp-=360.;}
} // RP flattening end
- delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
- while(delrp<-180.) {delrp+=360.;}
- while(delrp> 180.) {delrp-=360.;}
- fa_phirp_hit->Fill(phirp); // 1D histo
- fa_delrp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
- fa_cdelrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
-
- fa_tn_hit->Fill(T0MIN);
- fa_t0mn_hit->Fill((Float_t)NT0);
- fa_t0mn_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0);
- if (NT0 > 0) {
- fa_t0m_hit->Fill(t0m_hit);
- fa_t0m_b_hit->Fill(fMCEventHeader->GetB(),t0m_hit);
- }
+ if (NULL != fMCEventHeader) {
+ delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
+ while(delrp<-180.) {delrp+=360.;}
+ while(delrp> 180.) {delrp-=360.;}
+ fa_phirp_hit->Fill(phirp); // 1D histo
+ fa_delrp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
+ fa_cdelrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
+
+ fa_tn_hit->Fill(T0MIN);
+ fa_t0mn_hit->Fill((Float_t)NT0);
+ fa_t0mn_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0);
+ if (NT0 > 0) {
+ fa_t0m_hit->Fill(t0m_hit);
+ fa_t0m_b_hit->Fill(fMCEventHeader->GetB(),t0m_hit);
+ }
- fa_t0mn_f_hit->Fill((Float_t)NT0F);
- fa_t0mn_f_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0F);
- if (NT0F > 0) {
- fa_t0m_f_hit->Fill(t0mf_hit);
- fa_t0m_f_b_hit->Fill(fMCEventHeader->GetB(),t0mf_hit);
- }
- fa_t0mn_nf_hit->Fill((Float_t)NT0NF);
- fa_t0mn_nf_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0NF);
- if (NT0NF > 0) {
- fa_t0m_nf_hit->Fill(t0mnf_hit);
- fa_t0m_nf_b_hit->Fill(fMCEventHeader->GetB(),t0mnf_hit);
+ fa_t0mn_f_hit->Fill((Float_t)NT0F);
+ fa_t0mn_f_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0F);
+ if (NT0F > 0) {
+ fa_t0m_f_hit->Fill(t0mf_hit);
+ fa_t0m_f_b_hit->Fill(fMCEventHeader->GetB(),t0mf_hit);
+ }
+ fa_t0mn_nf_hit->Fill((Float_t)NT0NF);
+ fa_t0mn_nf_b_hit->Fill(fMCEventHeader->GetB(),(Float_t)NT0NF);
+ if (NT0NF > 0) {
+ fa_t0m_nf_hit->Fill(t0mnf_hit);
+ fa_t0m_nf_b_hit->Fill(fMCEventHeader->GetB(),t0mnf_hit);
+ }
}
// cout << "<I> CbmHadronAnalysis: Number of T0 particles: NTO = " << NT0 << endl;
////////////////////////////////////////////////////////////////////////////////////////////
@@ -2520,7 +2606,7 @@ void CbmHadronAnalysis::Exec(Option_t*)
Int_t NGTofTrack=0;
Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0;
- fa_mul_b_glo->Fill(fMCEventHeader->GetB(),nGlobTracks);
+ if (NULL != fMCEventHeader) fa_mul_b_glo->Fill(fMCEventHeader->GetB(),nGlobTracks);
Int_t NReas=0; //100; // activate reassignment of hits to global tracks
Int_t NRIt=0;
@@ -3485,13 +3571,14 @@ void CbmHadronAnalysis::Exec(Option_t*)
delrp=(phirp1-phirp2);
fa_phirps_glo->Fill(phirp1*RADDEG); // 1D histo
fa_phirps_glo->Fill(phirp2*RADDEG); // 1D histo
+ if (NULL != fMCEventHeader) {
// cout << " Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
- fa_cdrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
- delrp=delrp*RADDEG;
- if(delrp<-180.) delrp+=360.;
- if(delrp>180.) delrp-=360.;
- fa_drp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
-
+ fa_cdrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
+ delrp=delrp*RADDEG;
+ if(delrp<-180.) delrp+=360.;
+ if(delrp>180.) delrp-=360.;
+ fa_drp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
+ }
phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2); // full reaction plane
while(phirp<-180.) {phirp+=360.;}
while(phirp>180.) {phirp-=360.;}
@@ -3510,15 +3597,17 @@ void CbmHadronAnalysis::Exec(Option_t*)
phirp+=dphir*RADDEG;
while(phirp<-180.) {phirp+=360.;}
while(phirp>180.) {phirp-=360.;}
- } // RP flattening end
- delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
- while(delrp<-180.) {delrp+=360.;}
- while(delrp> 180.) {delrp-=360.;}
- fa_phirp_glo->Fill(phirp); // 1D histo
- fa_delrp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
- fa_cdelrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
-
- fa_mul_b_had->Fill(fMCEventHeader->GetB(),NGTofTrack);
+ } // RP flattening end
+ if (NULL != fMCEventHeader) {
+ delrp=phirp - RADDEG*fMCEventHeader->GetRotZ();
+ while(delrp<-180.) {delrp+=360.;}
+ while(delrp> 180.) {delrp-=360.;}
+ fa_phirp_glo->Fill(phirp); // 1D histo
+ fa_delrp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
+ fa_cdelrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
+
+ fa_mul_b_had->Fill(fMCEventHeader->GetB(),NGTofTrack);
+ }
}
// Hadron level
@@ -3652,8 +3741,10 @@ void CbmHadronAnalysis::ReconstructSecondaries()
const Double_t dTofSigT=0.08;
const Double_t dChiTofLim=3.;
- Int_t nStsHits=fStsHits->GetEntriesFast();
- Int_t nTrdHits=fTrdHits->GetEntriesFast();
+ nStsHits=0;
+ if (NULL != fStsHits ) nStsHits=fStsHits->GetEntriesFast();
+ Int_t nTrdHits=0;
+ if (NULL != fTrdHits ) nTrdHits=fTrdHits->GetEntriesFast();
LOG(DEBUG)<<"Secondaries from "<<nTofHits<<" TofHits, "<<nStsHits
<<" StsHits and "<<nTrdHits<<" TrdHits in event "<<iCandEv;
@@ -3764,31 +3855,33 @@ void CbmHadronAnalysis::ReconstructSecondaries()
for (Int_t i=0; i<nTofHits; i++) {
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // don't merge with fake beam counter
for (Int_t i2=0; i2<nTofHits; i2++){
if(i2 != i){
- CbmTofHit* pTofHit2 = (CbmTofHit*) fTofHits->At(i2);
- if(NULL == pTofHit2) continue;
- // Project to plane with smallest z coordinate
- if( pTofHit2->GetZ()<pTofHit->GetZ() ){ //invert order
- CbmTofHit* pTofHittmp= pTofHit;
- pTofHit = pTofHit2;
- pTofHit2 = pTofHittmp;
- }
- Double_t dPosZExp = pTofHit->GetZ() / pTofHit2->GetZ();
- Double_t dPosXExp = pTofHit2->GetX() * dPosZExp;
- Double_t dPosYExp = pTofHit2->GetY() * dPosZExp;
+ CbmTofHit* pTofHit2 = (CbmTofHit*) fTofHits->At(i2);
+ if(NULL == pTofHit2) continue;
+ // Project to plane with smallest z coordinate
+ if( pTofHit2->GetZ()<pTofHit->GetZ() ){ //invert order
+ CbmTofHit* pTofHittmp= pTofHit;
+ pTofHit = pTofHit2;
+ pTofHit2 = pTofHittmp;
+ }
+ Double_t dPosZExp = pTofHit->GetZ() / pTofHit2->GetZ();
+ Double_t dPosXExp = pTofHit2->GetX() * dPosZExp;
+ Double_t dPosYExp = pTofHit2->GetY() * dPosZExp;
Double_t dTimeExp = pTofHit2->GetTime() * dPosZExp;
- Double_t dChi2 = TMath::Power((dPosXExp- pTofHit->GetX())/dTofSigX,2)
- + TMath::Power((dPosYExp- pTofHit->GetY())/dTofSigY,2)
- + TMath::Power((dTimeExp- pTofHit->GetTime())/dTofSigT,2);
- Double_t dChi = TMath::Sqrt(dChi2)/3.;
- fhTofChi->Fill(dChi);
- if (dChi < dChiTofLim) { // merge info
+ Double_t dChi2 = TMath::Power((dPosXExp- pTofHit->GetX())/dTofSigX,2)
+ + TMath::Power((dPosYExp- pTofHit->GetY())/dTofSigY,2)
+ + TMath::Power((dTimeExp- pTofHit->GetTime())/dTofSigT,2);
+ Double_t dChi = TMath::Sqrt(dChi2)/3.;
+ fhTofChi->Fill(dChi);
+ if (dChi < dChiTofLim) { // merge info
pTofHit->SetTime((pTofHit->GetTime()+dTimeExp)*0.5); // update time
- fTofHits->Remove(pTofHit2);
- //pTofHit2->Delete(); // remove from TClonesArray
- LOG(DEBUG) << "Tof Hits "<<i<<" and "<<i2<<" merged " ;
- }
+ fTofHits->Remove(pTofHit2);
+ //pTofHit2->Delete(); // remove from TClonesArray
+ LOG(DEBUG) << "Tof Hits "<<i<<" and "<<i2<<" merged " ;
+ LOG(debug) << "Tof " << i << ", xyz " << pTofHit->GetX() << ", " << pTofHit->GetY() << ", " << pTofHit->GetZ();
+ }
}
}
}
@@ -3797,6 +3890,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
for (Int_t i=0; i<nTofHits; i++) {
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // skip fake beam counter
dStsDistMin[i] =1.E3;
dSts2DistMin[i]=1.E3;
for (Int_t l=0; l<NTrdStations; l++) dTrdDistMin[i][l]=1.E3;
@@ -3811,6 +3905,8 @@ void CbmHadronAnalysis::ReconstructSecondaries()
Double_t dDist2 = TMath::Power(pStsHit->GetX()-sPosXext,2) + TMath::Power(pStsHit->GetY()-sPosYext,2);
Double_t dDist = TMath::Sqrt(dDist2);
fhDperp->Fill(dDist);
+ LOG(debug) << "Sts " << j << ", xyz " << pStsHit->GetX() << ", " << pStsHit->GetY() << ", " << sPosZ;
+ LOG(debug) << "Tof " << i << ", xyz " << pTofHit->GetX() << ", " << pTofHit->GetY() << ", " << pTofHit->GetZ();
LOG(debug) << "Tof "<<i<<", Sts "<<j<<Form(" -> dist %6.3f, Min %6.3f ",dDist, dStsDistMin[i]);
if(dDist<fdDistPrimLim && dDist<dStsDistMin[i]) { // primary or proton candidate
@@ -3875,6 +3971,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
if(j<0) continue; // no STS assigned
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // skip fake beam counter
CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
@@ -3892,22 +3989,22 @@ void CbmHadronAnalysis::ReconstructSecondaries()
<< Form("Module 0x%08x, layer %d",pTrdHit->GetAddress(),CbmTrdAddress::GetLayerId(pTrdHit->GetAddress()))
<<" at z= "<< pTrdHit->GetZ()<<" dD = "<<dDtrans<<" < " << fdDistTRD;
fhDTRDprim->Fill(dDtrans);
- if(dDtrans<fdDistTRD && dDtrans<dTrdDistMin[i][iTrdLayer]) { // check if acceptable and take best match
- Int_t iMul=iTRD[i].size();
- if(dTrdDistMin[i][iTrdLayer]<1.E3){ // modify previous entry
- //find old entry in vector
- Int_t ll=0;
- for( ; ll<iMul; ll++)
- if( CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
- iTRD[i][ll]=l;
- }else { //add hit
- dTrdDistMin[i][iTrdLayer]=dDtrans;
- iTRD[i].resize( iMul + 1);
- iTRD[i][iMul]=l;
- }
- LOG(DEBUG)<<"assign TrdHit "<<l<<" to TofHit "<<i<<" in layer "<<iTrdLayer<<" with d = "<<dDtrans
- <<", TrdMul"<<iMul<<", dEdx = "<<pTrdHit->GetELoss();
- }
+ if(dDtrans<fdDistTRD && dDtrans<dTrdDistMin[i][iTrdLayer]) { // check if acceptable and take best match
+ Int_t iMul=iTRD[i].size();
+ if(dTrdDistMin[i][iTrdLayer]<1.E3){ // modify previous entry
+ //find old entry in vector
+ Int_t ll=0;
+ for( ; ll<iMul; ll++)
+ if( (Int_t)CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
+ iTRD[i][ll]=l;
+ }else { //add hit
+ dTrdDistMin[i][iTrdLayer]=dDtrans;
+ iTRD[i].resize( iMul + 1);
+ iTRD[i][iMul]=l;
+ }
+ LOG(DEBUG)<<"assign TrdHit "<<l<<" to TofHit "<<i<<" in layer "<<iTrdLayer<<" with d = "<<dDtrans
+ <<", TrdMul"<<iMul<<", dEdx = "<<pTrdHit->GetELoss();
+ }
}
}
//2.b - monitor TRD dEdx
@@ -3932,6 +4029,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
if(j<0) continue; // no STS assigned
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // skip fake beam counter
CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At(j);
Double_t dDx = pTofHit->GetX() - pStsHit->GetX();
Double_t dDy = pTofHit->GetY() - pStsHit->GetY();
@@ -3955,8 +4053,8 @@ void CbmHadronAnalysis::ReconstructSecondaries()
//find old entry in vector
Int_t ll=0;
for( ; ll<iMul; ll++)
- if( CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
- iTRD[i][ll]=l;
+ if( CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
+ iTRD[i][ll]=l;
}else { //add hit
dTrdDistMin[i][iTrdLayer]=dDtrans;
iTRD[i].resize( iMul + 1);
@@ -3986,6 +4084,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
for (Int_t i=0; i<nTofHits; i++) {
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // skip fake beam counter
if( iStsMin[i][0]>-1 && iStsMin[i][1]>-1 ) {
CbmStsHit* pStsHit = (CbmStsHit*) fStsHits->At( iStsMin[i][0] );
CbmStsHit* pSts2Hit = (CbmStsHit*)fStsHits->At( iStsMin[i][1] );
@@ -4032,6 +4131,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
Int_t kbest=-1;
CbmTofHit* pTofHit = (CbmTofHit*) fTofHits->At(i);
if(NULL == pTofHit) continue;
+ if(pTofHit->GetZ() == 0) continue; // skip fake beam counter
for (Int_t j=0; j<nStsHits; j++) {
LOG(debug) << "Tof "<<i<<", Sts "<<j<<Form(" ? sec cand %6.3f Min %6.3f ", dTofDistMin[j], fdDistPrimLim);
if( dTofDistMin[j] > fdDistPrimLim) { // Sts hit not in the primary class
@@ -4093,7 +4193,7 @@ void CbmHadronAnalysis::ReconstructSecondaries()
//find old entry in vector
Int_t ll=0;
for( ; ll<iMul; ll++)
- if( CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
+ if( (Int_t) CbmTrdAddress::GetLayerId(((CbmTrdHit*)fTrdHits->At(iTRD[i][ll]))->GetAddress()) == iTrdLayer ) break;
iTRD[i][ll]=l;
}else { //add hit
dTrdDistMin[i][iTrdLayer]=dDtrans;
diff --git a/analysis/PWGHAD/hadron/CbmHadronAnalysis.h b/analysis/PWGHAD/hadron/CbmHadronAnalysis.h
index ada9f66f2d2256fe0249a62ac4676a296e031291..b31740eb7a42a8e41450ae2d28f5ee1b037979a3 100644
--- a/analysis/PWGHAD/hadron/CbmHadronAnalysis.h
+++ b/analysis/PWGHAD/hadron/CbmHadronAnalysis.h
@@ -451,6 +451,7 @@ public:
virtual InitStatus Init();
virtual void Exec(Option_t* option);
+ virtual void ExecEvent(Option_t* option);
virtual void Finish();
void WriteHistogramms();
diff --git a/analysis/detectors/tof/CbmTofAnaTestbeam.cxx b/analysis/detectors/tof/CbmTofAnaTestbeam.cxx
index 915190de2851c4cddd2380854d939631860c6765..990df8c9671287b4528ea3118ba30b115e6957b9 100644
--- a/analysis/detectors/tof/CbmTofAnaTestbeam.cxx
+++ b/analysis/detectors/tof/CbmTofAnaTestbeam.cxx
@@ -18,6 +18,7 @@
#include "CbmMatch.h"
#include "CbmTrackMatchNew.h"
#include "CbmTofTracklet.h"
+#include "CbmTofTrackletTools.h"
#include "CbmEvent.h"
#include "CbmVertex.h"
#include "CbmTofPoint.h"
@@ -79,8 +80,7 @@ Double_t dDTD4Min=1.E8;
static Double_t StartAnalysisTime = 0.;
static Double_t dTLEvt=0.;
static Double_t StartSpillTime = -100.;
-//const Double_t SpillDuration = 20.; // in seconds
-const Double_t fdSpillBreak = 0.9;
+
Int_t iNspills=0;
static Double_t fdMemoryTime = 1.E12; // memory time in ns
@@ -405,6 +405,10 @@ CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose)
fhDutMul_Missed(NULL),
fhDutTIS_Found(NULL),
fhDutTIS_Missed(NULL),
+ fhDutTIR_Found(NULL),
+ fhDutTIR_Missed(NULL),
+ fhDutVel_Found(NULL),
+ fhDutVel_Missed(NULL),
fhDutDTLH_CluSize(NULL),
fhDutDTLH_Tot(NULL),
fhDutDTLH_Mul(NULL),
@@ -554,6 +558,7 @@ CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose)
fhSelHitTupleResidualXYT_Width(NULL),
fdMulDMax(0.),
fdSpillDuration(20.),
+ fdSpillBreak(0.9),
fdDTDia(0.),
fdDTD4MAX(0.),
fdMul0Max(0.),
@@ -610,6 +615,7 @@ CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose)
fEnableMatchPosScaling(kTRUE),
fFindTracks(NULL),
fClusterizer(NULL),
+ fTrackletTools(NULL),
fbMonteCarloComparison(kFALSE),
fbPointsInInputFile(kFALSE),
fbTracksInInputFile(kFALSE),
@@ -636,6 +642,7 @@ CbmTofAnaTestbeam::CbmTofAnaTestbeam(const char* name, Int_t verbose)
fbAttachDutHitToTracklet(kFALSE),
fbBestSelTrackletOnly(kFALSE),
fbUseSigCalib(kFALSE),
+ fiAnaMode(0),
fdMCSIGLIM(3.),
fdMCSIGT(100.),
fdMCSIGX(1.),
@@ -682,6 +689,7 @@ InitStatus CbmTofAnaTestbeam::Init()
fFindTracks = CbmTofFindTracks::Instance();
fClusterizer= CbmTofTestBeamClusterizer::Instance();
+ fTrackletTools = new CbmTofTrackletTools(); // initialize tools
if (NULL == fFindTracks) {
//fdTShift += fChannelInfoDut->GetZ()/30.; // in ns
@@ -786,8 +794,7 @@ void CbmTofAnaTestbeam::Exec(Option_t* opt)
{
CbmEvent* tEvent = dynamic_cast<CbmEvent*>(fEventsColl->At(iEvent));
LOG(debug) << "Process event "<<iEvent<<" with "<< tEvent->GetNofData(ECbmDataType::kTofHit)<<" hits from "
- << tEvent->GetNofData(ECbmDataType::kTofDigi) << ", "<<tEvent->GetNofData(ECbmDataType::kTofCalDigi)<< " digis "
- ;
+ << tEvent->GetNofData(ECbmDataType::kTofDigi) << ", "<<tEvent->GetNofData(ECbmDataType::kTofCalDigi)<< " digis ";
if(fTofDigisColl) fTofDigisColl->Clear("C");
//if(fTofDigisColl) fTofDigisColl->Delete();
@@ -795,42 +802,46 @@ void CbmTofAnaTestbeam::Exec(Option_t* opt)
if(fTofDigiMatchColl) fTofDigiMatchColl->Clear("C");
if(fTofTrackColl) fTofTrackColl->Clear("C");
-
Int_t iNbDigi=0;
+ /*
+ assert ( fTofDigisColl );
for (Int_t iDigi = 0; iDigi < tEvent->GetNofData(ECbmDataType::kTofCalDigi); iDigi++)
{
Int_t iDigiIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofCalDigi, iDigi));
- CbmTofDigi* tDigi = dynamic_cast<CbmTofDigi*>(fTofDigisCollIn->At(iDigiIndex));
- new((*fTofDigisColl)[iNbDigi++]) CbmTofDigi(*tDigi);
+ const CbmTofDigi* tDigi = fDigiMan->Get<CbmTofDigi>(iDigiIndex);
+// CbmTofDigi* tDigi = dynamic_cast<CbmTofDigi*>(fTofDigisCollIn->At(iDigiIndex));
+ assert (tDigi);
+ //LOG(INFO) << "Copy TofDigi " << iDigi << " from " << iDigiIndex << " to " << iNbDigi;
+ //new((*fTofDigisColl)[iNbDigi++]) CbmTofDigi(*tDigi); // does not work for tDigi, since no TObject
}
-
+ */
Int_t iNbHits=0;
for (Int_t iHit = 0; iHit < tEvent->GetNofData(ECbmDataType::kTofHit); iHit++)
{
Int_t iHitIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofHit, iHit));
CbmTofHit* tHit = dynamic_cast<CbmTofHit*>(fTofHitsCollIn->At(iHitIndex));
- new((*fTofHitsColl)[iNbHits]) CbmTofHit(*tHit);
+ new((*fTofHitsColl)[iNbHits]) CbmTofHit(*tHit);
- CbmMatch* tMatch = dynamic_cast<CbmMatch*>(fTofDigiMatchCollIn->At(iHitIndex));
- new((*fTofDigiMatchColl)[iNbHits]) CbmMatch(*tMatch);
+ CbmMatch* tMatch = dynamic_cast<CbmMatch*>(fTofDigiMatchCollIn->At(iHitIndex));
+ new((*fTofDigiMatchColl)[iNbHits]) CbmMatch(*tMatch);
- iNbHits++;
+ iNbHits++;
}
Int_t iNbTrks=0;
for (Int_t iTrk = 0; iTrk < tEvent->GetNofData(ECbmDataType::kTofTrack); iTrk++)
{
Int_t iTrkIndex = static_cast<Int_t>(tEvent->GetIndex(ECbmDataType::kTofTrack, iTrk));
- CbmTofTracklet* tTrk = dynamic_cast<CbmTofTracklet*>(fTofTrackCollIn->At(iTrkIndex));
- new((*fTofTrackColl)[iNbTrks++]) CbmTofTracklet(*tTrk);
+ CbmTofTracklet* tTrk = dynamic_cast<CbmTofTracklet*>(fTofTrackCollIn->At(iTrkIndex));
+ new((*fTofTrackColl)[iNbTrks++]) CbmTofTracklet(*tTrk);
}
ExecEvent(opt);
- if(iNbDigi) fTofDigisColl->Delete(); //Clear("C"); // Clear causes memory leak, FIXME
- if(iNbHits) fTofHitsColl->Clear("C");
+ if(iNbDigi) fTofDigisColl->Delete(); //Clear("C"); // Clear causes memory leak, FIXME
+ if(iNbHits) fTofHitsColl->Clear("C");
if(iNbHits) fTofDigiMatchColl->Delete(); //Clear("C");
- if(iNbTrks) fTofTrackColl->Delete(); //Clear("C");
+ if(iNbTrks) fTofTrackColl->Delete(); //Clear("C");
}
}
@@ -857,9 +868,9 @@ void CbmTofAnaTestbeam::ExecEvent(Option_t* /*option*/)
// if( 0 < fEvents )
if( 0 == ( fEvents%100000 ) && 0 < fEvents )
{
- cout << "-I- CbmTofAnaTestbeam::Exec : "
- << "event " << fEvents
- << " in "<<iNspills<< " spills processed." << endl;
+ LOG(info) << "CbmTofAnaTestbeam::Exec : "
+ << "event " << fEvents
+ << " in "<<iNspills<< " spills processed.";
}
fEvents += 1;
}
@@ -963,6 +974,19 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
FairRootManager *fManager = FairRootManager::Instance();
fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("Event"));
+ if(NULL == fEventsColl)
+ fEventsColl = dynamic_cast<TClonesArray*>(fManager->GetObject("CbmEvent"));
+
+ if(NULL == fEventsColl)
+ LOG(info) << "CbmEvent not found in input file, assume eventwise input" ;
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) {
+ LOG(error) << GetName() << ": No digi input!";
+ return kFALSE;
+ }
+
if (!fEventsColl) {
@@ -973,8 +997,7 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
fTofDigisColl = (TClonesArray *) fManager->GetObject("TofCalDigi");
if( NULL == fTofDigisColl){
- LOG(WARNING)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!! ... continuing with incomplete input "
- ;
+ LOG(WARNING)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!! ... continuing with incomplete input ";
// return kFALSE;
}// if( NULL == fTofDigisColl)
@@ -990,6 +1013,7 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
return kFALSE;
} // if( NULL == fTofDigiMatchColl)
+
fTofTrackColl = (TClonesArray *) fManager->GetObject("TofTracks");
if( NULL == fTofTrackColl) {
LOG(info)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofTracklet TClonesArray!!!";
@@ -1000,6 +1024,14 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
if( NULL == fTrbHeader) {
LOG(info)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofTrbHeader Object";
}
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) {
+ LOG(error) << GetName() << ": No input digis!";
+ //return kFALSE;
+ }else
+ LOG(info) << "DigiManager has Tof Digis";
if(fbMonteCarloComparison)
{
@@ -1123,9 +1155,8 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
fTofDigisCollIn = (TClonesArray *) fManager->GetObject("TofDigi");
*/
if( NULL == fTofDigisCollIn){
- LOG(WARNING)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!! ... continuing with incomplete input "
- ;
- return kFALSE;
+ LOG(WARNING)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the TofDigi TClonesArray!!! ... continuing with incomplete input ";
+ //return kFALSE;
}// if( NULL == fTofDigisColl)
fTofHitsCollIn = (TClonesArray *) fManager->GetObject("TofHit");
@@ -1136,8 +1167,10 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
fTofDigiMatchCollIn= (TClonesArray *) fManager->GetObject("TofDigiMatch");
if( NULL == fTofDigiMatchCollIn) {
- LOG(error)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the Match TClonesArray!!!";
- //return kFALSE;
+ fTofDigiMatchCollIn= (TClonesArray *) fManager->GetObject("TofCalDigiMatch");
+ if( NULL == fTofDigiMatchCollIn) {
+ LOG(error)<<"CbmTofAnaTestbeam::RegisterInputs => Could not get the Match TClonesArray!!!";
+ }
} // if( NULL == fTofDigiMatchColl)
fTofTrackCollIn = (TClonesArray *) fManager->GetObject("TofTracks");
@@ -1151,7 +1184,8 @@ Bool_t CbmTofAnaTestbeam::RegisterInputs()
if ( ! fDigiMan->IsPresent(ECbmModuleId::kTof) ) {
LOG(error) << GetName() << ": No input digis!";
return kFALSE;
- }
+ }else
+ LOG(info) << "DigiManager has Tof Digis";
// Create work Arrays
fTofDigisColl = new TClonesArray("CbmTofDigi",100);
@@ -1449,8 +1483,8 @@ Bool_t CbmTofAnaTestbeam::CreateHistos()
gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !
// define histos here
- Double_t TISmax = 10.; //fdSpillDuration;
- Double_t TISnbins = 100.; //50.;
+ Double_t TISmax = 11.; //fdSpillDuration;
+ Double_t TISnbins = 110.; //50.;
fhRate_all = new TH1F("hRate_all", "Event Rate; Rate (Hz); t (sec)", 1000, 0, 100);
fhTriggerPattern = new TH1I("tof_trb_trigger_pattern", "CTS trigger pattern", 16, 0, 16);
fhTriggerType = new TH1I("tof_trb_trigger_types", "CTS trigger types", 16, 0, 16);
@@ -1461,15 +1495,15 @@ Bool_t CbmTofAnaTestbeam::CreateHistos()
fhTIS_sel2 = new TH1F("TIS_sel2", "Time in Spill (sel2); t (sec)", TISnbins, 0, TISmax);
- Double_t TIRmax = 2000.; // Run duration
- Double_t TIRnbins = 200.;
+ Double_t TIRmax = 1800.; // Run duration
+ Double_t TIRnbins = 18000.;
fhTIR_all = new TH1F("TIR_all", "Time in Run (all); t (sec)", TIRnbins, 0, TIRmax);
fhTIR_sel = new TH1F("TIR_sel", "Time in Run (sel); t (sec)", TIRnbins, 0, TIRmax);
fhTIR_sel1 = new TH1F("TIR_sel1", "Time in Run (sel1); t (sec)", TIRnbins, 0, TIRmax);
fhTIR_sel2 = new TH1F("TIR_sel2", "Time in Run (sel2); t (sec)", TIRnbins, 0, TIRmax);
- Double_t TISmax2 = 10.;
- Double_t TISnbins2 = 1.E2;
+ Double_t TISmax2 = 11.;
+ Double_t TISnbins2 = 110.;
fhTIS_Nhit = new TH2F("TIS_Nhit", "Time in Spill (Nhit); t (sec); N_{hit}", TISnbins2, 0, TISmax2,25,0,50);
fhTIS_Ntrk = new TH2F("TIS_Ntrk", "Time in Spill (Ntrk); t (sec); N_{trk}", TISnbins2, 0, TISmax2,10,0,10);
@@ -1912,12 +1946,24 @@ Bool_t CbmTofAnaTestbeam::CreateHistos()
fhDutMul_Missed=new TH1F( Form("hDutMul_Missed_%03d",iDutId),
Form("hDutMul_Missed_%03d; Hit Multiplicity",iDutId),
32, 0., 32.);
- fhDutTIS_Found=new TH1F( Form("hDutTIS_Found_%03d",iDutId),
+ fhDutTIS_Found=new TH2F( Form("hDutTIS_Found_%03d",iDutId),
Form("hDutTIS_Found_%03d; Time in spill (s)",iDutId),
- TISnbins, 0, TISmax);
- fhDutTIS_Missed=new TH1F( Form("hDutTIS_Missed_%03d",iDutId),
+ TISnbins, 0, TISmax,9,0,9);
+ fhDutTIS_Missed=new TH2F( Form("hDutTIS_Missed_%03d",iDutId),
Form("hDutTIS_Missed_%03d; Time in spill (s)",iDutId),
- TISnbins, 0, TISmax);
+ TISnbins, 0, TISmax,9,0,9);
+ fhDutTIR_Found=new TH2F( Form("hDutTIR_Found_%03d",iDutId),
+ Form("hDutTIR_Found_%03d; Time in spill (s)",iDutId),
+ TIRnbins, 0, TIRmax,9,0,9);
+ fhDutTIR_Missed=new TH2F( Form("hDutTIR_Missed_%03d",iDutId),
+ Form("hDutTIR_Missed_%03d; Time in spill (s)",iDutId),
+ TIRnbins, 0, TIRmax,9,0,9);
+ fhDutVel_Found=new TH1F( Form("hDutVel_Found_%03d",iDutId),
+ Form("hDutVel_Found_%03d; velocity (cm/ns)",iDutId),
+ 50, 0., 50.);
+ fhDutVel_Missed=new TH1F( Form("hDutVel_Missed_%03d",iDutId),
+ Form("hDutVel_Missed_%03d; velocity (cm/ns)",iDutId),
+ 50, 0., 50.);
fhDutDTLH_CluSize=new TH2F( Form("hDutDTLH_CluSize_%03d",iDutId),
Form("hDutDTLH_CluSize_%03d; log(#DeltaT); CluSize",iDutId),
50, 0., 12., 10, 1., 11.);
@@ -3215,6 +3261,8 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
Double_t dDTSpill=dTDia-StartSpillTime;
StartSpillTime=dTDia;
*/
+ if(fdSpillBreak == 0.) StartSpillTime=dTAv;
+
Double_t dDTSpill=dTAv-StartSpillTime;
if( fDetIdMap.size() > 3 && dMulD>0 ){ // FIXME - hardwired constants
Double_t dDTLEvt=dTAv-dTLEvt;
@@ -3228,7 +3276,8 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
dDTSpill/1.E9,fEvents,dMulD,(Int_t)fDetIdMap.size(),dDTLEvt/1.E9);
}
}
- fhRate_all->Fill((dTAv-StartAnalysisTime)/1.E9,1./ fhRate_all->GetBinWidth(1));
+ Double_t dTIR=(dTAv-StartAnalysisTime)/1.E9;
+ fhRate_all->Fill(dTIR,1./ fhRate_all->GetBinWidth(1));
if (StartSpillTime<0) {
LOG(debug) << "SpillStartTime not available, abort Anatestbeam " ;
return kFALSE;
@@ -3286,7 +3335,7 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
if(fTrbHeader != NULL) fhTIS_all->Fill(fTrbHeader->GetTimeInSpill());
else fhTIS_all->Fill((dTAv-StartSpillTime)/1.E9);
- fhTIR_all->Fill((dTAv-StartAnalysisTime)/1.E9);
+ fhTIR_all->Fill(dTIR);
LOG(debug)<<Form(" FoundMatches: %d with first chi2s = %12.1f, %12.1f, %12.1f, %12.1f",iNbMatchedHits,
Chi2List[0],Chi2List[1],Chi2List[2],Chi2List[3])
@@ -3313,7 +3362,7 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
if(fTrbHeader != NULL) fhTIS_sel->Fill(fTrbHeader->GetTimeInSpill());
else fhTIS_sel->Fill((dTAv-StartSpillTime)/1.E9);
- fhTIR_sel->Fill((dTAv-StartAnalysisTime)/1.E9);
+ fhTIR_sel->Fill(dTIR);
fhTofD4sel->Fill(pHitRef->GetTime()-dTDia); // general normalisation
fhDTD4sel->Fill(dDTD4Min); // general normalisation
@@ -3410,7 +3459,7 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
if(fTrbHeader != NULL) fhTIS_sel2->Fill(fTrbHeader->GetTimeInSpill());
else fhTIS_sel2->Fill((dTAv-StartSpillTime)/1.E9);
- fhTIR_sel2->Fill((dTAv-StartAnalysisTime)/1.E9);
+ fhTIR_sel2->Fill(dTIR);
if(NULL != fClusterizer)
if(fClusterizer->fdMemoryTime>0) {
@@ -3736,10 +3785,9 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
fhSelMatchEfficiency->Fill(kTRUE, fiNAccRefTracks);
fhSelHitTupleMatchEfficiencyTIS->Fill(kTRUE, (dTAv - StartSpillTime)/1.E9);
-
if(fTrbHeader != NULL) fhTIS_sel1->Fill(fTrbHeader->GetTimeInSpill());
else fhTIS_sel1->Fill((dTAv-StartSpillTime)/1.E9);
- fhTIR_sel1->Fill((dTAv-StartAnalysisTime)/1.E9);
+ fhTIR_sel1->Fill(dTIR);
if(NULL != fClusterizer)
if(fClusterizer->fdMemoryTime>0) {
@@ -4639,7 +4687,11 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
Bool_t bGoodSelTrackletDutMatch(kFALSE);
auto itDutHitMatch = TrackletMatchedDutHitRedChiSq.find(iTrk);
-
+
+ Double_t dVel= 0.;
+ Double_t dTt=fTrackletTools->FitTt(pTrk,fiDutAddr);
+ if(dTt > 0.) dVel=1./dTt;
+
if(NStations == pTrk->GetNofHits())
{
if(fbAttachDutHitToTracklet)
@@ -4825,6 +4877,7 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
Double_t dDYB= pTrk->GetYdif(fiDutAddr, pHit); // ignore pHit in calc of reference
Double_t dDT = pHit->GetTime() - pTrk->GetFitT(pHit->GetZ()); // pTrk->GetTdif(fStationType[iSt]);
Double_t dDTB= pTrk->GetTdif(fiDutAddr, pHit); // ignore pHit in calc of reference
+
fhDutPullX->Fill(dDX);
fhDutPullXB->Fill(dDXB);
fhDutPullY->Fill(dDY);
@@ -4832,6 +4885,33 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
fhDutPullT->Fill(dDT);
fhDutPullTB->Fill(dDTB);
fhDutChi_Found->Fill(pTrk->GetChiSq());
+
+ Int_t iGet4=-1;
+ Double_t dTot=0.;
+ CbmMatch* digiMatch=(CbmMatch *)fTofDigiMatchColl->At(iDutHitIndex);
+ if( NULL != digiMatch )
+ if ( fDigiMan->IsPresent(ECbmModuleId::kTof) )
+ for (Int_t iLink=0; iLink<digiMatch->GetNofLinks(); iLink++){ // loop over digis
+ CbmLink L0 = digiMatch->GetLink(iLink);
+ Int_t iDigInd0=L0.GetIndex();
+ const CbmTofDigi* pDig0 = fDigiMan->Get<CbmTofDigi>(iDigInd0);
+ if(iGet4==-1) {
+ iGet4 = Int_t(pDig0->GetChannel())%8 + 1;
+ dTot=pDig0->GetTot();
+ }
+ else {
+ if(iGet4>0) {
+ if(pDig0->GetTot() > dTot) {
+ iGet4 = Int_t(pDig0->GetChannel())%8 + 1;
+ dTot=pDig0->GetTot();
+ }
+ //if ( Int_t(pDig0->GetChannel())%8 + 1 != iGet4 ) iGet4=0;
+ }
+ }
+ }
+ else LOG(error) << "no Tof Digis";
+ else LOG(error) << "no Tof Digi Match";
+
if(fbAttachDutHitToTracklet)
{
// fhDutChi_Match->Fill(vHitMap[iTrk].begin()->first);
@@ -4840,7 +4920,10 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
fhDutXY_Found->Fill(hitpos_local[0],hitpos_local[1]);
fhDutDTLH_Found->Fill(dDelTLH);
fhDutMul_Found->Fill( dMul4 );
- fhDutTIS_Found->Fill( dTiS );
+ fhDutTIS_Found->Fill( dTiS , iGet4);
+ fhDutTIR_Found->Fill( dTIR , iGet4);
+ fhDutVel_Found->Fill( dVel );
+
fhDutXYDX->Fill(hitpos_local[0],hitpos_local[1],dDX);
fhDutXYDY->Fill(hitpos_local[0],hitpos_local[1],dDY);
fhDutXYDT->Fill(hitpos_local[0],hitpos_local[1],dDTB);
@@ -5021,7 +5104,10 @@ Bool_t CbmTofAnaTestbeam::FillHistos()
fhDutChi_Missed->Fill(pTrk->GetChiSq());
fhDutXY_Missed->Fill(hitpos_local[0],hitpos_local[1]);
fhDutMul_Missed->Fill( dMul4 );
- fhDutTIS_Missed->Fill( dTiS );
+ Int_t iGet4=((Int_t)hitpos_local[0]+16)%8 +1;
+ fhDutTIS_Missed->Fill( dTiS, iGet4);
+ fhDutTIR_Missed->Fill( dTIR, iGet4);
+ fhDutVel_Missed->Fill( dVel );
fhDutDTLH_Missed->Fill( dDelTLH );
fhDutDTLH_Missed_TIS->Fill( dDelTLH, dTiS );
if(NULL != pLHit){
diff --git a/analysis/detectors/tof/CbmTofAnaTestbeam.h b/analysis/detectors/tof/CbmTofAnaTestbeam.h
index 435704c4f3f090ba66c3e5e1bbfe42a85e7a84dc..26ad67fc25cb019bacb1443637855cee9c089370 100644
--- a/analysis/detectors/tof/CbmTofAnaTestbeam.h
+++ b/analysis/detectors/tof/CbmTofAnaTestbeam.h
@@ -113,7 +113,8 @@ class CbmTofAnaTestbeam : public FairTask {
inline void SetPosYS2Sel (Double_t val) { fdPosYS2Sel = val;}
inline void SetPosYS2SelOff (Double_t val) { fdPosYS2SelOff = val;}
inline void SetSel2TOff (Double_t val) { fdSel2TOff = val;}
- inline void SetSpillDuration (Double_t val) { fdSpillDuration = val;}
+ inline void SetSpillDuration (Double_t val) { fdSpillDuration = val;}
+ inline void SetSpillBreak (Double_t val) { fdSpillBreak = val;}
inline void SetMulDMax (Double_t val) { fdMulDMax = val;}
inline void SetDTDia (Double_t val) { fdDTDia = val;}
inline void SetDTD4MAX (Double_t val) { fdDTD4MAX = val;}
@@ -199,6 +200,7 @@ class CbmTofAnaTestbeam : public FairTask {
inline void SetAttachDutHitToTracklet(Bool_t bval) { fbAttachDutHitToTracklet = bval; }
inline void SetBestSelTrackletOnly(Bool_t bval) { fbBestSelTrackletOnly = bval; }
inline void SetUseSigCalib(Bool_t bval) { fbUseSigCalib = bval; }
+ inline void SetAnaMode(Int_t ival) { fiAnaMode = ival; }
inline void SetMCSIGLIM ( Double_t val ) { fdMCSIGLIM = val; }
inline void SetMCSIGT ( Double_t val ) { fdMCSIGT = val; }
@@ -475,8 +477,12 @@ class CbmTofAnaTestbeam : public FairTask {
TH1 * fhDutDTLH_Missed;
TH1 * fhDutMul_Found;
TH1 * fhDutMul_Missed;
- TH1 * fhDutTIS_Found;
- TH1 * fhDutTIS_Missed;
+ TH2 * fhDutTIS_Found;
+ TH2 * fhDutTIS_Missed;
+ TH2 * fhDutTIR_Found;
+ TH2 * fhDutTIR_Missed;
+ TH1 * fhDutVel_Found;
+ TH1 * fhDutVel_Missed;
TH2 * fhDutDTLH_CluSize;
TH2 * fhDutDTLH_Tot;
TH2 * fhDutDTLH_Mul;
@@ -657,7 +663,8 @@ class CbmTofAnaTestbeam : public FairTask {
TH1 *fhSelTypeNNResidualY_Width; // 'calibration' histo
TH1 *fhSelHitTupleResidualXYT_Width;// 'calibration' histo
Double_t fdMulDMax; // max multiplicity in Diamond counter
- Double_t fdSpillDuration; // spill length in sec
+ Double_t fdSpillDuration; // min. spill length in sec
+ Double_t fdSpillBreak; // min. spill break in sec
Double_t fdDTDia; // max time difference between diamonds
Double_t fdDTD4MAX; // max time difference between reference & diamond
Double_t fdMul0Max; // max multiplicity in Dut
@@ -720,6 +727,7 @@ class CbmTofAnaTestbeam : public FairTask {
CbmTofFindTracks* fFindTracks; // Pointer to Task
CbmTofTestBeamClusterizer* fClusterizer; // Pointer to Task
+ CbmTofTrackletTools* fTrackletTools;// Pointer to Class
Bool_t fbMonteCarloComparison;
Bool_t fbPointsInInputFile;
@@ -749,6 +757,7 @@ class CbmTofAnaTestbeam : public FairTask {
Bool_t fbAttachDutHitToTracklet;
Bool_t fbBestSelTrackletOnly;
Bool_t fbUseSigCalib;
+ Int_t fiAnaMode;
Double_t fdMCSIGLIM;
Double_t fdMCSIGT;