Add more AODPid information from tracking detectors when filtering ESDs
[u/mrichter/AliRoot.git] / PWG4 / PartCorrBase / AliCaloTrackReader.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id:  $ */
16
17 //_________________________________________________________________________
18 // Base class for reading data: MonteCarlo, ESD or AOD, of PHOS EMCAL and 
19 // Central Barrel Tracking detectors (CTS).
20 // Not all MC particles/tracks/clusters are kept, some kinematical/fidutial restrictions are done.
21 // Mother class of : AliCaloTrackESDReader: Fills ESD data in 3 TObjArrays (PHOS, EMCAL, CTS)
22 //                 : AliCaloTrackMCReader: Fills Kinematics data in 3 TObjArrays (PHOS, EMCAL, CTS)
23 //                 : AliCaloTrackAODReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS) 
24 //                
25 //-- Author: Gustavo Conesa (LNF-INFN) 
26 //////////////////////////////////////////////////////////////////////////////
27
28
29 // --- ROOT system ---
30 #include <TObjArray.h>
31 #include "TFile.h"
32
33 //---- ANALYSIS system ----
34 #include "AliCaloTrackReader.h"
35 #include "AliFidutialCut.h"
36 #include "AliMCEvent.h"
37 #include "AliAODMCHeader.h"
38 #include "AliGenPythiaEventHeader.h"
39 #include "AliAODMCParticle.h"
40
41 ClassImp(AliCaloTrackReader)
42   
43   
44 //____________________________________________________________________________
45   AliCaloTrackReader::AliCaloTrackReader() : 
46     TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0), 
47     fFidutialCut(0x0), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
48     fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0),
49     fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
50     fEMCALCells(0x0), fPHOSCells(0x0),
51     fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
52     fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
53     fFillEMCALCells(0),fFillPHOSCells(0), 
54         fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
55         fSecondInputFileName(""),fSecondInputFirstEvent(0), 
56         fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0), 
57     fAODPHOSNormalInputEntries(0), fTrackStatus(0), 
58         fReadStack(kFALSE), fReadAODMCParticles(kFALSE)
59 {
60   //Ctor
61   
62   //Initialize parameters
63   InitParameters();
64 }
65
66 //____________________________________________________________________________
67 AliCaloTrackReader::AliCaloTrackReader(const AliCaloTrackReader & g) :   
68   TObject(g), fEventNumber(g.fEventNumber), fCurrentFileName(g.fCurrentFileName), 
69   fDataType(g.fDataType), fDebug(g.fDebug),
70   fFidutialCut(g.fFidutialCut),
71   fComparePtHardAndJetPt(g.fComparePtHardAndJetPt),
72   fPtHardAndJetPtFactor(g.fPtHardAndJetPtFactor),
73   fCTSPtMin(g.fCTSPtMin), fEMCALPtMin(g.fEMCALPtMin),fPHOSPtMin(g.fPHOSPtMin), 
74   fAODCTS(new TObjArray(*g.fAODCTS)),  
75   fAODEMCAL(new TObjArray(*g.fAODEMCAL)),
76   fAODPHOS(new TObjArray(*g.fAODPHOS)),
77   fEMCALCells(new TNamed(*g.fEMCALCells)),
78   fPHOSCells(new TNamed(*g.fPHOSCells)),
79   fInputEvent(g.fInputEvent), fOutputEvent(g.fOutputEvent), fMC(g.fMC),
80   fFillCTS(g.fFillCTS),fFillEMCAL(g.fFillEMCAL),fFillPHOS(g.fFillPHOS),
81   fFillEMCALCells(g.fFillEMCALCells),fFillPHOSCells(g.fFillPHOSCells),
82   fSecondInputAODTree(g.fSecondInputAODTree), 
83   fSecondInputAODEvent(g.fSecondInputAODEvent),
84   fSecondInputFileName(g.fSecondInputFileName), 
85   fSecondInputFirstEvent(g.fSecondInputFirstEvent),
86   fAODCTSNormalInputEntries(g.fAODCTSNormalInputEntries), 
87   fAODEMCALNormalInputEntries(g.fAODEMCALNormalInputEntries), 
88   fAODPHOSNormalInputEntries(g.fAODPHOSNormalInputEntries),
89   fTrackStatus(g.fTrackStatus),
90   fReadStack(g.fReadStack), fReadAODMCParticles(g.fReadAODMCParticles)
91 {
92   // cpy ctor
93   
94 }
95
96 //_________________________________________________________________________
97 AliCaloTrackReader & AliCaloTrackReader::operator = (const AliCaloTrackReader & source)
98 {
99   // assignment operator
100   
101   if(&source == this) return *this;
102   
103   fDataType    = source.fDataType ;
104   fDebug       = source.fDebug ;
105   fEventNumber = source.fEventNumber ;
106   fCurrentFileName = source.fCurrentFileName ;
107   fFidutialCut = source.fFidutialCut;
108         
109   fComparePtHardAndJetPt = source.fComparePtHardAndJetPt;
110   fPtHardAndJetPtFactor  = source.fPtHardAndJetPtFactor;
111         
112   fCTSPtMin    = source.fCTSPtMin ;
113   fEMCALPtMin  = source.fEMCALPtMin ;
114   fPHOSPtMin   = source.fPHOSPtMin ; 
115   
116   fAODCTS     = new TObjArray(*source.fAODCTS) ;
117   fAODEMCAL   = new TObjArray(*source.fAODEMCAL) ;
118   fAODPHOS    = new TObjArray(*source.fAODPHOS) ;
119   fEMCALCells = new TNamed(*source.fEMCALCells) ;
120   fPHOSCells  = new TNamed(*source.fPHOSCells) ;
121
122   fInputEvent  = source.fInputEvent;
123   fOutputEvent = source.fOutputEvent;
124   fMC          = source.fMC;
125   
126   fFillCTS        = source.fFillCTS;
127   fFillEMCAL      = source.fFillEMCAL;
128   fFillPHOS       = source.fFillPHOS;
129   fFillEMCALCells = source.fFillEMCALCells;
130   fFillPHOSCells  = source.fFillPHOSCells;
131
132   fSecondInputAODTree    = source.fSecondInputAODTree;
133   fSecondInputAODEvent   = source.fSecondInputAODEvent;
134   fSecondInputFileName   = source.fSecondInputFileName;
135   fSecondInputFirstEvent = source.fSecondInputFirstEvent;
136
137   fAODCTSNormalInputEntries   = source.fAODCTSNormalInputEntries; 
138   fAODEMCALNormalInputEntries = source.fAODEMCALNormalInputEntries; 
139   fAODPHOSNormalInputEntries  = source.fAODPHOSNormalInputEntries;
140         
141   fTrackStatus        = source.fTrackStatus;
142   fReadStack          = source.fReadStack;
143   fReadAODMCParticles = source.fReadAODMCParticles;     
144         
145   return *this;
146   
147 }
148
149 //_________________________________
150 AliCaloTrackReader::~AliCaloTrackReader() {
151   //Dtor
152   
153   if(fFidutialCut) delete fFidutialCut ;
154         
155   if(fAODCTS){
156     fAODCTS->Clear() ; 
157     delete fAODCTS ;
158   }
159   
160   if(fAODEMCAL){
161     fAODEMCAL->Clear() ; 
162     delete fAODEMCAL ;
163   }
164   
165   if(fAODPHOS){
166     fAODPHOS->Clear() ; 
167     delete fAODPHOS ;
168   }
169   
170   if(fEMCALCells){
171     fEMCALCells->Clear() ; 
172     delete fEMCALCells ;
173   }
174   
175   if(fPHOSCells){
176     fPHOSCells->Clear() ; 
177     delete fPHOSCells ;
178   }
179
180   if(fInputEvent)  delete fInputEvent ;
181   if(fOutputEvent) delete fOutputEvent ;
182   if(fMC)          delete fMC ;  
183         
184   if(fSecondInputAODTree){
185         fSecondInputAODTree->Clear();
186         delete fSecondInputAODTree;
187   }
188         
189   if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
190
191 }
192
193
194 //_________________________________________________________________________
195 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
196         // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
197         // Only for PYTHIA.
198         if(!fReadStack) return kTRUE; //Information not filtered to AOD
199         
200         if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
201                 TParticle * jet =  new TParticle;
202                 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
203                 Int_t nTriggerJets =  pygeh->NTriggerJets();
204                 Float_t ptHard = pygeh->GetPtHard();
205                 
206                 //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
207             Float_t tmpjet[]={0,0,0,0};
208                 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
209                         pygeh->TriggerJet(ijet, tmpjet);
210                         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
211                         //Compare jet pT and pt Hard
212                         //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
213                         if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
214                                 printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
215                                            nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
216                                 return kFALSE;
217                         }
218                 }
219         }
220         
221         return kTRUE ;
222         
223 }
224
225 //____________________________________________________________________________
226 AliStack* AliCaloTrackReader::GetStack() const {
227   //Return pointer to stack
228   if(fMC)
229     return fMC->Stack();
230   else{
231     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n"); 
232     return 0x0 ;
233   }
234 }
235
236 //____________________________________________________________________________
237 AliHeader* AliCaloTrackReader::GetHeader() const {
238   //Return pointer to header
239   if(fMC)
240     return fMC->Header();
241   else{
242     printf("AliCaloTrackReader::Header is not available\n"); 
243     return 0x0 ;
244   }
245 }
246 //____________________________________________________________________________
247 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
248   //Return pointer to Generated event header
249   if(fMC)
250     return fMC->GenEventHeader();
251   else{
252     printf("AliCaloTrackReader::GenEventHeader is not available\n"); 
253     return 0x0 ;
254   }
255 }
256
257 //____________________________________________________________________________
258 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
259         //Return list of particles in AOD. Do it for the corresponding input event.
260         if(fDataType == kAOD){
261          //Normal input AOD
262          if(input == 0) return (TClonesArray*)((AliAODEvent*)fInputEvent)->FindListObject("mcparticles");
263           //Second input AOD
264          else if(input == 1 && fSecondInputAODEvent) return (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");        
265          else {
266              printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input); 
267                  return 0x0;
268          }
269         }
270         else {
271                 printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n"); 
272                 return 0x0;
273         }
274 }
275
276 //____________________________________________________________________________
277 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
278         //Return MC header in AOD. Do it for the corresponding input event.
279         if(fDataType == kAOD){
280                 //Normal input AOD
281                 if(input == 0) return (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
282                 //Second input AOD
283                 else if(input == 1) return  (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader"); 
284                 else {
285                         printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
286                         return 0x0;
287                 }
288         }
289         else {
290                 printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
291                 return 0x0;
292         }
293 }
294
295 //_______________________________________________________________
296 void AliCaloTrackReader::Init()
297 {
298         //Init reader. Method to be called in AliAnaPartCorrMaker
299         
300         //Get the file with second input events if the filename is given
301         //Get the tree and connect the AODEvent. Only with AODs
302         
303         if(fReadStack && fReadAODMCParticles){
304                 printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
305                 fReadStack = kFALSE;
306                 fReadAODMCParticles = kFALSE;
307         }
308         
309         if(fReadStack && !fMC) {
310                 printf("AliCaloTrackReader::Init() - MC handler not available, switch off stack reading\n");
311                 fReadStack = kFALSE;
312         }
313         
314         if(fSecondInputFileName!=""){
315                 if(fDataType == kAOD){
316                         TFile * input2   = new TFile(fSecondInputFileName,"read");
317                         printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
318                         fSecondInputAODTree = (TTree*) input2->Get("aodTree");
319                         if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n", 
320                                                                                    fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
321                         else{
322                          printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
323                          abort();
324                         }
325                         fSecondInputAODEvent = new AliAODEvent;
326                         fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
327                         if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
328                                 printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n", 
329                                            fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
330                                 abort();
331                         }
332                 }
333                 else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
334         }
335         
336         
337 }
338 //_______________________________________________________________
339 void AliCaloTrackReader::InitParameters()
340 {
341  
342   //Initialize the parameters of the analysis.
343   fDataType   = kESD ;
344   fCTSPtMin   = 0.2 ;
345   fEMCALPtMin = 0.2 ;
346   fPHOSPtMin  = 0.2 ;
347
348   fFillEMCAL      = kTRUE;
349   fFillPHOS       = kTRUE;
350   fFillCTS        = kTRUE;
351   fFillEMCALCells = kFALSE;
352   fFillPHOSCells  = kFALSE;
353
354   fFidutialCut           = new AliFidutialCut();
355   fSecondInputFileName   = "" ;
356   fSecondInputFirstEvent = 0 ;
357   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
358   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
359         
360 }
361
362
363 //________________________________________________________________
364 void AliCaloTrackReader::Print(const Option_t * opt) const
365 {
366
367   //Print some relevant parameters set for the analysis
368   if(! opt)
369     return;
370
371   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
372   printf("Data type      : %d\n", fDataType) ;
373   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
374   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
375   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
376   printf("Use CTS         =     %d\n", fFillCTS) ;
377   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
378   printf("Use PHOS        =     %d\n", fFillPHOS) ;
379   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
380   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
381   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
382   if(fComparePtHardAndJetPt)
383           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
384         
385   if(fSecondInputFileName!="") {
386           printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
387           printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
388   }
389         
390          printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
391         
392   printf("    \n") ;
393
394
395 //___________________________________________________
396 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
397   //Fill the event counter and input lists that are needed, called by the analysis maker.
398   
399   fEventNumber = iEntry;
400   fCurrentFileName = TString(currentFileName);
401         
402   if((fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){
403     printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n");
404     abort();
405   }
406         
407   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard 
408   if(fComparePtHardAndJetPt && GetStack()) {
409                 if(!ComparePtHardAndJetPt()) return kFALSE ;
410   }
411
412   //In case of mixing events with other AOD file        
413   if(fDataType == kAOD && fSecondInputAODTree){
414          
415          if(fDebug > 1) 
416                  printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
417          if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
418                  if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent) 
419                          printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
420                  return kFALSE;
421          }
422           
423          //Get the Event
424          Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
425          if ( nbytes == 0 ) {//If nothing in AOD
426                  printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
427                  abort() ; 
428          }
429           
430   }
431   //if(iEntry > 10) return kFALSE;
432   if(fFillCTS)   FillInputCTS();
433   if(fFillEMCAL) FillInputEMCAL();
434   if(fFillPHOS)  FillInputPHOS();
435   if(fFillEMCALCells) FillInputEMCALCells();
436   if(fFillPHOSCells)  FillInputPHOSCells();
437
438   return kTRUE ;
439 }
440
441 //__________________________________________________
442 void AliCaloTrackReader::ResetLists() {
443   //  Reset lists, called by the analysis maker 
444
445   if(fAODCTS)   fAODCTS -> Clear();
446   if(fAODEMCAL) fAODEMCAL -> Clear();
447   if(fAODPHOS)  fAODPHOS -> Clear();
448   if(fEMCALCells) fEMCALCells -> Clear();
449   if(fPHOSCells)  fPHOSCells -> Clear();
450
451 }