FairRoot/PandaRoot
PndEmcAnalysis.cxx
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 // File and Version Information:
3 // $Id: //
4 // Description:
5 // Class PndEmcAnalysis. Module to take the ADC waveforms and produces digi.
6 //
7 // Software developed for the BaBar Detector at the SLAC B-Factory.
8 // Adapted for the PANDA experiment at GSI
9 //
10 // Author List:
11 // Phil Strother Original Author
12 // Dima Melnichuk - adaption for PANDA
13 // Copyright Information:
14 // Copyright (C) 1996 Imperial College
15 //
16 //----------------------------------------------------------------------
17 
18 #include "PndEmcAnalysis.h"
19 #include "PndEmcWaveform.h"
20 #include "PndEmcDigi.h"
21 #include "PndEmcBump.h"
22 #include "PndEmcDigiPar.h"
23 #include "PndEmcRecoPar.h"
24 #include "FairEventHeader.h"
25 #include "FairRootManager.h"
26 #include "FairRunAna.h"
27 #include "FairRuntimeDb.h"
28 #include "TClonesArray.h"
29 #include "TStopwatch.h"
30 #include "PndMCTrack.h"
31 #include <iostream>
34 
35 using std::cout;
36 using std::endl;
37 using std::fstream;
38 
39 static const Int_t NSubDetectors = 15;
41 
42 
43 struct key
44 {
45  key(Int_t hit, Double_t e) : iHit(hit), dE(e){}
46  Int_t iHit;
48  bool operator < (const key& rhs) const {
49  return dE < rhs.dE;
50  }
51  bool operator > (const key& rhs) const {
52  return dE > rhs.dE;
53  }
54 };
55 //struct IndexKey
56 //{
57 // IndexKey(Int_t oIdx, Int_t newIdx) : GenIdx(oIdx), EnergyIdx(newIdx){}
58 // Int_t GenIdx;
59 // Int_t EnergyIdx;
60 // bool operator < (const IndexKey& rhs) const {
61 // return GenIdx < rhs.GenIdx;
62 // }
63 // bool operator > (const IndexKey& rhs) const {
64 // return GenIdx > rhs.GenIdx;
65 // }
66 //};
67 //bool operator< (const key& lhs, const key& rhs) {
68 // return lhs.dE < rhs.dE;
69 //}
70 //bool operator> (const key& lhs, const key& rhs) {
71 // return lhs.dE > rhs.dE;
72 //}
73 //ostream& operator<<(ostream& os, const TVector3& mom) {
74 // os<<"("<<mom.Px()<<", "<<mom.Py()<<", "<<mom.Pz()<<")";
75 // os<<std::endl;
76 // return os;
77 //}
78 
79 
80 
82  fEmcHitsArray(0)
83  ,fWaveformArray(0)
84  ,fDigiArray(0)
85  ,fSharedDigiArray(0)
86  ,fClusterArray(0)
87  ,fBumpArray(0)
88  ,fMcTrackArray(0)
89  ,fEvtHeaderArray(0)
90  ,fVerbose(verbose)
91  ,fStoreRooTFile(storeFile)
92 {
93  fileName="PndEmcAnalysis.root";
94  fTimeOrderedDigi = kFALSE;
95  fSaveWave = kFALSE;
96  fSaveDigi = kFALSE;
97  fSaveBump = kFALSE;
98  fSaveTask = kFALSE;
99  fSaveHits = kFALSE;
100  fSaveMcTruth = kFALSE;
101  fSaveReco = kFALSE;
102  fSavePidCharged = kFALSE;
103  fSavePidNeutral = kFALSE;
104  fSaveJpsiTo3g = kFALSE;
105 }
106 
107 //--------------
108 // Destructor --
109 //--------------
110 //#define TIMEBASEDSIM
112 {
113 }
114 
115 
117 {
118  // Get RootManager
119  FairRootManager* ioman = FairRootManager::Instance();
120  if ( ! ioman )
121  {
122  cout << "-E- PndEmcAnalysis::Init: "
123  << "RootManager not instantiated!" << endl;
124  return kFATAL;
125  }
126 
127  fEmcHitsArray = (TClonesArray*) ioman->GetObject("EmcHit");
128  if ( ! fEmcHitsArray ) {
129  cout << "-W- PndEmcAnalysis::Init: "
130  << "No EmcHit array!" << endl;
131  fSaveHits = kFALSE;
132  }
133  // Get input array
134  if(fTimeOrderedDigi){
135  fWaveformArray = (TClonesArray*) ioman->GetObject("EmcSortedWaveform");
136  }else{
137  fWaveformArray = (TClonesArray*) ioman->GetObject("EmcWaveform");
138  }
139  //fEvtHeaderArray = (TClonesArray*) ioman->GetObject("EventHeader.");
140  if ( ! fWaveformArray ) {
141  cout << "-W- PndEmcAnalysis::Init: "
142  << "No PndEmcWaveform array!" << endl;
143  fSaveWave = kFALSE;
144  }
145 
146  fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigi");
147  if ( ! fDigiArray ) {
148  cout << "-W- PndEmcAnalysis::Init: "
149  << "No PndEmcDigi array!" << endl;
150  fSaveDigi = kFALSE;
151  }
152  fSharedDigiArray = (TClonesArray*) ioman->GetObject("EmcSharedDigi");
153  if ( ! fSharedDigiArray ) {
154  cout << "-W- PndEmcAnalysis::Init: "
155  << "No PndEmcSharedDigi array!" << endl;
156  fSaveDigi = kFALSE;
157  }
158  fMcTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
159  if ( ! fMcTrackArray ) {
160  cout << "-W- PndEmcAnalysis::Init: "
161  << "No MCTrack array!" << endl;
162  fSaveMcTruth = kFALSE;
163  }
164  fClusterArray = (TClonesArray*) ioman->GetObject("EmcCluster");
165  if ( ! fClusterArray ) {
166  cout << "-W- PndEmcAnalysis::Init: "
167  << "No PndEmcCluster array!" << endl;
168  //fSaveBump = kFALSE;
169  }
170  fBumpArray = (TClonesArray*) ioman->GetObject("EmcBump");
171  if ( ! fBumpArray ) {
172  cout << "-W- PndEmcAnalysis::Init: "
173  << "No PndEmcBump array!" << endl;
174  fSaveBump = kFALSE;
175  }
176  fChargedCand = (TClonesArray*)ioman->GetObject("PidChargedCand");
177  if ( ! fChargedCand ) {
178  cout << "-W- PndEmcAnalysis::Init: "
179  << "No PidChargedCand array!" << endl;
180  fSavePidCharged = kFALSE;
181  }
182  fNeutralCand = (TClonesArray*)ioman->GetObject("PidNeutralCand");
183  if ( ! fNeutralCand ) {
184  cout << "-W- PndEmcAnalysis::Init: "
185  << "No PidNeutralCand array!" << endl;
186  fSavePidNeutral = kFALSE;
187  }
188 
189 
190  fFunctor = new StopTime();
192 
193  if(fStoreRooTFile)
194  fRootFile = new TFile(fileName,"RECREATE");
195 
196  if(fSaveHits){
197  tHit = new TTree("hit", "hit data");
198  tHit->Branch("ie" ,&fEventNo ,"EventNO/I");
199  tHit->Branch("ht" ,&fHitTime ,"HitTime/D");
200  tHit->Branch("mod" ,&fModule ,"Module/I");
201  tHit->Branch("he" ,&fHitEnergy ,"HitEnergy/D");
202  tHit->Branch("theta" ,&fTheta ,"Theta/D");
203  tHit->Branch("phi" ,&fPhi ,"Phi/D");
204  tHit->Branch("np" ,&fNPoint ,"NumOfPoint/I");
205  tHit->Branch("xcor" ,&fXcor ,"XCoor/I");
206  tHit->Branch("ycor" ,&fYcor ,"XCoor/I");
207  }
208 
209  if(fSaveWave){
210  tWave = new TTree("wave", "wave data");
211  tWave->Branch("et" ,&fEvtTime ,"EventTime/D");
212  tWave->Branch("mod" ,&fModule ,"Module/I");
213  tWave->Branch("hitid" ,&fHitIndex ,"HitIndex/I");
214  tWave->Branch("wavet" ,&fWaveTime ,"WaveTime/D");
215  tWave->Branch("peak" ,&fWavePeak ,"WavePeak/D");
216  tWave->Branch("wl" ,&fWaveLen ,"WaveLength/I");
217  tWave->Branch("iw" ,&fWaveIdx ,"WaveIndex/I");
218  tWave->Branch("ie" ,&fEventNo ,"EventNO/I");
219  tWave->Branch("baseE" ,&fBaseline ,"BaseLine/D");
220  tWave->Branch("xcor" ,&fXcor ,"XCoor/I");
221  tWave->Branch("ycor" ,&fYcor ,"YCoor/I");
222  tWave->Branch("wpc" ,&fPileupCount ,"Pileup/I");
223  tWave->Branch("he" ,&fHitEnergy ,"HitEnergy/D");
224  }
225  if(fSaveDigi){
226  tDigi = new TTree("digi", "digi data");
227  tDigi->Branch("digit" ,&fDigiTime ,"DigiTime/D");
228  tDigi->Branch("et" ,&fEvtTime ,"EventTime/D");
229  tDigi->Branch("mod" ,&fModule ,"Module/I");
230  tDigi->Branch("wt" ,&fWaveTime ,"WaveTime/D");
231  tDigi->Branch("hite" ,&fHitEnergy ,"HitEnergy/D");
232  tDigi->Branch("e" ,&fDigiEnergy ,"Energy/D");
233  tDigi->Branch("r" ,&fDistance ,"Path/D");
234  tDigi->Branch("iw" ,&fWaveIdx ,"WaveIndex/I");
235  tDigi->Branch("ie" ,&fEventNo ,"EventNO/I");
236  tDigi->Branch("pos" ,&fPosition[0] ,"Position[3]/D");
237  }
238  if(fSaveBump){
239  tBump = new TTree("bump", "bump data");
240  tBump->Branch("et" ,&fEvtTime ,"EventTime/D");
241  tBump->Branch("ie" ,&fEventNo ,"EventNO/I");
242  tBump->Branch("mod" ,&fModule ,"Module/I");
243  tBump->Branch("bt1" ,&fBumpTime1 ,"WeightedTimeByTimeError/D");
244  tBump->Branch("bt2" ,&fBumpTime2 ,"WeightedTimeByTimeErrorAndEnergy/D");
245  tBump->Branch("bt3" ,&fBumpTime3 ,"WeightedTimeByEnergy/D");
246  tBump->Branch("seedt" ,&fSeedTime ,"SeedTime/D");
247  tBump->Branch("e" ,&fBumpEnergy ,"BumpEnergy/D");
248  tBump->Branch("pos" ,&fPosition[0] ,"Position[3]/D");
249  tBump->Branch("seedpos",&fSeedPosition[0] ,"SeedPosition[3]/D");
250  tBump->Branch("ndigi" ,&fNumDigi ,"NumOfDigi/I");
251  tBump->Branch("nbump" ,&fNumBump ,"NumofBump/I");
252  tBump->Branch("ic" ,&fClusterIdx ,"ClusterIndex/I");
253  tBump->Branch("nseed" ,&fNumOfSeed ,"NumOfSeed/I");
254  tBump->Branch("nmap" ,&fDigiMapSize ,"NumDigiMap/I");
255  }
256  if(fSaveTask){
257  tTask = new TTree("task", "task data");
258  tTask->Branch("ie" ,&fEventNo ,"EventNO/I");
259  tTask->Branch("nhits" ,&fNumOfHits ,"NumOfEmcHit/I");
260  tTask->Branch("nwave" ,&fNumWaveform ,"NumOfWaveform/I");
261  tTask->Branch("ndigi" ,&fNumDigi ,"NumOfDigi/I");
262  tTask->Branch("sdigi" ,&fNumSharedDigi ,"NumOfSharedDigi/I");
263  tTask->Branch("nclus" ,&fNumCluster ,"NumOfCluster/I");
264  tTask->Branch("nbump" ,&fNumBump ,"NumOfBump/I");
265  tTask->Branch("ntrack" ,&fNumMcTrack ,"NumOfMcTrack/I");
266  tTask->Branch("totde" ,&totDigiEnergy ,"TotalDigiEnergy/D");
267  tTask->Branch("totse" ,&totSharedDigiEnergy,"TotalSharedDigiEnergy/D");
268  }
269  if(fTimeOrderedDigi){
270  tEvtPileup= new TTree("evtpu","event pileup data");
271  tEvtPileup->Branch("evtNo" ,&fEventNo ,"evtNo/I");
272  tEvtPileup->Branch("totw" ,&fTotWave ,"totw/I");
273  tEvtPileup->Branch("pilew" ,&fPileupWave ,"pilew/I");
274 
275  tTaskPileup= new TTree("taskpu","task pileup data");
276  tTaskPileup->Branch("nw" ,&totNumOfWave ,"TotNumOfWave/I");
277  tTaskPileup->Branch("nsw" ,&totNumOfSubWave ,"TotNumOfSubWave/I");
278  tTaskPileup->Branch("nd" ,&totNumOfDigi ,"TotNumOfDigi/I");
279  tTaskPileup->Branch("nsd" ,&totNumOfSharedDigi ,"TotNumOfSharedDigi/I");
280  tTaskPileup->Branch("nevt" ,&totNumOfEvents ,"TotNumOfEvents/I");
281  tTaskPileup->Branch("npevt" ,&totNumOfPileupEvents,"TotNumOfPileupEvents/I");
282  tTaskPileup->Branch("match" ,&fNumOfMatch ,"NumOfMatch/I");
283  tTaskPileup->Branch("mismatch",&fNumOfDismatch ,"NumOfDismatch/I");
284  }
285  if(fSaveReco){
286  tRec = new TTree("ResA", "Reconstruction");
287  tRec->Branch("E1" ,&fEnergy[0] ,"E1[30]/D");
288  tRec->Branch("E1C" ,&fEnergyC1[0] ,"E1C[30]/D");
289  tRec->Branch("E2C" ,&fEnergyC2[0] ,"E2C[30]/D");
290  tRec->Branch("p1" ,&fp4[0] ,"p1[120]/D");
291  tRec->Branch("Mcp1" ,&fMcp4[0] ,"Mcp1[120]/D");
292  tRec->Branch("pdg" ,&fPDGCode[0] ,"pdg[30]/I");
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");
296  tRec->Branch("pos" ,&fPositionV[0] ,"Position[30]/D");
297  tRec->Branch("Mod1" ,&fModule1[0] ,"Mod1[30]/I");
298  tRec->Branch("Theta1",&fTheta1[0] ,"Theta1[30]/D");
299  tRec->Branch("Phi1" ,&fPhi1[0] ,"Phi1[30]/D");
300  tRec->Branch("ExtE1" ,&fRestEnergy1 ,"ExtE1/D");
301  tRec->Branch("ExtE2" ,&fRestEnergy2 ,"ExtE2/D");
302  tRec->Branch("ExtE3" ,&fRestEnergy3 ,"ExtE3/D");
303  tRec->Branch("NBump" ,&fNumBump ,"NBump/I");
304  }
305  if(fSaveMcTruth){
306  tMcTruth = new TTree("Res1A", "Mc Truth");
307  tMcTruth->Branch("Mcp1" ,&fMcp4[0] ,"Mcp1[120]/D");
308  tMcTruth->Branch("Theta1" ,&fTheta1[0] ,"Theta1[30]/D");
309  tMcTruth->Branch("Phi1" ,&fPhi1[0] ,"Phi1[30]/D");
310  //tMcTruth->Branch("Det1" ,&fDet1[0] ,"Det1[30]/I");
311  tMcTruth->Branch("isHit1" ,&fHitEmc1[0] ,"SizeOfMcTrackArray[30]/I");
312  tMcTruth->Branch("isHit2" ,&fHitEmc2[0] ,"NumHitOfPrimaryTrack[30]/I");
313  tMcTruth->Branch("isHit3" ,&fNPointV[0] ,"TotalNumHit[30]/I");
314  tMcTruth->Branch("isHit4" ,&fNumGT ,"NumOfSecondaryHit/I");
315  //tMcTruth->Branch("X1" ,&fX1 ,"X1/D");
316  //tMcTruth->Branch("Y1" ,&fY1 ,"Y1/D");
317  //tMcTruth->Branch("Z1" ,&fZ1 ,"Z1/D");
318  }
319  if(fSavePidCharged)
320  {
321  tCharged= new TTree("char", "charged cand");
322  tCharged->Branch("char" ,&fchar[0] ,"char[30]/I");
323  tCharged->Branch("p4" ,&fp4[0] ,"p4[120]/D");
324  tCharged->Branch("pos" ,&fpos[0] ,"pos[90]/D");
325  tCharged->Branch("mom" ,&fmom[0] ,"mom[30]/D");
326  tCharged->Branch("mcIdx" ,&fmcidx[0] ,"mctrackid[30]/I");
327 
328  }
329  if(fSavePidNeutral)
330  {
331  tNeutral= new TTree("neut", "neutral cand");
332  tNeutral->Branch("char" ,&fchar[0] ,"char[30]/I");
333  tNeutral->Branch("mcIdx" ,&fmcidx[0] ,"mctrackid[30]/I");
334  tNeutral->Branch("p4" ,&fp4[0] ,"p4[120]/D");
335  tNeutral->Branch("pos" ,&fpos[0] ,"pos[90]/D");
336  tNeutral->Branch("rawe" ,&frawe[0] ,"rawe[30]/D");
337  tNeutral->Branch("cale" ,&fcale[0] ,"cale[30]/D");
338  tNeutral->Branch("ndigi" ,&fndigi[0] ,"ndigi[30]/I");
339  tNeutral->Branch("mod" ,&fmod[0] ,"EmcModule[30]/I");
340  tNeutral->Branch("cluIdx" ,&fcluIdx[0] ,"clusterIdx[30]/I");
341  tNeutral->Branch("e25" ,&fe25[0] ,"e25[30]/D");
342  tNeutral->Branch("e9" ,&fe9[0] ,"e9[30]/D");
343  tNeutral->Branch("nbump" ,&fnbump[0] ,"nbump[30]/I");
344  }
345  if(fSaveJpsiTo3g){
346  tAna3g = new TTree("Ana", "Jpsi-->3photon");
347  tAna3g->Branch("ng" ,&fNumberOfGoodPhoton ,"GoodPhoton/I");
348  tAna3g->Branch("eE" ,&fExternalEnergy ,"ExtralEnergy/D");
349  tAna3g->Branch("4p1" ,&f4p1[0] ,"Photon1[4]/D");
350  tAna3g->Branch("mc4p1" ,&fMc4p1[0] ,"Truth1[4]/D");
351  tAna3g->Branch("pos1" ,&fpos1[0] ,"pos1[3]/D");
352  tAna3g->Branch("e1" ,&fenergy1 ,"e1/D");
353  tAna3g->Branch("rawe1" ,&femcraw1 ,"rawe1/D");
354  tAna3g->Branch("nd1" ,&fndigi1 ,"NumOfDigi1/I");
355  tAna3g->Branch("mcTrk1" ,&fMcTrack1 ,"McTrackID1/I");
356  tAna3g->Branch("mod1" ,&fmod1 ,"mod1/I");
357  tAna3g->Branch("e251" ,&fe251 ,"e251/D");
358  tAna3g->Branch("lat1" ,&flat1 ,"lat1/D");
359  tAna3g->Branch("z201" ,&fz201 ,"z201/D");
360  tAna3g->Branch("z531" ,&fz531 ,"z531/D");
361 
362  tAna3g->Branch("4p2" ,&f4p2[0] ,"Photon2[4]/D");
363  tAna3g->Branch("mc4p2" ,&fMc4p2[0] ,"Truth2[4]/D");
364  tAna3g->Branch("pos2" ,&fpos2[0] ,"pos2[3]/D");
365  tAna3g->Branch("e2" ,&fenergy2 ,"e2/D");
366  tAna3g->Branch("rawe2" ,&femcraw2 ,"rawe2/D");
367  tAna3g->Branch("nd2" ,&fndigi2 ,"NumOfDigi2/I");
368  tAna3g->Branch("mcTrk2" ,&fMcTrack2 ,"McTrackID2/I");
369  tAna3g->Branch("mod2" ,&fmod2 ,"mod2/I");
370  tAna3g->Branch("e252" ,&fe252 ,"e252/D");
371  tAna3g->Branch("lat2" ,&flat2 ,"lat2/D");
372  tAna3g->Branch("z202" ,&fz202 ,"z202/D");
373  tAna3g->Branch("z532" ,&fz532 ,"z532/D");
374 
375  tAna3g->Branch("4p3" ,&f4p3[0] ,"Photon3[4]/D");
376  tAna3g->Branch("mc4p3" ,&fMc4p3[0] ,"Truth3[4]/D");
377  tAna3g->Branch("pos3" ,&fpos3[0] ,"pos3[3]/D");
378  tAna3g->Branch("e3" ,&fenergy3 ,"e3/D");
379  tAna3g->Branch("rawe3" ,&femcraw3 ,"rawe3/D");
380  tAna3g->Branch("nd3" ,&fndigi3 ,"NumOfDigi3/I");
381  tAna3g->Branch("mcTrk3" ,&fMcTrack3 ,"McTrackID3/I");
382  tAna3g->Branch("mod3" ,&fmod3 ,"mod3/I");
383  tAna3g->Branch("e253" ,&fe253 ,"e253/D");
384  tAna3g->Branch("lat3" ,&flat3 ,"lat3/D");
385  tAna3g->Branch("z203" ,&fz203 ,"z203/D");
386  tAna3g->Branch("z533" ,&fz533 ,"z533/D");
387 
388  }
389 
390 
391  evtset.clear();
392  pevtset.clear();
393  evtMap.clear();
394 
395  totNumOfWave = 0;
396  totNumOfSubWave = 0;
397  totNumOfDigi = 0;
398  totNumOfSharedDigi = 0;
399  totNumOfBump = 0;
400  totNumOfEvents = 0;
402  fNumOfMatch = 0;
403  fNumOfDismatch = 0;
404 
405  //vector resize
406  // const Int_t Num = 1;
407  // fEnergy.resize(Num);
408  // fEnergyC1.resize(Num);
409  // fEnergyC2.resize(Num);
410  // fPositionV.resize(3*Num);
411  // fp4.resize(4*Num);
412  // fMcp4.resize(4*Num);
413  // fZ201.resize(Num);
414  // fZ531.resize(Num);
415  // fLat1.resize(Num);
416  // fModule1.resize(Num);
417  // fTheta1.resize(Num);
418  // fPhi1.resize(Num);
419  // fDet1.resize(Num);
420  // fHitEmc1.resize(Num);
421  // fHitEmc2.resize(Num);
422  // fNPointV.resize(Num);
423 
424  cout << "-I- PndEmcAnalysis: Intialization successfull" << endl;
425  return kSUCCESS;
426 }
427 
428 void PndEmcAnalysis::Exec(Option_t*)
429 {
430  TStopwatch timer;
431  if (fVerbose>2){
432  timer.Start();
433  }
434 
435 
436  fEvtTime = FairRootManager::Instance()->GetEventTime();
437  cout<<"=======================PndEmcAnalysis======================="<<endl;
438  if(fSaveHits){
439  Int_t nHits = fEmcHitsArray->GetEntriesFast();
440  for(Int_t iHit =0; iHit < nHits; ++iHit){
441  PndEmcHit* theHit = (PndEmcHit*)fEmcHitsArray->At(iHit);
442  fHitTime = theHit->GetTime();
443  fModule = theHit->GetModule();
444  fHitEnergy = theHit->GetEnergy();
445  fTheta = theHit->GetTheta();
446  fPhi = theHit->GetPhi();
447  fNPoint = theHit->GetPointList().size();
448  fXcor = theHit->GetXPad();
449  fYcor = theHit->GetYPad();
450  tHit->Fill();
451  }
452  if(fVerbose>0)
453  cout<<"Save "<<nHits<<" PndEmcHit done"<<endl;
454  }
455  if(fTimeOrderedDigi){
456  Double_t time_length = 40.;//99.98%
457  if(FairRunAna::Instance()->IsTimeStamp()){
458  fWaveformArray = FairRootManager::Instance()->GetData("EmcSortedWaveform"
459  , fFunctor
460  , FairRootManager::Instance()->GetEventTime() + time_length);
461  }
462  }
463 
464  Int_t nWaveforms = fWaveformArray->GetEntriesFast();
465  totNumOfWave += nWaveforms;
466  if(fSaveWave){
467  for (Int_t iWaveform=0; iWaveform<nWaveforms; iWaveform++)
468  {
469  PndEmcWaveform* theWaveform = (PndEmcWaveform*) fWaveformArray->At(iWaveform);
470  if(fTimeOrderedDigi){
471  const std::vector<Int_t>& evtList = theWaveform->GetEvtList();
472  for(size_t i=0; i< evtList.size(); ++i){
473  evtset.insert(evtList[i]);
474  ++ evtMap[evtList[i]].first;//all waveforms counter
475  if(evtList.size() > 1){
476  pevtset.insert(evtList[i]);
477  ++ evtMap[evtList[i]].second;//pileup waveforms counter
478  }
479  }
480  }
481 
482  fPileupCount= theWaveform->GetPileupCount();
484  fHitIndex = theWaveform->GetHitIndex();
485  fModule = theWaveform->GetModule();
486  fWaveTime = theWaveform->GetTimeStamp();
487  fWavePeak = theWaveform->Max();
488  fWaveLen = theWaveform->GetWaveformLength();
489  fWaveIdx = iWaveform;
490  fBaseline = theWaveform->GetBaseline();
491  fXcor = theWaveform->GetTCI()->XCoord();
492  fYcor = theWaveform->GetTCI()->YCoord();
493  if(!fTimeOrderedDigi)
494  fHitEnergy = ((PndEmcHit*)fEmcHitsArray->At(fHitIndex))->GetEnergy();
495  tWave->Fill();
496  }
497  if(fVerbose>0)
498  cout<<"Save "<<nWaveforms<<" PndEmcWaveform done"<<endl;
499  }
500 
501  totDigiEnergy = 0.;
502  totSharedDigiEnergy = 0.;
503  Int_t nDigi = fDigiArray->GetEntriesFast();
504  Int_t nSharedDigi = fSharedDigiArray->GetEntriesFast();
505  if(!fTimeOrderedDigi){
506  totNumOfDigi += nDigi;
507  totNumOfSharedDigi += nSharedDigi;
508  }
509  if(fSaveDigi){
510  for (Int_t iDigi=0; iDigi<nDigi; iDigi++)
511  {
512  PndEmcDigi* theDigi = (PndEmcDigi*) fDigiArray->At(iDigi);
513  totDigiEnergy += theDigi->GetEnergy();
514  fDigiTime = theDigi->GetTimeStamp();
515  fModule = theDigi->GetModule();
516  fWaveIdx = theDigi->GetHitIndex();
517  if(!fTimeOrderedDigi){
518  fWaveTime = ((PndEmcWaveform*)fWaveformArray->At(fWaveIdx))->GetTimeStamp();
519  fHitEnergy = ((PndEmcHit*)fEmcHitsArray->At(fWaveIdx))->GetEnergy();
520  }
521  fDigiEnergy = theDigi->GetEnergy();
522  fDistance = theDigi->where().Mag();
523  fPosition = theDigi->where();
524  tDigi->Fill();
525  }
526  for (Int_t iDigi=0; iDigi<nSharedDigi; iDigi++)
527  {
528  PndEmcDigi* theDigi = (PndEmcDigi*) fSharedDigiArray->At(iDigi);
529  totSharedDigiEnergy += theDigi->GetEnergy();
530  }
531  if(fVerbose>0)
532  cout<<"Save "<<nDigi<<" PndEmcDigi done"<<endl;
533  }
534  Int_t nBump = fBumpArray->GetEntriesFast();
535  totNumOfBump += nBump;
536  std::vector<key> goodHits;
537  if(fSaveBump){
539  PndEmcClusterCalibrator::MakeEmcClusterCalibrator(2); //PndEmcAbsClusterCalibrator * calibrator2= //[R.K.03/2017] unused variable
540  for (Int_t iBump=0; iBump<nBump; iBump++)
541  {
542  PndEmcBump* theBump = (PndEmcBump*) fBumpArray->At(iBump);
543  if((fEnergy1C=calibrator1->Energy(theBump)) > 0.050){//50MeV
544  goodHits.push_back(key(iBump, fEnergy1C));
545  }
546  fModule = theBump->GetModule();
547  fBumpEnergy = theBump->energy();
548  fPosition = theBump->where();
549  fNumDigi = theBump->NumberOfDigis();
550  fNumBump = theBump->NBumps();
551  fClusterIdx = theBump->GetClusterIndex();
552  fNumOfSeed = theBump->LocalMaxMap().size();
553  fDigiMapSize = theBump->MemberDigiMap().size();
554  const std::vector<Int_t>& digis = theBump->DigiList();
555  if(fTimeOrderedDigi){
556  totNumOfDigi += digis.size();
557  totNumOfSharedDigi += digis.size();
558  }
559  PndEmcDigi* seedDigi(0);
560  Double_t maxE(-1.);
561  for(size_t id=0; id < digis.size(); ++id){
562  PndEmcDigi* theDigi = (PndEmcDigi*)fSharedDigiArray->At(digis[id]);
563  if(theDigi->GetEnergy() > maxE){
564  maxE = theDigi->GetEnergy();
565  seedDigi = theDigi;
566  }
567  theDigi->fEvtNo == fEventNo ? fNumOfMatch++ : fNumOfDismatch ++;
568  }
569  fSeedTime = seedDigi == 0? 0 :/* digiCalibrator.CalibrationEvtTimeByDigi(seedDigi, kFALSE);*/seedDigi->GetTimeStamp();
570  fBumpTime1 = theBump->GetTimeStamp();
571  //fBumpTime2 = theBump->GetTimeStamp1();
572  fSeedPosition = seedDigi->where();
573  //fBumpTime3 = theBump->GetTimeStamp3();
574  tBump->Fill();
575  }
576  if(fVerbose>0)
577  cout<<"Save "<<nBump<<" PndEmcBump done"<<endl;
578  }
579  if(fSaveReco){
580  std::vector<key> goodTrack;
581  //std::vector<IndexKey> indicator;
582  sort(goodHits.begin(), goodHits.end(), std::greater<key>());
583  if(goodHits.size() >=1){
584  //reset
585  for(Int_t idx=0; idx < NElement; ++idx){
587  = fZ201[idx] = fZ531[idx] = fLat1[idx] = fModule1[idx] = 0;
588  fPositionV[idx*3] = fPositionV[idx*3+1] = fPositionV[idx*3+2] = 0;
589  fp4[idx*4] = fp4[idx*4+1] = fp4[idx*4+2] = fp4[idx*4+3] = 0;
590  }
591  if(fMcTrackArray){
592  Int_t nTrack = fMcTrackArray->GetEntriesFast() > NElement ? NElement : fMcTrackArray->GetEntriesFast();
593  for(Int_t i=0; i < nTrack; ++i){
595  if(p1->IsGeneratorCreated()){
596  goodTrack.push_back(key(i, p1->Get4Momentum().E()));
597  //TLorentzVector p4(p1->GetMomentum(), p1->GetE());
598  //memcpy(&fMcp4[i*4], &p4[0], 4*sizeof(Double_t));
599  }else
600  break;
601  }
602  sort(goodTrack.begin(), goodTrack.end(), std::greater<key>());
603  for(size_t i=0; i< goodTrack.size(); ++i){
604  PndMCTrack* p1 = (PndMCTrack*)fMcTrackArray->At(goodTrack[i].iHit);
605  TLorentzVector p4(p1->GetMomentum(), p1->Get4Momentum().E());
606  memcpy(&fMcp4[i*4], &p4[0], 4*sizeof(Double_t));
607  fPDGCode[i] = p1->GetPdgCode();
608  }
609  }
610 
613  Int_t NEL = goodHits.size() > NElement ? NElement : goodHits.size();
615  for(size_t i=NEL;i<goodHits.size();++i){
616  PndEmcBump* theBump = (PndEmcBump*) fBumpArray->At(goodHits[i].iHit);
617  fRestEnergy1 += theBump->energy();
618  fRestEnergy2 += calibrator1->Energy(theBump);
619  fRestEnergy3 += calibrator1->Energy(theBump);
620  }
621  fNumBump = goodHits.size();
622  Int_t locIdx =0;
623 
624  for(Int_t idx=0; idx < NEL ; ++idx){
625  PndEmcBump* theBump = (PndEmcBump*)fBumpArray->At(goodHits[idx].iHit);
626  //PndEmcBump* theBump = (PndEmcBump*)fBumpArray->At(index);
627  //if(idx < goodTrack.size())
628  // locIdx = goodTrack[idx].iHit;
629  //else
630  // locIdx = idx;
631  locIdx = idx;
632  fEnergy[locIdx] = theBump->energy();
633  fEnergyC1[locIdx] = calibrator1->Energy(theBump)/1.018;
634  fEnergyC2[locIdx] = calibrator2->Energy(theBump);
635  TVector3 pos(calibrator1->Where(theBump));
636  memcpy(&fPositionV[locIdx*3], &pos[0], 3*sizeof(Double_t));
637  fTheta1[locIdx] = pos.Theta();
638  fPhi1[locIdx] = pos.Phi();
639  fZ201[locIdx] = theBump->Z20();
640  fZ531[locIdx] = theBump->Z53();
641  fLat1[locIdx] = theBump->LatMom();
642  fModule1[locIdx] = theBump->GetModule();
643  TLorentzVector p4(fEnergyC1[locIdx]*TMath::Sin(fTheta1[locIdx])*TMath::Cos(fPhi1[locIdx])
644  , fEnergyC1[locIdx]*TMath::Sin(fTheta1[locIdx])*TMath::Sin(fPhi1[locIdx])
645  , fEnergyC1[locIdx]*TMath::Cos(fTheta1[locIdx]), fEnergyC1[locIdx]);
646  memcpy(&fp4[locIdx*4], &p4[0], 4*sizeof(Double_t));
647  }
648  tRec->Fill();
649  if(fVerbose>0)
650  cout<<"Save "<<NEL<<" reconstruction data done"<<endl;
651  }
652  }
653  if(fSaveMcTruth){
654  Bool_t isFill = kFALSE;
655  Int_t nTrack = fMcTrackArray->GetEntriesFast();// > NElement ? NElement : fMcTrackArray->GetEntriesFast();
656  std::vector< std::map<Int_t, Int_t> > trackMap;//top node, and its tree
657  fNumGT = 0;
658  for(Int_t idx=0; idx < NElement; ++idx){
660  = fModule1[idx] = fHitEmc1[idx] = fHitEmc2[idx]= 0;
661  fPositionV[idx*3] = fPositionV[idx*3+1] = fPositionV[idx*3+2] = 0;
662  fMcp4[idx*4] = fMcp4[idx*4+1] = fMcp4[idx*4+2] = fMcp4[idx*4+3] = 0;
663  }
664  for(Int_t itrack = 0; itrack < nTrack; ++itrack){
665  PndMCTrack* pt1 =((PndMCTrack*)fMcTrackArray->At(itrack));
666  if(!pt1) break;
667  if(pt1->IsGeneratorCreated()){
668  pt1->Print(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);
672  isFill = kTRUE;
673  TLorentzVector p4(pt1->GetMomentum(), pt1->Get4Momentum().E());
674  memcpy(&fMcp4[itrack*4], &p4[0], 4*sizeof(Double_t));
675  fTheta1[itrack] = p4.Theta();
676  fPhi1[itrack] = p4.Phi();
677  fHitEmc1[itrack] = fMcTrackArray->GetEntriesFast();
678  fHitEmc2[itrack] = pt1->GetNPoints(kEMC);
679  fNPointV[itrack] = pt1->GetNPoints(kEMC);
680  //fDet1[itrack] = -1;
681  //for(Int_t i=0;i<15;++i){
682  // if(detIds[i] == kEMC){
683  // if(pt1->GetNPoints(detIds[i])>0)
684  // {fDet1[itrack] = detIds[i];break;}
685  // }else{
686  // if(pt1->GetNPoints(detIds[i])>0)
687  // {fDet1[itrack] = detIds[i];break;}
688  // }
689  //}
690  }else{
691  if( itrack < 6){
692  pt1->Print(itrack);
693  cout<<"MC points #"<<pt1->GetNPoints(kEMC)<<", #"<<trackMap.size()<<endl;
694  }
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()));
699  fNPointV[im] += pt1->GetNPoints(kEMC);
700  break;
701  }
702  }
703  if(pt1->GetNPoints(kEMC) > 0) ++ fNumGT;
704  }
705  }
706  if(isFill)
707  tMcTruth->Fill();
708  if(fVerbose>0)
709  cout<<"Save MCTruth data done"<<endl;
710  }
711  if(fSavePidCharged){
712  for(Int_t idx=0; idx < NElement; ++idx){
713  fchar[idx] = fmcidx[idx] = 0;
714  fmom[idx] = 0.;
715  fpos[idx*3] = fpos[idx*3+1] = fpos[idx*3+2] = 0;
716  fp4[idx*4] = fp4[idx*4+1] = fp4[idx*4+2] = fp4[idx*4+3] = 0;
717  }
718  Int_t nCand = fChargedCand->GetEntriesFast() > NElement ? NElement : fChargedCand->GetEntriesFast();
719  for(Int_t ic=0;ic < nCand; ++ic){
720  PndPidCandidate* theCand = (PndPidCandidate*)fChargedCand->At(ic);
721  fchar[ic] = theCand->GetCharge();
722  fmom[ic] = theCand->GetMomentum().Mag();
723  fmcidx[ic] = theCand->GetMcIndex();
724  TVector3 pos = theCand->GetPosition();
725  memcpy(&fpos[ic*3], &pos[0], 3*sizeof(Double_t));
726  TLorentzVector p4 = theCand->GetLorentzVector();
727  memcpy(&fp4[ic*4], &p4[0], 4*sizeof(Double_t));
728  }
729  tCharged->Fill();
730  if(fVerbose>0)
731  cout<<"Save "<<nCand<<" charged candidates done"<<endl;
732  }
733  if(fSavePidNeutral){
734  Int_t nCand = fNeutralCand->GetEntriesFast() > NElement ? NElement : fNeutralCand->GetEntriesFast();
735  for(Int_t idx=0; idx < NElement; ++idx){
736  fchar[idx] = fmcidx[idx] = fnbump[idx] = fndigi[idx] = fmod[idx] = fcluIdx[idx] = 0;
737  frawe[idx] = fcale[idx] = fe25[idx] = fe9[idx] = 0;
738  fpos[idx*3] = fpos[idx*3+1] = fpos[idx*3+2] = 0;
739  fp4[idx*4] = fp4[idx*4+1] = fp4[idx*4+2] = fp4[idx*4+3] = 0;
740  }
741  for(Int_t ic=0;ic < nCand; ++ic){
742  PndPidCandidate* theCand = (PndPidCandidate*)fNeutralCand->At(ic);
743  fchar[ic] = theCand->GetCharge();
744  frawe[ic] = theCand->GetEmcRawEnergy();
745  fcale[ic] = theCand->GetEmcCalEnergy();
746  fndigi[ic]= theCand->GetEmcNumberOfCrystals();
747  fnbump[ic] = theCand->GetEmcNumberOfBumps();
748  fmod[ic] = theCand->GetEmcModule();
749  fcluIdx[ic] = theCand->GetEmcIndex();
750  fe25[ic] = theCand->GetEmcClusterE25();
751  fe9[ic] = theCand->GetEmcClusterE9();
752  TVector3 pos = theCand->GetPosition();
753  memcpy(&fpos[ic*3], &pos[0], 3*sizeof(Double_t));
754  TLorentzVector p4 = theCand->GetLorentzVector();
755  memcpy(&fp4[ic*4], &p4[0], 4*sizeof(Double_t));
756  }
757  tNeutral->Fill();
758  if(fVerbose>0)
759  cout<<"Save "<<nCand<<" neutral candidates done"<<endl;
760  }
761  if(fSaveJpsiTo3g && fChargedCand->GetEntriesFast() == 0){
762  Int_t nCand = fNeutralCand->GetEntriesFast();
763  if(fVerbose>0)
764  cout<<"begin write #"<<nCand<<" candidates."<<endl;
765  std::vector<key> goodCand;
766  for(Int_t ic=0;ic<nCand;++ic){
767  PndPidCandidate* theCand = (PndPidCandidate*)fNeutralCand->At(ic);
768  if(theCand->GetEnergy() > 0.02)
769  goodCand.push_back(key(ic, theCand->GetEnergy()));
770  }
771  if(goodCand.size() >= 3){
772  sort(goodCand.begin(), goodCand.end(), std::greater<key>());
773  PndPidCandidate* cand1 = (PndPidCandidate*)fNeutralCand->At(goodCand[0].iHit);
774  fenergy1 = cand1->GetEnergy();
775  femcraw1 = cand1->GetEmcRawEnergy();
776  fndigi1 = cand1->GetEmcNumberOfCrystals();
777  fmod1 = cand1->GetEmcModule();
778  fMcTrack1 = cand1->GetMcIndex();
779  fe251 = cand1->GetEmcClusterE25();
780  flat1 = cand1->GetEmcClusterLat();
781  fz201 = cand1->GetEmcClusterZ20();
782  fz531 = cand1->GetEmcClusterZ53();
783  TLorentzVector p4 = cand1->GetLorentzVector();
784  TVector3 pos = cand1->GetPosition();
785  memcpy(&f4p1[0], &p4[0], 4*sizeof(Double_t));
786  memcpy(&fpos1[0], &pos[0], 3*sizeof(Double_t));
787  //f4p1 = cand1->GetLorentzVector();
788  //fpos1 = cand1->GetPosition();
789  PndPidCandidate* cand2 = (PndPidCandidate*)fNeutralCand->At(goodCand[1].iHit);
790  fenergy2 = cand2->GetEnergy();
791  femcraw2 = cand2->GetEmcRawEnergy();
792  fndigi2 = cand2->GetEmcNumberOfCrystals();
793  fmod2 = cand2->GetEmcModule();
794  fe252 = cand2->GetEmcClusterE25();
795  flat2 = cand2->GetEmcClusterLat();
796  fz202 = cand2->GetEmcClusterZ20();
797  fz532 = cand2->GetEmcClusterZ53();
798  fMcTrack2 = cand2->GetMcIndex();
799  p4 = cand2->GetLorentzVector();
800  pos = cand2->GetPosition();
801  memcpy(&f4p2[0], &p4[0], 4*sizeof(Double_t));
802  memcpy(&fpos2[0], &pos[0], 3*sizeof(Double_t));
803  //f4p2 = cand2->GetLorentzVector();
804  //fpos2 = cand2->GetPosition();
805  PndPidCandidate* cand3 = (PndPidCandidate*)fNeutralCand->At(goodCand[2].iHit);
806  fenergy3 = cand3->GetEnergy();
807  femcraw3 = cand3->GetEmcRawEnergy();
808  fndigi3 = cand3->GetEmcNumberOfCrystals();
809  fmod3 = cand3->GetEmcModule();
810  fe253 = cand3->GetEmcClusterE25();
811  flat3 = cand3->GetEmcClusterLat();
812  fz203 = cand3->GetEmcClusterZ20();
813  fz533 = cand3->GetEmcClusterZ53();
814  fMcTrack3 = cand3->GetMcIndex();
815  p4 = cand3->GetLorentzVector();
816  pos = cand3->GetPosition();
817  memcpy(&f4p3[0], &p4[0], 4*sizeof(Double_t));
818  memcpy(&fpos3[0], &pos[0], 3*sizeof(Double_t));
819  //f4p3 = cand3->GetLorentzVector();
820  //fpos3 = cand3->GetPosition();
821  fNumberOfGoodPhoton = goodCand.size();
822  fExternalEnergy = 0.;
823  for(size_t jc=3;jc<goodCand.size();++jc){
824  PndPidCandidate* cand = (PndPidCandidate*)fNeutralCand->At(goodCand[jc].iHit);
825  fExternalEnergy += cand->GetEnergy();
826  }
827  memset(&fMc4p1[0], 0, 4*sizeof(Double_t));
828  memset(&fMc4p2[0], 0, 4*sizeof(Double_t));
829  memset(&fMc4p3[0], 0, 4*sizeof(Double_t));
830  if(fMcTrackArray){
831  if(fMcTrack1 >= 0){
833  p4 = TLorentzVector(p1->GetMomentum(), p1->Get4Momentum().E());
834  memcpy(&fMc4p1[0], &p4[0], 4*sizeof(Double_t));
835  }
836  if(fMcTrack2 >= 0){
838  p4 = TLorentzVector(p2->GetMomentum(), p2->Get4Momentum().E());
839  memcpy(&fMc4p2[0], &p4[0], 4*sizeof(Double_t));
840  }
841  if(fMcTrack3 >= 0){
843  p4 = TLorentzVector(p3->GetMomentum(), p3->Get4Momentum().E());
844  memcpy(&fMc4p3[0], &p4[0], 4*sizeof(Double_t));
845  }
846  }
847  tAna3g->Fill();
848  }
849  }
850 
851 
852  if(fSaveTask){
853  fNumOfHits = fEmcHitsArray->GetEntriesFast();
854  fNumWaveform = fWaveformArray->GetEntriesFast();
855  fNumDigi = fDigiArray->GetEntriesFast();
856  fNumSharedDigi= fSharedDigiArray->GetEntriesFast();
857  fNumCluster = fClusterArray->GetEntriesFast();
858  fNumBump = fBumpArray->GetEntriesFast();
859  fNumMcTrack = fMcTrackArray->GetEntriesFast();
860  tTask->Fill();
861  if(fVerbose>0)
862  cout<<"Save task data done"<<endl;
863  }
864 
865  if (fVerbose>2){
866  timer.Stop();
867  Double_t rtime = timer.RealTime();
868  Double_t ctime = timer.CpuTime();
869  cout << "PndEmcAnalysis, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
870  }
871  //fEvtHeaderArray->Delete();
872  fEventNo ++;
873  cout<<"================================================================="<<endl;
874 }
875 
877 
878  // Get run and runtime database
879  FairRun* run = FairRun::Instance();
880  if ( ! run ) Fatal("SetParContainers", "No analysis run");
881 
882  FairRuntimeDb* db = run->GetRuntimeDb();
883  if ( ! db ) Fatal("SetParContainers", "No runtime database");
884 
885  // Get Emc digitisation parameter container
886  //fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
887 
888  // Get Emc reconstruction parameter container
889  //fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");
890 
891  // Get Emc FPGA parameter container
892  //fFpgaPar = (PndEmcFpgaPar*) db->getContainer("PndEmcFpgaPar");
893 }
894 
896 {
897 
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;
902  std::cout<<"Total subwaveforms# "<<totNumOfSubWave<<std::endl;
903  std::cout<<"Total fpga digis# "<<totNumOfDigi<<std::endl;
904  std::cout<<"Total shared digis# "<<totNumOfSharedDigi<<std::endl;
905  std::cout<<"Total bumps# "<<totNumOfBump<<std::endl;
906  //std::cout<<"Digis to events they belong to#"<<matchCount<<std::endl;
907  //std::cout<<"Digis to events they dont belong to#"<<unmatchCount<<std::endl;
908  std::cout<<"*********************************************************"<<std::endl;
909  if(fTimeOrderedDigi){
910  std::map<Int_t, std::pair<Int_t,Int_t> >::iterator it = evtMap.begin();
911  for(; it != evtMap.end(); ++it){
912  fEventNo = (*it).first;
913  fTotWave = (*it).second.first;
914  fPileupWave = (*it).second.second;
915  tEvtPileup->Fill();
916  }
917  totNumOfEvents = evtset.size();
918  totNumOfPileupEvents = pevtset.size();
919  tTaskPileup->Fill();
920  }
921 
922 
923  if(fStoreRooTFile){
924  TFile* gOldFile = gFile;
925  gFile = fRootFile;
926  gFile->cd();
927  if(fSaveHits){
928  cout<<"write "<<tHit->GetTitle()<<endl;
929  tHit->Write();
930  }
931  if(fSaveWave){
932  cout<<"write "<<tWave->GetTitle()<<endl;
933  tWave->Write();
934  }
935  if(fSaveDigi){
936  cout<<"write "<<tDigi->GetTitle()<<endl;
937  tDigi->Write();
938  }
939  if(fSaveBump){
940  cout<<"write "<<tBump->GetTitle()<<endl;
941  tBump->Write();
942  }
943  if(fSaveTask){
944  cout<<"write "<<tTask->GetTitle()<<endl;
945  tTask->Write();
946  }
947  if(fSaveReco){
948  cout<<"write "<<tRec->GetTitle()<<endl;
949  tRec->Write();
950  }
951  if(fSaveMcTruth){
952  cout<<"write "<<tMcTruth->GetTitle()<<endl;
953  tMcTruth->Write();
954  }
955  if(fTimeOrderedDigi){
956  cout<<"write "<<tTaskPileup->GetTitle()<<endl;
957  cout<<"write "<<tEvtPileup->GetTitle()<<endl;
958  tTaskPileup->Write();
959  tEvtPileup->Write();
960  }
961  if(fSavePidCharged){
962  cout<<"write "<<tCharged->GetTitle()<<endl;
963  tCharged->Write();
964  }
965  if(fSavePidNeutral){
966  cout<<"write "<<tNeutral->GetTitle()<<endl;
967  tNeutral->Write();
968  }
969  if(fSaveJpsiTo3g){
970  cout<<"write "<<tAna3g->GetTitle()<<endl;
971  tAna3g->Write();
972  }
973 
974  gFile = gOldFile;
975  }
976 
977  evtset.clear();
978  pevtset.clear();
979  evtMap.clear();
980 }
981 
key(Int_t hit, Double_t e)
Double_t fBumpEnergy
Double_t GetEmcClusterE25() const
TVector3 pos
bool operator<(const key &rhs) const
Double_t fRestEnergy1
TClonesArray * fEmcHitsArray
int fVerbose
Definition: poormantracks.C:24
Int_t GetCharge() const
Bool_t fSavePidNeutral
Short_t GetXPad() const
Definition: PndEmcHit.cxx:87
Int_t fHitEmc1[NElement]
Double_t fExternalEnergy
Int_t run
Definition: autocutx.C:47
Float_t GetPhi() const
Definition: PndEmcHit.h:57
TVector3 where() const
virtual Double_t GetEnergy() const
Definition: PndEmcDigi.cxx:296
represents the reconstructed hit of one emc crystal
Definition: PndEmcDigi.h:40
Int_t fchar[NElement]
Double_t totDigiEnergy
Int_t i
Definition: run_full.C:25
Float_t GetEmcRawEnergy() const
Short_t GetModule() const
Int_t fcluIdx[NElement]
Int_t fNPointV[NElement]
Double_t fpos[NElement]
Int_t totNumOfPileupEvents
Int_t GetNPoints(DetectorId detId) const
Definition: PndMCTrack.cxx:120
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
#define verbose
static const Int_t NElement
TClonesArray * fMcTrackArray
Float_t GetEmcCalEnergy() const
std::map< Int_t, std::pair< Int_t, Int_t > > evtMap
static T Sin(const T &x)
Definition: PndCAMath.h:42
std::set< Int_t > pevtset
Int_t GetEmcNumberOfCrystals() const
TLorentzVector GetLorentzVector() const
TClonesArray * fBumpArray
Int_t GetHitIndex() const
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
Double_t GetBaseline() 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
Definition: PndEmcCluster.h:40
Double_t GetEmcClusterE9() const
virtual InitStatus Init()
Int_t fDet1[NElement]
Double_t fe9[NElement]
TClonesArray * fClusterArray
Double_t fcale[NElement]
Double_t fBumpTime3
Double_t fBaseline
BinaryFunctor * fFunctor
Float_t GetTheta() const
Definition: PndEmcHit.h:56
static T Cos(const T &x)
Definition: PndCAMath.h:43
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
Double_t fe25[NElement]
TLorentzVector Get4Momentum() const
Definition: PndMCTrack.cxx:102
bool operator>(const key &rhs) const
Short_t GetModule() const
Definition: PndEmcDigi.h:103
Double_t fEnergyC1[NElement]
Double_t fEnergyC2[NElement]
Int_t GetMcIndex() const
Double_t fmom[NElement]
Int_t GetEmcIndex() const
Double_t fMc4p3[4]
Double_t fHitEnergy
Double_t fEnergy1C
int idx[MAX]
Definition: autocutx.C:38
Double_t fDigiTime
std::vector< PndEmcPoint * > & GetPointList()
Definition: PndEmcHit.h:68
PndEmcAnalysis(Int_t verbose=0, Bool_t storeRootFile=kTRUE)
Double_t fSeedTime
virtual ~PndEmcAnalysis()
Double_t fMcp4[4 *NElement]
Short_t GetYPad() const
Definition: PndEmcHit.cxx:120
Double_t f4p3[4]
int nHits
Definition: RiemannTest.C:16
std::set< Int_t > evtset
Double_t f4p1[4]
Double_t Z53() const
Definition: PndEmcCluster.h:69
Double_t Z20() const
Definition: PndEmcCluster.h:67
Double_t
Double_t fpos3[3]
TClonesArray * fNeutralCand
Double_t f4p2[4]
PndEmcTwoCoordIndex * GetTCI() const
virtual Double_t GetEnergy() const
Definition: PndEmcHit.h:54
Double_t fBumpTime1
Bool_t fSaveJpsiTo3g
Int_t fNumberOfGoodPhoton
const std::map< Int_t, Int_t > & LocalMaxMap() const
Definition: PndEmcCluster.h:42
Double_t GetEmcClusterZ53() const
const std::vector< Int_t > & GetEvtList() const
TStopwatch timer
Definition: hit_dirc.C:51
Double_t fpos2[3]
Double_t GetEmcClusterLat() const
Int_t fnbump[NElement]
Double_t LatMom() const
Definition: PndEmcCluster.h:71
void Print(Int_t iTrack=0) const
Definition: PndMCTrack.cxx:95
static PndEmcAbsClusterCalibrator * MakeEmcClusterCalibrator(Int_t method, Int_t version=1)
TClonesArray * fDigiArray
Int_t GetEmcNumberOfBumps() const
Double_t GetEmcClusterZ20() const
Double_t dE
Double_t fTheta1[NElement]
Double_t fRestEnergy3
Int_t fHitEmc2[NElement]
Int_t GetWaveformLength() const
virtual void FinishTask()
Double_t fZ201[NElement]
Double_t fPositionV[3 *NElement]
const std::map< Int_t, Int_t > & MemberDigiMap() const
Definition: PndEmcCluster.h:43
Double_t frawe[NElement]
virtual TVector3 Where(PndEmcCluster *clust, Int_t pid=22)=0
represents a simulated waveform in an emc crystal
Double_t ctime
Definition: hit_dirc.C:114
TPad * p2
Definition: hist-t7.C:117
virtual Double_t GetTime() const
Definition: PndEmcHit.h:55
TClonesArray * fSharedDigiArray
Double_t fp4[4 *NElement]
represents the deposited energy of one emc crystal from simulation
Definition: PndEmcHit.h:26
Int_t NumberOfDigis() const
Int_t fEvtNo
Definition: PndEmcDigi.h:118
Int_t NBumps() const
TVector3 fPosition
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)
Definition: hit.C:1
Int_t GetClusterIndex()
Definition: PndEmcBump.h:54
Double_t fPhi1[NElement]
Int_t GetHitIndex()
Definition: PndEmcDigi.h:110
virtual void Exec(Option_t *opt)
TVector3 GetPosition() const
Int_t fmod[NElement]
Int_t GetPileupCount() const
Int_t iHit
Double_t fMc4p1[4]
virtual Double_t energy() const
ClassImp(PndAnaContFact)
static const Int_t NSubDetectors
TPad * p1
Definition: hist-t7.C:116
const TVector3 & where() const
Definition: PndEmcDigi.h:111
TVector3 fSeedPosition
Short_t GetModule() const
Double_t fWavePeak
TClonesArray * fChargedCand
Double_t fDistance
Double_t fEnergy[NElement]
Bool_t fTimeOrderedDigi
set to kTRUE to use the time ordering of the output data.
Int_t fPDGCode[NElement]
Bool_t fSavePidCharged
Double_t fRestEnergy2
Int_t fmcidx[NElement]
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
Double_t rtime
Definition: hit_dirc.C:113
TVector3 GetMomentum() const
TClonesArray * fWaveformArray
Int_t fModule1[NElement]
Short_t GetModule() const
Definition: PndEmcHit.h:58
Int_t totNumOfSharedDigi
Int_t GetEmcModule() const
Double_t fZ531[NElement]
Double_t fDigiEnergy
Int_t fndigi[NElement]
represents a reconstructed (splitted) emc cluster
Definition: PndEmcBump.h:34
Double_t fMc4p2[4]
Double_t fBumpTime2
Double_t totSharedDigiEnergy
Double_t fpos1[3]
Double_t fWaveTime
Double_t fLat1[NElement]