24 #include "FairEventHeader.h"
25 #include "FairRootManager.h"
26 #include "FairRunAna.h"
27 #include "FairRuntimeDb.h"
28 #include "TClonesArray.h"
29 #include "TStopwatch.h"
40 static DetectorId detIds[
NSubDetectors] = {
kRICH,
kDRC,
kDSK,
kEMC,
kGEM,
kLUMI,
kMDT,
kMVD,
kRPC,
kSTT,
kFTOF,
kTOF,
kFTS,
kHYPG,
kHYP };
91 ,fStoreRooTFile(storeFile)
119 FairRootManager* ioman = FairRootManager::Instance();
122 cout <<
"-E- PndEmcAnalysis::Init: "
123 <<
"RootManager not instantiated!" << endl;
129 cout <<
"-W- PndEmcAnalysis::Init: "
130 <<
"No EmcHit array!" << endl;
135 fWaveformArray = (TClonesArray*) ioman->GetObject(
"EmcSortedWaveform");
141 cout <<
"-W- PndEmcAnalysis::Init: "
142 <<
"No PndEmcWaveform array!" << endl;
146 fDigiArray = (TClonesArray*) ioman->GetObject(
"EmcDigi");
148 cout <<
"-W- PndEmcAnalysis::Init: "
149 <<
"No PndEmcDigi array!" << endl;
154 cout <<
"-W- PndEmcAnalysis::Init: "
155 <<
"No PndEmcSharedDigi array!" << endl;
160 cout <<
"-W- PndEmcAnalysis::Init: "
161 <<
"No MCTrack array!" << endl;
164 fClusterArray = (TClonesArray*) ioman->GetObject(
"EmcCluster");
166 cout <<
"-W- PndEmcAnalysis::Init: "
167 <<
"No PndEmcCluster array!" << endl;
170 fBumpArray = (TClonesArray*) ioman->GetObject(
"EmcBump");
172 cout <<
"-W- PndEmcAnalysis::Init: "
173 <<
"No PndEmcBump array!" << endl;
176 fChargedCand = (TClonesArray*)ioman->GetObject(
"PidChargedCand");
178 cout <<
"-W- PndEmcAnalysis::Init: "
179 <<
"No PidChargedCand array!" << endl;
182 fNeutralCand = (TClonesArray*)ioman->GetObject(
"PidNeutralCand");
184 cout <<
"-W- PndEmcAnalysis::Init: "
185 <<
"No PidNeutralCand array!" << endl;
197 tHit =
new TTree(
"hit",
"hit data");
203 tHit->Branch(
"phi" ,&
fPhi ,
"Phi/D");
210 tWave =
new TTree(
"wave",
"wave data");
226 tDigi =
new TTree(
"digi",
"digi data");
239 tBump =
new TTree(
"bump",
"bump data");
244 tBump->Branch(
"bt2" ,&
fBumpTime2 ,
"WeightedTimeByTimeErrorAndEnergy/D");
257 tTask =
new TTree(
"task",
"task data");
270 tEvtPileup=
new TTree(
"evtpu",
"event pileup data");
275 tTaskPileup=
new TTree(
"taskpu",
"task pileup data");
286 tRec =
new TTree(
"ResA",
"Reconstruction");
290 tRec->Branch(
"p1" ,&
fp4[0] ,
"p1[120]/D");
291 tRec->Branch(
"Mcp1" ,&
fMcp4[0] ,
"Mcp1[120]/D");
293 tRec->Branch(
"Z201" ,&
fZ201[0] ,
"Z201[30]/D");
294 tRec->Branch(
"Z531" ,&
fZ531[0] ,
"Z531[30]/D");
295 tRec->Branch(
"Lat1" ,&
fLat1[0] ,
"Lat1[30]/D");
298 tRec->Branch(
"Theta1",&
fTheta1[0] ,
"Theta1[30]/D");
299 tRec->Branch(
"Phi1" ,&
fPhi1[0] ,
"Phi1[30]/D");
306 tMcTruth =
new TTree(
"Res1A",
"Mc Truth");
321 tCharged=
new TTree(
"char",
"charged cand");
331 tNeutral=
new TTree(
"neut",
"neutral cand");
346 tAna3g =
new TTree(
"Ana",
"Jpsi-->3photon");
349 tAna3g->Branch(
"4p1" ,&
f4p1[0] ,
"Photon1[4]/D");
362 tAna3g->Branch(
"4p2" ,&
f4p2[0] ,
"Photon2[4]/D");
375 tAna3g->Branch(
"4p3" ,&
f4p3[0] ,
"Photon3[4]/D");
424 cout <<
"-I- PndEmcAnalysis: Intialization successfull" << endl;
436 fEvtTime = FairRootManager::Instance()->GetEventTime();
437 cout<<
"=======================PndEmcAnalysis======================="<<endl;
440 for(Int_t iHit =0; iHit <
nHits; ++iHit){
453 cout<<
"Save "<<nHits<<
" PndEmcHit done"<<endl;
457 if(FairRunAna::Instance()->IsTimeStamp()){
458 fWaveformArray = FairRootManager::Instance()->GetData(
"EmcSortedWaveform"
460 , FairRootManager::Instance()->GetEventTime() + time_length);
467 for (Int_t iWaveform=0; iWaveform<nWaveforms; iWaveform++)
471 const std::vector<Int_t>& evtList = theWaveform->
GetEvtList();
472 for(
size_t i=0;
i< evtList.size(); ++
i){
475 if(evtList.size() > 1){
498 cout<<
"Save "<<nWaveforms<<
" PndEmcWaveform done"<<endl;
510 for (Int_t iDigi=0; iDigi<nDigi; iDigi++)
526 for (Int_t iDigi=0; iDigi<nSharedDigi; iDigi++)
532 cout<<
"Save "<<nDigi<<
" PndEmcDigi done"<<endl;
536 std::vector<key> goodHits;
540 for (Int_t iBump=0; iBump<nBump; iBump++)
554 const std::vector<Int_t>& digis = theBump->
DigiList();
561 for(
size_t id=0;
id < digis.size(); ++id){
569 fSeedTime = seedDigi == 0? 0 :seedDigi->GetTimeStamp();
577 cout<<
"Save "<<nBump<<
" PndEmcBump done"<<endl;
580 std::vector<key> goodTrack;
582 sort(goodHits.begin(), goodHits.end(), std::greater<key>());
583 if(goodHits.size() >=1){
593 for(Int_t
i=0;
i < nTrack; ++
i){
602 sort(goodTrack.begin(), goodTrack.end(), std::greater<key>());
603 for(
size_t i=0;
i< goodTrack.size(); ++
i){
613 Int_t NEL = goodHits.size() > NElement ? NElement : goodHits.size();
615 for(
size_t i=NEL;
i<goodHits.size();++
i){
635 TVector3
pos(calibrator1->
Where(theBump));
646 memcpy(&
fp4[locIdx*4], &p4[0], 4*
sizeof(
Double_t));
650 cout<<
"Save "<<NEL<<
" reconstruction data done"<<endl;
656 std::vector< std::map<Int_t, Int_t> > trackMap;
664 for(Int_t itrack = 0; itrack < nTrack; ++itrack){
669 std::map<Int_t, Int_t> newMap;
670 newMap.insert(std::pair<Int_t,Int_t>(itrack, pt1->
GetMotherID()));
671 trackMap.push_back(newMap);
676 fPhi1[itrack] = p4.Phi();
693 cout<<
"MC points #"<<pt1->
GetNPoints(
kEMC)<<
", #"<<trackMap.size()<<endl;
695 for(
size_t im=0; im < trackMap.size(); ++im){
696 std::map<Int_t, Int_t>& theMap = trackMap[im];
697 if(theMap.find(pt1->
GetMotherID()) != theMap.end()){
698 theMap.insert(std::pair<Int_t,Int_t>(itrack, pt1->
GetMotherID()));
709 cout<<
"Save MCTruth data done"<<endl;
719 for(Int_t ic=0;ic < nCand; ++ic){
731 cout<<
"Save "<<nCand<<
" charged candidates done"<<endl;
741 for(Int_t ic=0;ic < nCand; ++ic){
759 cout<<
"Save "<<nCand<<
" neutral candidates done"<<endl;
764 cout<<
"begin write #"<<nCand<<
" candidates."<<endl;
765 std::vector<key> goodCand;
766 for(Int_t ic=0;ic<nCand;++ic){
771 if(goodCand.size() >= 3){
772 sort(goodCand.begin(), goodCand.end(), std::greater<key>());
823 for(
size_t jc=3;jc<goodCand.size();++jc){
862 cout<<
"Save task data done"<<endl;
869 cout <<
"PndEmcAnalysis, Real time " << rtime <<
" s, CPU time " << ctime <<
" s" << endl;
873 cout<<
"================================================================="<<endl;
879 FairRun*
run = FairRun::Instance();
880 if ( ! run ) Fatal(
"SetParContainers",
"No analysis run");
882 FairRuntimeDb* db = run->GetRuntimeDb();
883 if ( ! db ) Fatal(
"SetParContainers",
"No runtime database");
898 std::cout<<
"*********************************************************"<<std::endl;
899 std::cout<<
"PndEmcAnalysis::FinishTask"<<std::endl;
900 std::cout<<
"*********************************************************"<<std::endl;
901 std::cout<<
"Total waveforms# "<<
totNumOfWave<<std::endl;
903 std::cout<<
"Total fpga digis# "<<
totNumOfDigi<<std::endl;
908 std::cout<<
"*********************************************************"<<std::endl;
910 std::map<Int_t, std::pair<Int_t,Int_t> >::iterator it =
evtMap.begin();
911 for(; it !=
evtMap.end(); ++it){
924 TFile* gOldFile = gFile;
928 cout<<
"write "<<
tHit->GetTitle()<<endl;
932 cout<<
"write "<<
tWave->GetTitle()<<endl;
936 cout<<
"write "<<
tDigi->GetTitle()<<endl;
940 cout<<
"write "<<
tBump->GetTitle()<<endl;
944 cout<<
"write "<<
tTask->GetTitle()<<endl;
948 cout<<
"write "<<
tRec->GetTitle()<<endl;
952 cout<<
"write "<<
tMcTruth->GetTitle()<<endl;
962 cout<<
"write "<<
tCharged->GetTitle()<<endl;
966 cout<<
"write "<<
tNeutral->GetTitle()<<endl;
970 cout<<
"write "<<
tAna3g->GetTitle()<<endl;
key(Int_t hit, Double_t e)
Double_t GetEmcClusterE25() const
bool operator<(const key &rhs) const
TClonesArray * fEmcHitsArray
virtual Double_t GetEnergy() const
represents the reconstructed hit of one emc crystal
Float_t GetEmcRawEnergy() const
Short_t GetModule() const
Int_t totNumOfPileupEvents
Int_t GetNPoints(DetectorId detId) const
static const Int_t NElement
TClonesArray * fMcTrackArray
Float_t GetEmcCalEnergy() const
std::map< Int_t, std::pair< Int_t, Int_t > > evtMap
std::set< Int_t > pevtset
Int_t GetEmcNumberOfCrystals() const
TLorentzVector GetLorentzVector() const
TClonesArray * fBumpArray
TVector3 GetMomentum() const
static DetectorId detIds[NSubDetectors]
Double_t GetEnergy() const
virtual Double_t Energy(PndEmcCluster *clust, Int_t pid=22)=0
const std::vector< Int_t > & DigiList() const
Double_t GetEmcClusterE9() const
virtual InitStatus Init()
TClonesArray * fClusterArray
Bool_t IsGeneratorCreated(void) const
TLorentzVector Get4Momentum() const
bool operator>(const key &rhs) const
Short_t GetModule() const
Double_t fEnergyC1[NElement]
Double_t fEnergyC2[NElement]
Int_t GetEmcIndex() const
std::vector< PndEmcPoint * > & GetPointList()
PndEmcAnalysis(Int_t verbose=0, Bool_t storeRootFile=kTRUE)
virtual ~PndEmcAnalysis()
Double_t fMcp4[4 *NElement]
TClonesArray * fNeutralCand
virtual Double_t GetEnergy() const
Int_t fNumberOfGoodPhoton
const std::map< Int_t, Int_t > & LocalMaxMap() const
Double_t GetEmcClusterZ53() const
Double_t GetEmcClusterLat() const
void Print(Int_t iTrack=0) const
static PndEmcAbsClusterCalibrator * MakeEmcClusterCalibrator(Int_t method, Int_t version=1)
TClonesArray * fDigiArray
Int_t GetEmcNumberOfBumps() const
Double_t GetEmcClusterZ20() const
Double_t fTheta1[NElement]
virtual void FinishTask()
Double_t fPositionV[3 *NElement]
const std::map< Int_t, Int_t > & MemberDigiMap() const
virtual TVector3 Where(PndEmcCluster *clust, Int_t pid=22)=0
virtual Double_t GetTime() const
TClonesArray * fSharedDigiArray
Double_t fp4[4 *NElement]
represents the deposited energy of one emc crystal from simulation
Int_t NumberOfDigis() const
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
virtual void Exec(Option_t *opt)
TVector3 GetPosition() const
virtual Double_t energy() const
static const Int_t NSubDetectors
const TVector3 & where() const
TClonesArray * fChargedCand
Double_t fEnergy[NElement]
Bool_t fTimeOrderedDigi
set to kTRUE to use the time ordering of the output data.
Int_t GetMotherID() const
TVector3 GetMomentum() const
TClonesArray * fWaveformArray
Short_t GetModule() const
Int_t GetEmcModule() const
represents a reconstructed (splitted) emc cluster
Double_t totSharedDigiEnergy