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