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>
23 // implemented standard geometry (AliEMCALGeometry) instead of dummy one (AliJetDummyGeo)
24 // moved geometry definition in AliJetReader
25 // marco.bregant@subatech.in2p3.fr
26 //-------------------------------------------------------------------------
28 // --- Standard library ---
29 #include <Riostream.h>
31 // --- ROOT system ---
34 #include <TStopwatch.h>
35 #include <TLorentzVector.h>
40 #include <TGeoManager.h>
42 #include <TRefArray.h>
44 #include <TProcessID.h>
47 // --- AliRoot header files ---
48 #include "AliJetESDReader.h"
49 #include "AliJetESDReaderHeader.h"
50 #include "AliESDEvent.h"
51 #include "AliVEvent.h"
53 #include "AliESDtrack.h"
54 #include "AliEMCALGeometry.h"
55 #include "AliEMCALEMCGeometry.h"
56 #include "AliJetESDFillUnitArrayTracks.h"
57 #include "AliJetESDFillUnitArrayEMCalDigits.h"
58 #include "AliJetHadronCorrectionv1.h"
59 #include "AliJetUnitArray.h"
60 #include "AliAnalysisTask.h"
61 #include "AliOADBContainer.h"
63 ClassImp(AliJetESDReader)
65 AliJetESDReader::AliJetESDReader():
75 fApplyElectronCorrection(kFALSE),
76 fApplyMIPCorrection(kTRUE),
77 fApplyFractionHadronicCorrection(kFALSE),
78 fFractionHadronicCorrection(0.3),
94 //____________________________________________________________________________
95 AliJetESDReader::~AliJetESDReader()
113 //____________________________________________________________________________
114 void AliJetESDReader::OpenInputFiles()
116 // chain for the ESDs
117 fChain = new TChain("esdTree");
119 // get directory and pattern name from the header
120 const char* dirName=fReaderHeader->GetDirectory();
121 const char* pattern=fReaderHeader->GetPattern();
123 // Add files matching patters to the chain
124 void *dir = gSystem->OpenDirectory(dirName);
125 const char *name = 0x0;
126 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
128 while ((name = gSystem->GetDirEntry(dir))){
129 if (a>=nesd) continue;
131 if (strstr(name,pattern)){
132 printf("Adding %s\n",name);
133 fChain->AddFile(Form("%s/%s/AliESDs.root",dirName,name));
138 gSystem->FreeDirectory(dir);
140 fESD = new AliESDEvent();
141 fESD->ReadFromTree(fChain);
143 int nMax = fChain->GetEntries();
145 printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
147 // set number of events in header
148 if (fReaderHeader->GetLastEvent() == -1)
149 fReaderHeader->SetLastEvent(nMax);
151 Int_t nUsr = fReaderHeader->GetLastEvent();
152 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
157 //__________________________________________________________
158 void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
161 // Set flag to apply MIP correction fApplyMIPCorrection
162 // - exclusive with fApplyFractionHadronicCorrection
165 fApplyMIPCorrection = val;
166 if(fApplyMIPCorrection == kTRUE)
168 SetApplyFractionHadronicCorrection(kFALSE);
169 printf("Enabling MIP Correction \n");
173 printf("Disabling MIP Correction \n");
177 //__________________________________________________________
178 void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
181 // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
182 // - exclusive with fApplyMIPCorrection
185 fApplyFractionHadronicCorrection = val;
186 if(fApplyFractionHadronicCorrection == kTRUE)
188 SetApplyMIPCorrection(kFALSE);
189 printf("Enabling Fraction Hadronic Correction \n");
193 printf("Disabling Fraction Hadronic Correction \n");
197 //__________________________________________________________
198 void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
201 // Set value to fFractionHadronicCorrection (default is 0.3)
202 // apply EMC hadronic correction fApplyFractionHadronicCorrection
203 // - exclusive with fApplyMIPCorrection
206 fFractionHadronicCorrection = val;
207 if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
209 SetApplyFractionHadronicCorrection(kTRUE);
210 printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
214 SetApplyFractionHadronicCorrection(kFALSE);
219 //____________________________________________________________________________
220 void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /*aod*/, const TObject* /*mc*/) {
222 fESD = (AliESDEvent*) esd;
225 //____________________________________________________________________________
226 Bool_t AliJetESDReader::FillMomentumArray()
228 // Fill momentum array
235 // clear momentum array
237 fDebug = fReaderHeader->GetDebug();
238 fOpt = fReaderHeader->GetDetector();
244 // get number of tracks in event (for the loop)
245 nt = fESD->GetNumberOfTracks();
246 printf("Fill Momentum Array %5d \n", nt);
248 // temporary storage of signal and pt cut flag
249 Int_t* sflag = new Int_t[nt];
250 Int_t* cflag = new Int_t[nt];
252 // get cuts set by user
253 Float_t ptMin = fReaderHeader->GetPtCut();
254 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
255 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
257 //loop over tracks in ESD
258 for (Int_t it = 0; it < nt; it++) {
259 AliESDtrack *track = fESD->GetTrack(it);
260 UInt_t status = track->GetStatus();
263 track->GetPxPyPz(mom);
264 p3.SetXYZ(mom[0],mom[1],mom[2]);
266 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
267 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
268 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
269 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
270 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
271 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
273 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
275 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
277 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
279 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
283 // set the signal flags
284 fSignalFlag.Set(goodTrack,sflag);
285 fCutFlag.Set(goodTrack,cflag);
294 //____________________________________________________________________________
295 void AliJetESDReader::CreateTasks(TChain* tree)
298 // For reader task initialization
301 fDebug = fReaderHeader->GetDebug();
302 fDZ = fReaderHeader->GetDZ();
305 // Init EMCAL geometry, if needed
308 else Info(" SetEMCALGeometry","was already done.. it's called just once !!");
312 // Create and init unit array
315 // Create global reader task for analysis
316 fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
317 // Create a task for to fill the charged particle information
318 fFillUAFromTracks = new AliJetESDFillUnitArrayTracks();
319 fFillUAFromTracks->SetReaderHeader(fReaderHeader);
320 fFillUAFromTracks->SetGeom(fGeom);
321 fFillUAFromTracks->SetTPCGrid(fTpcGrid);
322 fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
324 { // Calo dead zones inclusion
325 fFillUAFromTracks->SetGrid0(fGrid0);
326 fFillUAFromTracks->SetGrid1(fGrid1);
327 fFillUAFromTracks->SetGrid2(fGrid2);
328 fFillUAFromTracks->SetGrid3(fGrid3);
329 fFillUAFromTracks->SetGrid4(fGrid4);
331 fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
332 fFillUAFromTracks->SetHadCorrector(fHadCorr);
333 // Create a task for to fill the neutral particle information
334 fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
335 fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
336 fFillUAFromEMCalDigits->SetGeom(fGeom);
337 fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
338 fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
339 fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
340 fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
341 fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
342 // Add the task to global task
343 fFillUnitArray->Add(fFillUAFromTracks);
344 fFillUnitArray->Add(fFillUAFromEMCalDigits);
345 fFillUAFromTracks->SetActive(kFALSE);
346 fFillUAFromEMCalDigits->SetActive(kFALSE);
348 cout << "Tasks instantiated at that stage ! " << endl;
349 cout << "You can loop over events now ! " << endl;
353 //____________________________________________________________________________
354 Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
357 // Reader task execussion
361 fRefArray = refArray;
363 // clear momentum array
366 fDebug = fReaderHeader->GetDebug();
367 fOpt = fReaderHeader->GetDetector();
373 // TPC only or Digits+TPC or Clusters+TPC
374 if(fOpt%2==!0 && fOpt!=0){
375 fFillUAFromTracks->SetESD(fESD);
376 fFillUAFromTracks->SetActive(kTRUE);
377 fFillUAFromTracks->SetUnitArray(fUnitArray);
378 fFillUAFromTracks->SetRefArray(fRefArray);
379 fFillUAFromTracks->SetProcId(fProcId);
380 // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
381 fFillUAFromTracks->Exec("tpc"); // Temporarily added
383 fNumCandidate = fFillUAFromTracks->GetMult();
384 fNumCandidateCut = fFillUAFromTracks->GetMultCut();
388 // Digits only or Digits+TPC
389 if(fOpt>=2 && fOpt<=3){
390 fFillUAFromEMCalDigits->SetESD(fESD);
391 fFillUAFromEMCalDigits->SetActive(kTRUE);
392 fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
393 fFillUAFromEMCalDigits->SetRefArray(fRefArray);
394 fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
395 fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
396 fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
397 fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added
398 fNumCandidate = fFillUAFromEMCalDigits->GetMult();
399 fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
402 // fFillUnitArray->ExecuteTask(); // Temporarily commented
408 //____________________________________________________________________________
409 void AliJetESDReader::InitParameters()
411 // Initialise parameters
412 fOpt = fReaderHeader->GetDetector();
413 // fHadCorr = 0; // For hadron correction
415 if(fOpt==0 || fOpt==1)
416 fECorrection = 0; // For electron correction
417 else fECorrection = 1; // For electron correction
419 fNumUnits = fGeom->GetEMCGeometry()->GetNCells(); // Number of cells in EMCAL
420 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
423 //____________________________________________________________________________
424 void AliJetESDReader::InitUnitArray()
427 // Create and Initialise unit arrays
430 Int_t nElements = fTpcGrid->GetNEntries();
431 Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
432 if(fArrayInitialised) fUnitArray->Delete();
434 if(fTpcGrid->GetGridType()==0)
435 { // Fill the following quantities :
436 // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
437 // detector flag, in/out jet, pt cut, mass, cluster ID)
438 for(Int_t nBin = 1; nBin < nElements+1; nBin++)
440 fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
441 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
442 Deta = fTpcGrid->GetDeta();
443 Dphi = fTpcGrid->GetDphi();
444 new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
448 if(fTpcGrid->GetGridType()==1)
451 Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
455 // Define a grid of cell for the gaps between SM
456 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
457 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
458 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
459 fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
460 fGrid0->SetGridType(0);
461 fGrid0->SetMatrixIndexes();
462 fGrid0->SetIndexIJ();
463 n0 = fGrid0->GetNEntries();
464 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
465 fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
466 fGrid1->SetGridType(0);
467 fGrid1->SetMatrixIndexes();
468 fGrid1->SetIndexIJ();
469 n1 = fGrid1->GetNEntries();
470 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
471 fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
472 fGrid2->SetGridType(0);
473 fGrid2->SetMatrixIndexes();
474 fGrid2->SetIndexIJ();
475 n2 = fGrid2->GetNEntries();
476 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
477 fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
478 fGrid3->SetGridType(0);
479 fGrid3->SetMatrixIndexes();
480 fGrid3->SetIndexIJ();
481 n3 = fGrid3->GetNEntries();
482 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
483 fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
484 fGrid4->SetGridType(0);
485 fGrid4->SetMatrixIndexes();
486 fGrid4->SetIndexIJ();
487 n4 = fGrid4->GetNEntries();
489 nGaps = n0+n1+n2+n3+n4;
493 for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
497 fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
498 // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
499 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
500 Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
501 Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
502 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
505 if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
506 fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
507 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
508 Deta = fTpcGrid->GetDeta();
509 Dphi = fTpcGrid->GetDphi();
510 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
514 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
515 if(nBin<fNumUnits+nElements+n0)
517 Float_t phi0 = eta = 0.;
518 fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
519 Deta = fGrid0->GetDeta();
520 Dphi = fGrid0->GetDphi();
521 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
523 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
525 Float_t phi0 = eta = 0.;
526 fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta);
527 Deta = fGrid1->GetDeta();
528 Dphi = fGrid1->GetDphi();
529 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
531 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
533 Float_t phi0 = eta = 0.;
534 fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
535 Deta = fGrid2->GetDeta();
536 Dphi = fGrid2->GetDphi();
537 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
539 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
541 Float_t phi0 = eta = 0.;
542 fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
543 Deta = fGrid3->GetDeta();
544 Dphi = fGrid3->GetDphi();
545 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
547 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
549 Float_t phi0 = eta = 0.;
550 fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
551 Deta = fGrid4->GetDeta();
552 Dphi = fGrid4->GetDphi();
553 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
559 } // end loop on nBin
560 } // end grid type == 1
561 fArrayInitialised = 1;