1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-------------------------------------------------------------------------
18 // ESD reader for jet analysis
19 // Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
20 // Magali Estienne <magali.estienne@subatech.in2p3.fr>
21 //-------------------------------------------------------------------------
23 // --- Standard library ---
24 #include <Riostream.h>
26 // --- ROOT system ---
28 #include <TStopwatch.h>
29 #include <TLorentzVector.h>
34 #include <TGeoManager.h>
36 #include <TRefArray.h>
38 #include <TProcessID.h>
41 // --- AliRoot header files ---
42 #include "AliJetESDReader.h"
43 #include "AliJetESDReaderHeader.h"
44 #include "AliESDEvent.h"
46 #include "AliESDtrack.h"
47 #include "AliJetDummyGeo.h"
48 #include "AliJetESDFillUnitArrayTracks.h"
49 #include "AliJetESDFillUnitArrayEMCalDigits.h"
50 #include "AliJetHadronCorrectionv1.h"
51 #include "AliJetUnitArray.h"
52 #include "AliAnalysisTask.h"
54 ClassImp(AliJetESDReader)
56 AliJetESDReader::AliJetESDReader():
67 fApplyElectronCorrection(kFALSE),
68 fApplyMIPCorrection(kTRUE),
69 fApplyFractionHadronicCorrection(kFALSE),
70 fFractionHadronicCorrection(0.3),
86 //____________________________________________________________________________
87 AliJetESDReader::~AliJetESDReader()
105 //____________________________________________________________________________
106 void AliJetESDReader::OpenInputFiles()
108 // chain for the ESDs
109 fChain = new TChain("esdTree");
111 // get directory and pattern name from the header
112 const char* dirName=fReaderHeader->GetDirectory();
113 const char* pattern=fReaderHeader->GetPattern();
115 // Add files matching patters to the chain
116 void *dir = gSystem->OpenDirectory(dirName);
117 const char *name = 0x0;
118 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
120 while ((name = gSystem->GetDirEntry(dir))){
121 if (a>=nesd) continue;
123 if (strstr(name,pattern)){
124 printf("Adding %s\n",name);
125 fChain->AddFile(Form("%s/%s/AliESDs.root",dirName,name));
130 gSystem->FreeDirectory(dir);
132 fESD = new AliESDEvent();
133 fESD->ReadFromTree(fChain);
135 int nMax = fChain->GetEntries();
137 printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
139 // set number of events in header
140 if (fReaderHeader->GetLastEvent() == -1)
141 fReaderHeader->SetLastEvent(nMax);
143 Int_t nUsr = fReaderHeader->GetLastEvent();
144 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
149 //__________________________________________________________
150 void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
153 // Set flag to apply MIP correction fApplyMIPCorrection
154 // - exclusive with fApplyFractionHadronicCorrection
157 fApplyMIPCorrection = val;
158 if(fApplyMIPCorrection == kTRUE)
160 SetApplyFractionHadronicCorrection(kFALSE);
161 printf("Enabling MIP Correction \n");
165 printf("Disabling MIP Correction \n");
169 //__________________________________________________________
170 void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
173 // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
174 // - exclusive with fApplyMIPCorrection
177 fApplyFractionHadronicCorrection = val;
178 if(fApplyFractionHadronicCorrection == kTRUE)
180 SetApplyMIPCorrection(kFALSE);
181 printf("Enabling Fraction Hadronic Correction \n");
185 printf("Disabling Fraction Hadronic Correction \n");
189 //__________________________________________________________
190 void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
193 // Set value to fFractionHadronicCorrection (default is 0.3)
194 // apply EMC hadronic correction fApplyFractionHadronicCorrection
195 // - exclusive with fApplyMIPCorrection
198 fFractionHadronicCorrection = val;
199 if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
201 SetApplyFractionHadronicCorrection(kTRUE);
202 printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
206 SetApplyFractionHadronicCorrection(kFALSE);
211 //____________________________________________________________________________
212 void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /*aod*/, const TObject* /*mc*/) {
214 fESD = (AliESDEvent*) esd;
217 //____________________________________________________________________________
218 Bool_t AliJetESDReader::FillMomentumArray()
220 // Fill momentum array
227 // clear momentum array
229 fDebug = fReaderHeader->GetDebug();
230 fOpt = fReaderHeader->GetDetector();
236 // get number of tracks in event (for the loop)
237 nt = fESD->GetNumberOfTracks();
238 printf("Fill Momentum Array %5d \n", nt);
240 // temporary storage of signal and pt cut flag
241 Int_t* sflag = new Int_t[nt];
242 Int_t* cflag = new Int_t[nt];
244 // get cuts set by user
245 Float_t ptMin = fReaderHeader->GetPtCut();
246 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
247 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
249 //loop over tracks in ESD
250 for (Int_t it = 0; it < nt; it++) {
251 AliESDtrack *track = fESD->GetTrack(it);
252 UInt_t status = track->GetStatus();
255 track->GetPxPyPz(mom);
256 p3.SetXYZ(mom[0],mom[1],mom[2]);
258 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
259 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
260 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
261 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
262 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
263 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
265 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
267 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
269 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
271 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
275 // set the signal flags
276 fSignalFlag.Set(goodTrack,sflag);
277 fCutFlag.Set(goodTrack,cflag);
286 //____________________________________________________________________________
287 void AliJetESDReader::CreateTasks(TChain* tree)
290 // For reader task initialization
293 fDebug = fReaderHeader->GetDebug();
294 fDZ = fReaderHeader->GetDZ();
297 // Init EMCAL geometry
301 // Create and init unit array
304 // Create global reader task for analysis
305 fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
306 // Create a task for to fill the charged particle information
307 fFillUAFromTracks = new AliJetESDFillUnitArrayTracks();
308 fFillUAFromTracks->SetReaderHeader(fReaderHeader);
309 fFillUAFromTracks->SetGeom(fGeom);
310 fFillUAFromTracks->SetTPCGrid(fTpcGrid);
311 fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
313 { // Calo dead zones inclusion
314 fFillUAFromTracks->SetGrid0(fGrid0);
315 fFillUAFromTracks->SetGrid1(fGrid1);
316 fFillUAFromTracks->SetGrid2(fGrid2);
317 fFillUAFromTracks->SetGrid3(fGrid3);
318 fFillUAFromTracks->SetGrid4(fGrid4);
320 fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
321 fFillUAFromTracks->SetHadCorrector(fHadCorr);
322 // Create a task for to fill the neutral particle information
323 fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
324 fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
325 fFillUAFromEMCalDigits->SetGeom(fGeom);
326 fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
327 fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
328 fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
329 fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
330 fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
331 // Add the task to global task
332 fFillUnitArray->Add(fFillUAFromTracks);
333 fFillUnitArray->Add(fFillUAFromEMCalDigits);
334 fFillUAFromTracks->SetActive(kFALSE);
335 fFillUAFromEMCalDigits->SetActive(kFALSE);
337 cout << "Tasks instantiated at that stage ! " << endl;
338 cout << "You can loop over events now ! " << endl;
342 //____________________________________________________________________________
343 Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
346 // Reader task execussion
350 fRefArray = refArray;
352 // clear momentum array
355 fDebug = fReaderHeader->GetDebug();
356 fOpt = fReaderHeader->GetDetector();
362 // TPC only or Digits+TPC or Clusters+TPC
363 if(fOpt%2==!0 && fOpt!=0){
364 fFillUAFromTracks->SetESD(fESD);
365 fFillUAFromTracks->SetActive(kTRUE);
366 fFillUAFromTracks->SetUnitArray(fUnitArray);
367 fFillUAFromTracks->SetRefArray(fRefArray);
368 fFillUAFromTracks->SetProcId(fProcId);
369 // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
370 fFillUAFromTracks->Exec("tpc"); // Temporarily added
372 fNumCandidate = fFillUAFromTracks->GetMult();
373 fNumCandidateCut = fFillUAFromTracks->GetMultCut();
377 // Digits only or Digits+TPC
378 if(fOpt>=2 && fOpt<=3){
379 fFillUAFromEMCalDigits->SetESD(fESD);
380 fFillUAFromEMCalDigits->SetActive(kTRUE);
381 fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
382 fFillUAFromEMCalDigits->SetRefArray(fRefArray);
383 fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
384 fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
385 fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
386 fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added
387 fNumCandidate = fFillUAFromEMCalDigits->GetMult();
388 fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
391 // fFillUnitArray->ExecuteTask(); // Temporarily commented
396 //____________________________________________________________________________
397 Bool_t AliJetESDReader::SetEMCALGeometry()
400 // Set the EMCal Geometry
403 if (!fTree->GetFile())
406 TString geomFile(fTree->GetFile()->GetName());
407 geomFile.ReplaceAll("AliESDs", "geometry");
409 // temporary workaround for PROOF bug #18505
410 geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
411 if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
413 // Define EMCAL geometry to be able to read ESDs
414 fGeom = AliJetDummyGeo::GetInstance();
416 fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
418 // To be setted to run some AliEMCALGeometry functions
419 TGeoManager::Import(geomFile);
420 fGeom->GetTransformationForSM();
421 printf("\n EMCal Geometry set ! \n");
426 //____________________________________________________________________________
427 void AliJetESDReader::InitParameters()
429 // Initialise parameters
430 fOpt = fReaderHeader->GetDetector();
431 // fHadCorr = 0; // For hadron correction
433 if(fOpt==0 || fOpt==1)
434 fECorrection = 0; // For electron correction
435 else fECorrection = 1; // For electron correction
437 fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
438 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
441 //____________________________________________________________________________
442 void AliJetESDReader::InitUnitArray()
445 // Create and Initialise unit arrays
448 Int_t nElements = fTpcGrid->GetNEntries();
449 Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
450 if(fArrayInitialised) fUnitArray->Delete();
452 if(fTpcGrid->GetGridType()==0)
453 { // Fill the following quantities :
454 // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
455 // detector flag, in/out jet, pt cut, mass, cluster ID)
456 for(Int_t nBin = 1; nBin < nElements+1; nBin++)
458 fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
459 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
460 Deta = fTpcGrid->GetDeta();
461 Dphi = fTpcGrid->GetDphi();
462 new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
466 if(fTpcGrid->GetGridType()==1)
469 Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
473 // Define a grid of cell for the gaps between SM
474 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
475 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
476 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
477 fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
478 fGrid0->SetGridType(0);
479 fGrid0->SetMatrixIndexes();
480 fGrid0->SetIndexIJ();
481 n0 = fGrid0->GetNEntries();
482 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
483 fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
484 fGrid1->SetGridType(0);
485 fGrid1->SetMatrixIndexes();
486 fGrid1->SetIndexIJ();
487 n1 = fGrid1->GetNEntries();
488 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
489 fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
490 fGrid2->SetGridType(0);
491 fGrid2->SetMatrixIndexes();
492 fGrid2->SetIndexIJ();
493 n2 = fGrid2->GetNEntries();
494 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
495 fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
496 fGrid3->SetGridType(0);
497 fGrid3->SetMatrixIndexes();
498 fGrid3->SetIndexIJ();
499 n3 = fGrid3->GetNEntries();
500 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
501 fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
502 fGrid4->SetGridType(0);
503 fGrid4->SetMatrixIndexes();
504 fGrid4->SetIndexIJ();
505 n4 = fGrid4->GetNEntries();
507 nGaps = n0+n1+n2+n3+n4;
511 for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
515 fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
516 // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
517 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
518 Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
519 Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
520 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
523 if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
524 fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
525 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
526 Deta = fTpcGrid->GetDeta();
527 Dphi = fTpcGrid->GetDphi();
528 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
532 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
533 if(nBin<fNumUnits+nElements+n0)
535 Float_t phi0 = eta = 0.;
536 fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
537 Deta = fGrid0->GetDeta();
538 Dphi = fGrid0->GetDphi();
539 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
541 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
543 Float_t phi0 = eta = 0.;
544 fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta);
545 Deta = fGrid1->GetDeta();
546 Dphi = fGrid1->GetDphi();
547 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
549 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
551 Float_t phi0 = eta = 0.;
552 fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
553 Deta = fGrid2->GetDeta();
554 Dphi = fGrid2->GetDphi();
555 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
557 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
559 Float_t phi0 = eta = 0.;
560 fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
561 Deta = fGrid3->GetDeta();
562 Dphi = fGrid3->GetDphi();
563 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
565 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
567 Float_t phi0 = eta = 0.;
568 fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
569 Deta = fGrid4->GetDeta();
570 Dphi = fGrid4->GetDphi();
571 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
577 } // end loop on nBin
578 } // end grid type == 1
579 fArrayInitialised = 1;