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