0774eaa23a90aaa9bfb45751001b033a59436cbb
[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/fiducial 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 //                 : AliCaloTrackReader: Fills AOD data in 3 TObjArrays (PHOS, EMCAL, CTS) 
24 //              
25 // This part is commented: Mixing analysis can be done, input AOD with events
26 // is opened in the AliCaloTrackReader::Init()
27
28 //-- Author: Gustavo Conesa (LNF-INFN) 
29 //////////////////////////////////////////////////////////////////////////////
30
31
32 // --- ROOT system ---
33 #include "TFile.h"
34
35 //---- ANALYSIS system ----
36 #include "AliCaloTrackReader.h"
37 #include "AliMCEvent.h"
38 #include "AliAODMCHeader.h"
39 #include "AliGenPythiaEventHeader.h"
40 #include "AliVEvent.h"
41 #include "AliAODEvent.h"
42 #include "AliVTrack.h"
43 #include "AliVParticle.h"
44 #include "AliMixedEvent.h"
45 #include "AliESDtrack.h"
46 #include "AliEMCALRecoUtils.h"
47 #include "AliESDtrackCuts.h"
48
49 ClassImp(AliCaloTrackReader)
50   
51   
52 //____________________________________________________________________________
53   AliCaloTrackReader::AliCaloTrackReader() : 
54     TObject(), fEventNumber(-1), fCurrentFileName(""),fDataType(0), fDebug(0), 
55     fFiducialCut(0x0), fCheckFidCut(kFALSE), fComparePtHardAndJetPt(kFALSE), fPtHardAndJetPtFactor(7),
56     fCTSPtMin(0), fEMCALPtMin(0),fPHOSPtMin(0), fAODBranchList(new TList ),
57     fAODCTS(new TObjArray()), fAODEMCAL(new TObjArray()), fAODPHOS(new TObjArray()),
58     fEMCALCells(0x0), fPHOSCells(0x0),
59     fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
60     fFillCTS(0),fFillEMCAL(0),fFillPHOS(0),
61     fFillEMCALCells(0),fFillPHOSCells(0), 
62 //    fSecondInputAODTree(0x0), fSecondInputAODEvent(0x0),
63 //    fSecondInputFileName(""),fSecondInputFirstEvent(0), 
64 //    fAODCTSNormalInputEntries(0), fAODEMCALNormalInputEntries(0), 
65 //    fAODPHOSNormalInputEntries(0), 
66     fTrackStatus(0),   fESDtrackCuts(0), fTrackMult(0), fTrackMultEtaCut(0.8),
67     fReadStack(kFALSE), fReadAODMCParticles(kFALSE), 
68     fDeltaAODFileName("deltaAODPartCorr.root"),fFiredTriggerClassName(""),
69     fAnaLED(kFALSE),fTaskName(""),fCaloUtils(0x0), 
70     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), 
71     fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),fCaloFilterPatch(kFALSE),
72     fEMCALClustersListName("")
73 {
74   //Ctor
75   
76   //Initialize parameters
77   InitParameters();
78 }
79
80 //_________________________________
81 AliCaloTrackReader::~AliCaloTrackReader() {
82   //Dtor
83   
84   if(fFiducialCut) delete fFiducialCut ;
85         
86   if(fAODBranchList){
87     fAODBranchList->Delete();
88     delete fAODBranchList ;
89   }  
90   
91   if(fAODCTS){
92     if(fDataType!=kMC)fAODCTS->Clear() ; 
93     else              fAODCTS->Delete() ; 
94     delete fAODCTS ;
95   }
96   
97   if(fAODEMCAL){
98     if(fDataType!=kMC)fAODEMCAL->Clear("C") ; 
99     else              fAODEMCAL->Delete() ; 
100     delete fAODEMCAL ;
101   }
102   
103   if(fAODPHOS){
104     if(fDataType!=kMC)fAODPHOS->Clear("C") ; 
105     else              fAODPHOS->Delete() ; 
106     delete fAODPHOS ;
107   }
108   
109 //  if(fEMCALCells){
110 //    delete fEMCALCells ;
111 //  }
112 //  
113 //  if(fPHOSCells){
114 //    delete fPHOSCells ;
115 //  }
116
117   if(fVertex){
118     for (Int_t i = 0; i < fNMixedEvent; i++) {
119       delete [] fVertex[i] ;
120
121     }
122     delete [] fVertex ;
123         }
124
125   if(fESDtrackCuts)   delete fESDtrackCuts;
126   
127 //  Pointers not owned, done by the analysis frame
128 //  if(fInputEvent)  delete fInputEvent ;
129 //  if(fOutputEvent) delete fOutputEvent ;
130 //  if(fMC)          delete fMC ;  
131         
132 //  if(fSecondInputAODTree){
133 //    fSecondInputAODTree->Clear();
134 //    delete fSecondInputAODTree;
135 //  }
136 //      
137 //  if(fSecondInputAODEvent) delete fSecondInputAODEvent ;
138         
139   //  Pointer not owned, deleted by maker
140   //if (fCaloUtils) delete fCaloUtils ;
141
142 }
143
144 //_________________________________________________________________________
145 Bool_t AliCaloTrackReader::ComparePtHardAndJetPt(){
146   // Check the event, if the requested ptHard is much larger than the jet pT, then there is a problem.
147   // Only for PYTHIA.
148   if(!fReadStack) return kTRUE; //Information not filtered to AOD
149   
150   if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader")){
151     TParticle * jet =  0;
152     AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
153     Int_t nTriggerJets =  pygeh->NTriggerJets();
154     Float_t ptHard = pygeh->GetPtHard();
155     
156     //if(fDebug > 1) printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %f\n",nTriggerJets, ptHard);
157     Float_t tmpjet[]={0,0,0,0};
158     for(Int_t ijet = 0; ijet< nTriggerJets; ijet++){
159       pygeh->TriggerJet(ijet, tmpjet);
160       jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
161       //Compare jet pT and pt Hard
162       //if(fDebug > 1) printf("AliMCAnalysisUtils:: %d pycell jet pT %f\n",ijet, jet->Pt());
163       if(jet->Pt() > fPtHardAndJetPtFactor * ptHard) {
164         printf("AliMCAnalysisUtils::PythiaEventHeader: Njets: %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
165                nTriggerJets, ptHard, jet->Pt(), fPtHardAndJetPtFactor);
166         return kFALSE;
167       }
168     }
169     if(jet) delete jet; 
170   }
171   
172   return kTRUE ;
173   
174 }
175
176 //____________________________________________________________________________
177 AliStack* AliCaloTrackReader::GetStack() const {
178   //Return pointer to stack
179   if(fMC)
180     return fMC->Stack();
181   else{
182     if(fDebug > 1) printf("AliCaloTrackReader::GetStack() - Stack is not available\n"); 
183     return 0x0 ;
184   }
185 }
186
187 //____________________________________________________________________________
188 AliHeader* AliCaloTrackReader::GetHeader() const {
189   //Return pointer to header
190   if(fMC)
191     return fMC->Header();
192   else{
193     printf("AliCaloTrackReader::Header is not available\n"); 
194     return 0x0 ;
195   }
196 }
197 //____________________________________________________________________________
198 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const {
199   //Return pointer to Generated event header
200   if(fMC)
201     return fMC->GenEventHeader();
202   else{
203     printf("AliCaloTrackReader::GenEventHeader is not available\n"); 
204     return 0x0 ;
205   }
206 }
207
208 //____________________________________________________________________________
209 TClonesArray* AliCaloTrackReader::GetAODMCParticles(Int_t input) const {
210   //Return list of particles in AOD. Do it for the corresponding input event.
211   
212   TClonesArray * rv = NULL ; 
213   if(fDataType == kAOD){
214  
215     if(input == 0){
216       //Normal input AOD
217       AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
218       if(evt)
219         rv = (TClonesArray*)evt->FindListObject("mcparticles");
220       else  
221         printf("AliCaloTrackReader::GetAODMCParticles() - wrong AOD input index? %d, or non existing tree? \n",input); 
222       
223     } //else if(input == 1 && fSecondInputAODEvent){ //Second input AOD   
224 //      rv = (TClonesArray*) fSecondInputAODEvent->FindListObject("mcparticles");       
225 //    } 
226     
227   } else {
228       printf("AliCaloTrackReader::GetAODMCParticles() - Input are not AODs\n"); 
229   }
230   
231   return rv ; 
232 }
233
234 //____________________________________________________________________________
235 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader(Int_t input) const {
236   //Return MC header in AOD. Do it for the corresponding input event.
237   AliAODMCHeader *mch = NULL;
238   if(fDataType == kAOD){
239     //Normal input AOD
240     if(input == 0) {
241       mch = (AliAODMCHeader*)((AliAODEvent*)fInputEvent)->FindListObject("mcheader");
242     }
243     //          //Second input AOD
244     //          else if(input == 1){ 
245     //       mch = (AliAODMCHeader*) fSecondInputAODEvent->FindListObject("mcheader");
246     //  }
247     else {
248       printf("AliCaloTrackReader::GetAODMCHeader() - wrong AOD input index, %d\n",input);
249     }
250   }
251   else {
252     printf("AliCaloTrackReader::GetAODMCHeader() - Input are not AODs\n");
253   }
254   
255   return mch;
256 }
257
258 //_______________________________________________________________
259 void AliCaloTrackReader::Init()
260 {
261   //Init reader. Method to be called in AliAnaPartCorrMaker
262   
263   //Get the file with second input events if the filename is given
264   //Get the tree and connect the AODEvent. Only with AODs
265   
266   if(fReadStack && fReadAODMCParticles){
267     printf("AliCaloTrackReader::Init() - Cannot access stack and mcparticles at the same time, change them \n");
268     fReadStack = kFALSE;
269     fReadAODMCParticles = kFALSE;
270   }
271   
272 //      if(fSecondInputFileName!=""){
273 //              if(fDataType == kAOD){
274 //                      TFile * input2   = new TFile(fSecondInputFileName,"read");
275 //                      printf("AliCaloTrackReader::Init() - Second input file opened: %s, size %d \n", input2->GetName(), (Int_t) input2->GetSize());
276 //                      fSecondInputAODTree = (TTree*) input2->Get("aodTree");
277 //                      if(fSecondInputAODTree) printf("AliCaloTrackReader::Init() - Second input tree opened: %s, entries %d \n", 
278 //                                                                                 fSecondInputAODTree->GetName(), (Int_t) fSecondInputAODTree->GetEntries());
279 //                      else{
280 //                       printf("AliCaloTrackReader::Init() - Second input tree not available, STOP \n");
281 //                       abort();
282 //                      }
283 //                      fSecondInputAODEvent = new AliAODEvent;
284 //                      fSecondInputAODEvent->ReadFromTree(fSecondInputAODTree);
285 //                      if(fSecondInputFirstEvent >= fSecondInputAODTree->GetEntriesFast()){
286 //                              printf("AliCaloTrackReader::Init() - Requested first event of second input %d, is larger than number of events %d, STOP\n", 
287 //                                         fSecondInputFirstEvent, (Int_t) fSecondInputAODTree->GetEntriesFast());
288 //                              abort();
289 //                      }
290 //              }
291 //              else printf("AliCaloTrackReader::Init() - Second input not added, reader is not AOD\n");
292 //      }
293 }
294
295 //_______________________________________________________________
296 void AliCaloTrackReader::InitParameters()
297 {
298   //Initialize the parameters of the analysis.
299   fDataType   = kESD ;
300   fCTSPtMin   = 0.2 ;
301   fEMCALPtMin = 0.2 ;
302   fPHOSPtMin  = 0.2 ;
303
304   //Do not filter the detectors input by default.
305   fFillEMCAL      = kFALSE;
306   fFillPHOS       = kFALSE;
307   fFillCTS        = kFALSE;
308   fFillEMCALCells = kFALSE;
309   fFillPHOSCells  = kFALSE;
310
311   //fSecondInputFileName   = "" ;
312   //fSecondInputFirstEvent = 0 ;
313   fReadStack             = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
314   fReadAODMCParticles    = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
315   fDeltaAODFileName      = "deltaAODPartCorr.root";
316   fFiredTriggerClassName      = "";
317                 
318   fAnaLED = kFALSE;
319   
320   //We want tracks fitted in the detectors:
321   //fTrackStatus=AliESDtrack::kTPCrefit;
322   //fTrackStatus|=AliESDtrack::kITSrefit;  
323   
324   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
325
326   fV0ADC[0] = 0;   fV0ADC[1] = 0; 
327   fV0Mul[0] = 0;   fV0Mul[1] = 0; 
328
329 }
330
331 //________________________________________________________________
332 void AliCaloTrackReader::Print(const Option_t * opt) const
333 {
334
335   //Print some relevant parameters set for the analysis
336   if(! opt)
337     return;
338
339   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
340   printf("Task name      : %s\n", fTaskName.Data()) ;
341   printf("Data type      : %d\n", fDataType) ;
342   printf("CTS Min pT     : %2.1f GeV/c\n", fCTSPtMin) ;
343   printf("EMCAL Min pT   : %2.1f GeV/c\n", fEMCALPtMin) ;
344   printf("PHOS Min pT    : %2.1f GeV/c\n", fPHOSPtMin) ;
345   printf("Use CTS         =     %d\n", fFillCTS) ;
346   printf("Use EMCAL       =     %d\n", fFillEMCAL) ;
347   printf("Use PHOS        =     %d\n", fFillPHOS) ;
348   printf("Use EMCAL Cells =     %d\n", fFillEMCALCells) ;
349   printf("Use PHOS  Cells =     %d\n", fFillPHOSCells) ;
350   printf("Track status    =     %d\n", (Int_t) fTrackStatus) ;
351   printf("Track Mult Eta Cut =  %d\n", (Int_t) fTrackMultEtaCut) ;
352   printf("Write delta AOD =     %d\n", fWriteOutputDeltaAOD) ;
353
354   if(fComparePtHardAndJetPt)
355           printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
356         
357 //  if(fSecondInputFileName!="") {
358 //        printf("Second Input File Name     =     %s\n", fSecondInputFileName.Data()) ;
359 //        printf("Second Input First Event   =     %d\n", fSecondInputFirstEvent) ;
360 //  }
361         
362   printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
363   printf("Delta AOD File Name =     %s\n", fDeltaAODFileName.Data()) ;
364   printf("    \n") ;
365
366
367 //___________________________________________________
368 Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) {
369   //Fill the event counter and input lists that are needed, called by the analysis maker.
370
371   fEventNumber = iEntry;
372   fCurrentFileName = TString(currentFileName);
373   if(!fInputEvent) {
374           if(fDebug >= 0) printf("AliCaloTrackReader::FillInputEvent() - Input event not available, skip event analysis\n");
375           return kFALSE;
376   }
377   //Select events only fired by a certain trigger configuration if it is provided
378   Int_t eventType = 0;
379   if(fInputEvent->GetHeader())
380           eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
381   if( fFiredTriggerClassName  !="" && !fAnaLED){
382     if(eventType!=7)
383       return kFALSE; //Only physics event, do not use for simulated events!!!
384     if(fDebug > 0) 
385       printf("AliCaloTrackReader::FillInputEvent() - FiredTriggerClass <%s>, selected class <%s>, compare name %d\n",
386              GetFiredTriggerClasses().Data(),fFiredTriggerClassName.Data(), GetFiredTriggerClasses().Contains(fFiredTriggerClassName));
387     if( !GetFiredTriggerClasses().Contains(fFiredTriggerClassName) ) return kFALSE;
388   }
389   else if(fAnaLED){
390 //        kStartOfRun =       1,    // START_OF_RUN
391 //        kEndOfRun =         2,    // END_OF_RUN
392 //        kStartOfRunFiles =  3,    // START_OF_RUN_FILES
393 //        kEndOfRunFiles =    4,    // END_OF_RUN_FILES
394 //        kStartOfBurst =     5,    // START_OF_BURST
395 //        kEndOfBurst =       6,    // END_OF_BURST
396 //        kPhysicsEvent =     7,    // PHYSICS_EVENT
397 //        kCalibrationEvent = 8,    // CALIBRATION_EVENT
398 //        kFormatError =      9,    // EVENT_FORMAT_ERROR
399 //        kStartOfData =      10,   // START_OF_DATA
400 //        kEndOfData =        11,   // END_OF_DATA
401 //        kSystemSoftwareTriggerEvent   = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
402 //        kDetectorSoftwareTriggerEvent = 13  // DETECTOR_SOFTWARE_TRIGGER_EVENT
403          
404           if(eventType!=7 && fDebug > 1 )printf("AliCaloTrackReader::FillInputEvent() - DO LED, Event Type <%d>, 8 Calibration \n",  eventType);
405           if(eventType!=8)return kFALSE;
406   }
407                 
408   //In case of analysis of events with jets, skip those with jet pt > 5 pt hard 
409   if(fComparePtHardAndJetPt && GetStack()) {
410     if(!ComparePtHardAndJetPt()) return kFALSE ;
411   }
412
413   //In case of mixing events with other AOD file        
414  // if(fDataType == kAOD && fSecondInputAODTree){
415 //       
416 //    if(fDebug > 1) 
417 //      printf("AliCaloTrackReader::FillInputEvent() - Get event %d from second input AOD file \n", iEntry+fSecondInputFirstEvent);
418 //    if(fSecondInputAODTree->GetEntriesFast() <= iEntry+fSecondInputFirstEvent) {
419 //      if(fSecondInputAODTree->GetEntriesFast() == iEntry+fSecondInputFirstEvent) 
420 //                       printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry);
421 //      return kFALSE;
422 //    }
423 //    
424 //    //Get the Event
425 //    Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent);
426 //    if ( nbytes == 0 ) {//If nothing in AOD
427 //      printf("AliCaloTrackReader::FillInputEvent() - Nothing in Second AOD input, STOP\n");
428 //      abort() ; 
429 //    }
430 //    
431 //  }
432         
433   //Fill Vertex array
434   
435   FillVertexArray();
436   
437   //Fill the arrays with cluster/tracks/cells data
438    if(fFillEMCALCells) 
439     FillInputEMCALCells();
440   if(fFillPHOSCells)  
441     FillInputPHOSCells();
442         
443   if(fFillCTS){   
444     FillInputCTS();
445     //Accept events with at least one track
446     if(fTrackMult == 0) return kFALSE;
447   }
448   
449   //In case of data produced with calo filter, some information stored in non usual places
450   if(IsCaloFilterPatchOn()){
451     fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
452     //printf("Track multiplicity %d \n",fTrackMult);
453     if(fTrackMult == 0) return kFALSE;
454   }
455   
456   if(fFillEMCAL) 
457     FillInputEMCAL();
458   if(fFillPHOS)  
459     FillInputPHOS();
460
461   FillInputVZERO();
462         
463   return kTRUE ;
464 }
465
466 //__________________________________________________
467 void AliCaloTrackReader::ResetLists() {
468   //  Reset lists, called by the analysis maker 
469
470   if(fAODCTS)     fAODCTS     -> Clear();
471   if(fAODEMCAL)   fAODEMCAL   -> Clear("C");
472   if(fAODPHOS)    fAODPHOS    -> Clear("C");
473 //  if(fEMCALCells) fEMCALCells -> Clear("");
474 //  if(fPHOSCells)  fPHOSCells  -> Clear("");
475
476   fV0ADC[0] = 0;   fV0ADC[1] = 0; 
477   fV0Mul[0] = 0;   fV0Mul[1] = 0; 
478
479 }
480
481 //____________________________________________________________________________
482 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)  
483 {
484   fInputEvent  = input;
485   fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ; 
486   if (fMixedEvent) {
487     fNMixedEvent = fMixedEvent->GetNumberOfEvents() ; 
488   }
489
490   //Delete previous vertex
491   if(fVertex){
492     for (Int_t i = 0; i < fNMixedEvent; i++) {
493       delete [] fVertex[i] ; 
494     }
495     delete [] fVertex ;
496   }
497   
498   fVertex = new Double_t*[fNMixedEvent] ; 
499   for (Int_t i = 0; i < fNMixedEvent; i++) {
500     fVertex[i] = new Double_t[3] ; 
501     fVertex[i][0] = 0.0 ; 
502     fVertex[i][1] = 0.0 ; 
503     fVertex[i][2] = 0.0 ; 
504   }
505 }
506
507 //____________________________________________________________________________
508 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const {
509   //Return vertex position to be used for single event analysis
510   vertex[0]=fVertex[0][0];  
511   vertex[1]=fVertex[0][1];  
512   vertex[2]=fVertex[0][2];
513 }
514
515 //____________________________________________________________________________
516 void AliCaloTrackReader::GetVertex(Double_t vertex[3], const Int_t evtIndex) const {
517   //Return vertex position for mixed event, recover the vertex in a particular event.
518   
519   //Int_t evtIndex = 0; // for single events only one vertex stored in position 0, default value
520   //if (fMixedEvent && clusterID >=0) {
521   //  evtIndex=GetMixedEvent()->EventIndexForCaloCluster(clusterID) ; 
522   //}
523   
524   vertex[0]=fVertex[evtIndex][0];  vertex[1]=fVertex[evtIndex][1];  vertex[2]=fVertex[evtIndex][2];
525   
526 }
527 //
528
529
530 //____________________________________________________________________________
531 void AliCaloTrackReader::FillVertexArray() {
532   
533   //Fill data member with vertex
534   //In case of Mixed event, multiple vertices
535   
536   //Delete previous vertex
537   if(fVertex){
538     for (Int_t i = 0; i < fNMixedEvent; i++) {
539       delete [] fVertex[i] ; 
540     }
541     delete [] fVertex ;  
542   }
543   
544   fVertex = new Double_t*[fNMixedEvent] ; 
545   for (Int_t i = 0; i < fNMixedEvent; i++) {
546     fVertex[i] = new Double_t[3] ; 
547     fVertex[i][0] = 0.0 ; 
548     fVertex[i][1] = 0.0 ; 
549     fVertex[i][2] = 0.0 ; 
550   }          
551   
552   if (!fMixedEvent) { //Single event analysis
553     if(fDataType!=kMC){
554
555       if(fInputEvent->GetPrimaryVertex()){
556         fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]); 
557       }
558       else {
559         printf("AliCaloTrackReader::FillVertexArray() - NULL primary vertex\n");
560         fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
561       }//Primary vertex pointer do not exist
562       
563     } else {//MC read event 
564       fVertex[0][0]=0.;   fVertex[0][1]=0.;   fVertex[0][2]=0.;
565     }
566       
567     if(fDebug > 1)
568       printf("AliCaloTrackReader::FillVertexArray() - Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]);
569
570   } else { // MultiEvent analysis
571     for (Int_t iev = 0; iev < fNMixedEvent; iev++) {
572       if (fMixedEvent->GetVertexOfEvent(iev))
573         fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
574       else { // no vertex found !!!!
575         AliWarning("No vertex found");
576       }
577
578       if(fDebug > 1)
579         printf("AliCaloTrackReader::FillVertexArray() - Multi Event %d Vertex : %f,%f,%f\n",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]);
580
581     }
582   }
583   
584 }
585
586 //____________________________________________________________________________
587 void AliCaloTrackReader::FillInputCTS() {
588   //Return array with Central Tracking System (CTS) tracks
589   
590   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputCTS()\n");
591   
592   Int_t nTracks   = fInputEvent->GetNumberOfTracks() ;
593   Double_t p[3];
594   fTrackMult = 0;
595   Int_t nstatus = 0;
596   for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
597     AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
598
599     //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
600     if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) 
601       continue ;
602     
603     nstatus++;
604     
605     if(fDataType==kESD && !fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
606     
607     // Track filter selection
608     //if (fTrackFilter) {
609           //  selectInfo = fTrackFilter->IsSelected(esdTrack);
610           //  if (!selectInfo && !(esd->GetPrimaryVertex())->UsesTrack(esdTrack->GetID())) continue;
611    // }
612     
613     //Count the tracks in eta < 0.9
614     //printf("Eta %f cut  %f\n",TMath::Abs(track->Eta()),fTrackMultEtaCut);
615     if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
616     
617     track->GetPxPyPz(p) ;
618     TLorentzVector momentum(p[0],p[1],p[2],0);
619     
620     if(fCTSPtMin < momentum.Pt()){
621       
622       if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) 
623         continue;
624       
625       if(fDebug > 2 && momentum.Pt() > 0.1) 
626         printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
627                momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
628       
629       if (fMixedEvent) {
630         track->SetID(itrack);
631       }
632       
633       fAODCTS->Add(track);        
634        
635     }//Pt and Fiducial cut passed. 
636   }// track loop
637         
638   //fAODCTSNormalInputEntries = fAODCTS->GetEntriesFast();
639   if(fDebug > 1) 
640     printf("AliCaloTrackReader::FillInputCTS()   - aod entries %d, input tracks %d, pass status %d, multipliticy %d\n", fAODCTS->GetEntriesFast(), nTracks, nstatus, fTrackMult);//fAODCTSNormalInputEntries);
641   
642     //  //If second input event available, add the clusters.
643     //  if(fSecondInputAODTree && fSecondInputAODEvent){
644     //    nTracks   = fSecondInputAODEvent->GetNumberOfTracks() ;
645     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS()   - Add second input tracks, entries %d\n", nTracks) ;
646     //    for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
647     //            AliAODTrack * track = ((AliAODEvent*)fSecondInputAODEvent)->GetTrack(itrack) ; // retrieve track from esd
648     //            
649     //            //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
650     //            if (fTrackStatus && !((track->GetStatus() & fTrackStatus) == fTrackStatus)) continue ;
651     //            
652     //            track->GetPxPyPz(p) ;
653     //            TLorentzVector momentum(p[0],p[1],p[2],0);
654     //            
655     //            if(fCTSPtMin < momentum.Pt()){
656     //
657     //                    if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"CTS")) continue;
658     //
659     //                    if(fDebug > 2 && momentum.Pt() > 0.1) printf("AliCaloTrackReader::FillInputCTS() - Selected tracks E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
660     //                                                                 momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
661     //                    
662     //                    fAODCTS->Add(track);
663     //                    
664     //            }//Pt and Fiducial cut passed. 
665     //    }// track loop
666     //    
667     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputCTS()   - aod normal entries %d, after second input %d\n", fAODCTSNormalInputEntries, fAODCTS->GetEntriesFast());
668     //  }       //second input loop
669     //  
670 }
671
672 //____________________________________________________________________________
673 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, const Int_t iclus) {
674   //Fill the EMCAL data in the array, do it
675   
676   Int_t vindex = 0 ;  
677   if (fMixedEvent) 
678     vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
679   
680   //Check if the cluster contains any bad channel and if close to calorimeter borders
681   if(GetCaloUtils()->ClusterContainsBadChannel("EMCAL",clus->GetCellsAbsId(), clus->GetNCells())) 
682     return;
683   if(!GetCaloUtils()->CheckCellFiducialRegion(clus, (AliVCaloCells*)fInputEvent->GetEMCALCells(), fInputEvent, vindex)) 
684     return;
685   
686   TLorentzVector momentum ;
687   
688   clus->GetMomentum(momentum, fVertex[vindex]);      
689   
690   if(fEMCALPtMin < momentum.Pt()){
691     
692     if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) 
693       return;
694     
695     if(fDebug > 2 && momentum.E() > 0.1) 
696       printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
697              momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
698     
699     //Float_t pos[3];
700     //clus->GetPosition(pos);
701     //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
702     
703     //Recalibrate the cluster energy 
704     if(GetCaloUtils()->IsRecalibrationOn()) {
705       Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, GetEMCALCells());
706       clus->SetE(energy);
707       //printf("Recalibrated Energy %f\n",clus->E());  
708       GetCaloUtils()->RecalculateClusterShowerShapeParameters(GetEMCALCells(),clus);
709       GetCaloUtils()->RecalculateClusterPID(clus);
710       
711     }
712     
713     //Recalculate distance to bad channels, if new list of bad channels provided
714     GetCaloUtils()->RecalculateClusterDistanceToBadChannel(GetEMCALCells(),clus);
715     
716     //Recalculate cluster position
717     if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
718       GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus); 
719       //clus->GetPosition(pos);
720       //printf("After  Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
721     }
722     
723     //Correct non linearity
724     if(GetCaloUtils()->IsCorrectionOfClusterEnergyOn()){
725       GetCaloUtils()->CorrectClusterEnergy(clus) ;
726       //printf("Linearity Corrected Energy %f\n",clus->E());  
727     }          
728     
729     if (fMixedEvent) 
730       clus->SetID(iclus) ; 
731     
732     fAODEMCAL->Add(clus);       
733   }
734 }
735
736 //____________________________________________________________________________
737 void AliCaloTrackReader::FillInputEMCAL() {
738   //Return array with EMCAL clusters in aod format
739   
740   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputEMCAL()\n");
741   
742   //Loop to select clusters in fiducial cut and fill container with aodClusters
743   if(fEMCALClustersListName==""){
744     Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
745     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
746       AliVCluster * clus = 0;
747       if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
748         if (IsEMCALCluster(clus)){          
749           FillInputEMCALAlgorithm(clus, iclus);
750         }//EMCAL cluster
751       }// cluster exists
752     }// cluster loop
753   }//Get the clusters from the input event
754   else {
755     TClonesArray * clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
756     if(!clusterList){
757       printf("AliCaloTrackReader::FillInputEMCAL() - Wrong name of list with clusters? <%s>\n",fEMCALClustersListName.Data());
758       return;
759     }
760     Int_t nclusters = clusterList->GetEntriesFast();
761     for (Int_t iclus =  0; iclus <  nclusters; iclus++) {
762       AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
763       //printf("E %f\n",clus->E());
764       if (clus) FillInputEMCALAlgorithm(clus, iclus);
765       else printf("AliCaloTrackReader::FillInputEMCAL() - Null cluster in list!\n");
766     }// cluster loop
767   }
768   
769   //Recalculate track matching
770   GetCaloUtils()->RecalculateClusterTrackMatching(fInputEvent);
771   
772   //fAODEMCALNormalInputEntries = fAODEMCAL->GetEntriesFast();
773   if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod entries %d\n",  fAODEMCAL->GetEntriesFast());//fAODEMCALNormalInputEntries);
774   
775     //If second input event available, add the clusters.
776     //  if(fSecondInputAODTree && fSecondInputAODEvent){
777     //    GetSecondInputAODVertex(v);
778     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
779     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - Add second input clusters, entries %d\n", nclusters) ;
780     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
781     //                  AliAODCaloCluster * clus = 0;
782     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
783     //                          if (clus->IsEMCAL()){
784     //                                  TLorentzVector momentum ;
785     //                                  clus->GetMomentum(momentum, v);      
786     //                                  
787     //                                  if(fEMCALPtMin < momentum.Pt()){
788     //
789     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"EMCAL")) continue;
790     //
791     //                                          if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
792     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
793     //                                    fAODEMCAL->Add(clus); 
794     //                                  }//Pt and Fiducial cut passed.
795     //                          }//EMCAL cluster
796     //                  }// cluster exists
797     //          }// cluster loop
798     //          
799     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputEMCAL() - aod normal entries %d, after second input %d\n", fAODEMCALNormalInputEntries, fAODEMCAL->GetEntriesFast());
800     //
801     //  } //second input loop
802 }
803
804   //____________________________________________________________________________
805 void AliCaloTrackReader::FillInputPHOS() {
806   //Return array with PHOS clusters in aod format
807   
808   if(fDebug > 2 ) printf("AliCaloTrackReader::FillInputPHOS()\n");
809           
810   //Loop to select clusters in fiducial cut and fill container with aodClusters
811   Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
812   for (Int_t iclus = 0; iclus < nclusters; iclus++) {
813     AliVCluster * clus = 0;
814     if ( (clus = fInputEvent->GetCaloCluster(iclus)) ) {
815       if (IsPHOSCluster(clus)){
816         //Check if the cluster contains any bad channel and if close to calorimeter borders
817         Int_t vindex = 0 ;  
818         if (fMixedEvent) 
819           vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
820         if( GetCaloUtils()->ClusterContainsBadChannel("PHOS",clus->GetCellsAbsId(), clus->GetNCells())) 
821           continue;
822         if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells(), fInputEvent, vindex)) 
823           continue;
824         
825         TLorentzVector momentum ;
826         
827         clus->GetMomentum(momentum, fVertex[vindex]);      
828         
829         if(fPHOSPtMin < momentum.Pt()){
830           
831           if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) 
832             continue;
833           
834           if(fDebug > 2 && momentum.E() > 0.1) 
835             printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
836                    momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
837           
838             //Recalibrate the cluster energy 
839           if(GetCaloUtils()->IsRecalibrationOn()) {
840             Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
841             clus->SetE(energy);
842           }
843           
844           if (fMixedEvent) {
845             clus->SetID(iclus) ; 
846           }              
847           
848           fAODPHOS->Add(clus);  
849           
850         }//Pt and Fiducial cut passed.
851       }//PHOS cluster
852     }//cluster exists
853   }//esd cluster loop
854   
855   //fAODPHOSNormalInputEntries = fAODPHOS->GetEntriesFast() ;
856   if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod entries %d\n",  fAODPHOS->GetEntriesFast());//fAODPHOSNormalInputEntries);
857   
858     //If second input event available, add the clusters.
859     //  if(fSecondInputAODTree && fSecondInputAODEvent){  
860     //    GetSecondInputAODVertex(v);
861     //    nclusters = ((AliAODEvent*)fSecondInputAODEvent)->GetNumberOfCaloClusters();
862     //    if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - Add second input clusters, entries %d\n", nclusters);
863     //          for (Int_t iclus =  0; iclus < nclusters; iclus++) {
864     //                  AliAODCaloCluster * clus = 0;
865     //                  if ( (clus = ((AliAODEvent*)fSecondInputAODEvent)->GetCaloCluster(iclus)) ) {
866     //                          if (clus->IsPHOS()){
867     //                                  TLorentzVector momentum ;
868     //                                  clus->GetMomentum(momentum, v);      
869     //                                  
870     //                                  if(fPHOSPtMin < momentum.Pt()){
871     //
872     //                                          if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(momentum,"PHOS")) continue;
873     //
874     //                                          if(fDebug > 2 && momentum.E() > 0.1) printf("AliCaloTrackReader::FillInputPHOS() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
875     //                                                                                                                                  momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
876     //                                          fAODPHOS->Add(clus);    
877     //                                  }//Pt and Fiducial cut passed.
878     //                          }//PHOS cluster
879     //                  }// cluster exists
880     //          }// cluster loop
881     //          if(fDebug > 1) printf("AliCaloTrackReader::FillInputPHOS()  - aod normal entries %d, after second input %d\n", fAODPHOSNormalInputEntries, fAODPHOS->GetEntriesFast());
882     //  }       //second input loop
883   
884 }
885
886 //____________________________________________________________________________
887 void AliCaloTrackReader::FillInputEMCALCells() {
888     //Return array with EMCAL cells in aod format
889   
890   fEMCALCells = fInputEvent->GetEMCALCells(); 
891   
892 }
893
894 //____________________________________________________________________________
895 void AliCaloTrackReader::FillInputPHOSCells() {
896     //Return array with PHOS cells in aod format
897   
898   fPHOSCells = fInputEvent->GetPHOSCells(); 
899   
900 }
901
902
903 //____________________________________________________________________________
904 Bool_t AliCaloTrackReader::IsEMCALCluster(AliVCluster* cluster) const {
905   // Check if it is a cluster from EMCAL. For old AODs cluster type has
906   // different number and need to patch here
907     
908   if(fDataType==kAOD && fOldAOD)
909   {
910     if (cluster->GetType() == 2) return kTRUE;
911     else                         return kFALSE;
912   }
913   else 
914   {
915     return cluster->IsEMCAL();
916   }
917
918 }
919
920 //____________________________________________________________________________
921 Bool_t AliCaloTrackReader::IsPHOSCluster(AliVCluster * cluster) const {
922   //Check if it is a cluster from PHOS.For old AODs cluster type has
923   // different number and need to patch here
924   
925   if(fDataType==kAOD && fOldAOD)
926   {
927     Int_t type = cluster->GetType();
928     if (type == 0 || type == 1) return kTRUE;
929     else                        return kFALSE;
930   }
931   else 
932   {
933     return cluster->IsPHOS();
934   }
935   
936 }
937