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 "AliJetFillUnitArrayTracks.h"
49 #include "AliJetFillUnitArrayEMCalDigits.h"
50 #include "AliJetUnitArray.h"
51 #include "AliAnalysisTask.h"
53 ClassImp(AliJetESDReader)
55 AliJetESDReader::AliJetESDReader():
89 //____________________________________________________________________________
90 AliJetESDReader::~AliJetESDReader()
108 //____________________________________________________________________________
109 void AliJetESDReader::OpenInputFiles()
111 // chain for the ESDs
112 fChain = new TChain("esdTree");
114 // get directory and pattern name from the header
115 const char* dirName=fReaderHeader->GetDirectory();
116 const char* pattern=fReaderHeader->GetPattern();
118 // Add files matching patters to the chain
119 void *dir = gSystem->OpenDirectory(dirName);
120 const char *name = 0x0;
121 int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
123 while ((name = gSystem->GetDirEntry(dir))){
124 if (a>=nesd) continue;
126 if (strstr(name,pattern)){
127 printf("Adding %s\n",name);
129 sprintf(path,"%s/%s/AliESDs.root",dirName,name);
130 fChain->AddFile(path);
135 gSystem->FreeDirectory(dir);
137 fESD = new AliESDEvent();
138 fESD->ReadFromTree(fChain);
140 int nMax = fChain->GetEntries();
142 printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
144 // set number of events in header
145 if (fReaderHeader->GetLastEvent() == -1)
146 fReaderHeader->SetLastEvent(nMax);
148 Int_t nUsr = fReaderHeader->GetLastEvent();
149 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
154 //____________________________________________________________________________
155 void AliJetESDReader::SetInputEvent(TObject* esd, TObject* /*aod*/, TObject* /*mc*/) {
157 fESD = (AliESDEvent*) esd;
160 //____________________________________________________________________________
161 Bool_t AliJetESDReader::FillMomentumArray()
163 // Fill momentum array
170 // clear momentum array
172 fDebug = fReaderHeader->GetDebug();
173 fOpt = fReaderHeader->GetDetector();
179 // get number of tracks in event (for the loop)
180 nt = fESD->GetNumberOfTracks();
181 printf("Fill Momentum Array %5d \n", nt);
183 // temporary storage of signal and pt cut flag
184 Int_t* sflag = new Int_t[nt];
185 Int_t* cflag = new Int_t[nt];
187 // get cuts set by user
188 Float_t ptMin = fReaderHeader->GetPtCut();
189 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
190 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
192 //loop over tracks in ESD
193 for (Int_t it = 0; it < nt; it++) {
194 AliESDtrack *track = fESD->GetTrack(it);
195 UInt_t status = track->GetStatus();
198 track->GetPxPyPz(mom);
199 p3.SetXYZ(mom[0],mom[1],mom[2]);
201 if ((status & AliESDtrack::kTPCrefit) == 0) continue; // quality check
202 if ((status & AliESDtrack::kITSrefit) == 0) continue; // quality check
203 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
204 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
205 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
206 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
208 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
210 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
212 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
214 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
218 // set the signal flags
219 fSignalFlag.Set(goodTrack,sflag);
220 fCutFlag.Set(goodTrack,cflag);
229 //____________________________________________________________________________
230 void AliJetESDReader::CreateTasks(TChain* tree)
233 // For reader task initialization
236 fDebug = fReaderHeader->GetDebug();
237 fDZ = fReaderHeader->GetDZ();
240 // Init EMCAL geometry
244 // Create and init unit array
247 // Create global reader task for analysis
248 fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
249 // Create a task for to fill the charged particle information
250 fFillUAFromTracks = new AliJetFillUnitArrayTracks();
251 fFillUAFromTracks->SetReaderHeader(fReaderHeader);
252 fFillUAFromTracks->SetGeom(fGeom);
253 fFillUAFromTracks->SetTPCGrid(fTpcGrid);
254 fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
256 { // Calo dead zones inclusion
257 fFillUAFromTracks->SetGrid0(fGrid0);
258 fFillUAFromTracks->SetGrid1(fGrid1);
259 fFillUAFromTracks->SetGrid2(fGrid2);
260 fFillUAFromTracks->SetGrid3(fGrid3);
261 fFillUAFromTracks->SetGrid4(fGrid4);
263 fFillUAFromTracks->SetHadCorrection(fHCorrection);
264 fFillUAFromTracks->SetHadCorrector(fHadCorr);
265 // Create a task for to fill the neutral particle information
266 fFillUAFromEMCalDigits = new AliJetFillUnitArrayEMCalDigits();
267 fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
268 fFillUAFromEMCalDigits->SetGeom(fGeom);
269 fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
270 fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
271 fFillUAFromEMCalDigits->SetEleCorrection(fECorrection);
272 // Add the task to global task
273 fFillUnitArray->Add(fFillUAFromTracks);
274 fFillUnitArray->Add(fFillUAFromEMCalDigits);
275 fFillUAFromTracks->SetActive(kFALSE);
276 fFillUAFromEMCalDigits->SetActive(kFALSE);
278 cout << "Tasks instantiated at that stage ! " << endl;
279 cout << "You can loop over events now ! " << endl;
283 //____________________________________________________________________________
284 Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
287 // Reader task execussion
291 fRefArray = refArray;
292 vector<Float_t> vtmp(3);
294 // clear momentum array
297 fDebug = fReaderHeader->GetDebug();
298 fOpt = fReaderHeader->GetDetector();
304 // TPC only or Digits+TPC or Clusters+TPC
305 if(fOpt%2==!0 && fOpt!=0){
306 fFillUAFromTracks->SetESD(fESD);
307 fFillUAFromTracks->SetActive(kTRUE);
308 fFillUAFromTracks->SetUnitArray(fUnitArray);
309 fFillUAFromTracks->SetRefArray(fRefArray);
310 fFillUAFromTracks->SetProcId(fProcId);
311 // fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
312 fFillUAFromTracks->Exec("tpc"); // Temporarily added
314 fNumCandidate = fFillUAFromTracks->GetMult();
315 fNumCandidateCut = fFillUAFromTracks->GetMultCut();
319 // Digits only or Digits+TPC
320 if(fOpt>=2 && fOpt<=3){
321 fFillUAFromEMCalDigits->SetESD(fESD);
322 fFillUAFromEMCalDigits->SetActive(kTRUE);
323 fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
324 fFillUAFromEMCalDigits->SetRefArray(fRefArray);
325 fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
326 fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
327 fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
328 fFillUAFromEMCalDigits->Exec("digits"); // Temporarily added
329 fNumCandidate = fFillUAFromEMCalDigits->GetMult();
330 fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
333 // fFillUnitArray->ExecuteTask(); // Temporarily commented
338 //____________________________________________________________________________
339 Bool_t AliJetESDReader::SetEMCALGeometry()
342 // Set the EMCal Geometry
345 if (!fTree->GetFile())
348 TString geomFile(fTree->GetFile()->GetName());
349 geomFile.ReplaceAll("AliESDs", "geometry");
351 // temporary workaround for PROOF bug #18505
352 geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
353 if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
355 // Define EMCAL geometry to be able to read ESDs
356 fGeom = AliJetDummyGeo::GetInstance();
358 fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
360 // To be setted to run some AliEMCALGeometry functions
361 TGeoManager::Import(geomFile);
362 fGeom->GetTransformationForSM();
363 printf("\n EMCal Geometry set ! \n");
368 //____________________________________________________________________________
369 void AliJetESDReader::InitParameters()
371 // Initialise parameters
372 fOpt = fReaderHeader->GetDetector();
373 fHadCorr = 0; // For hadron correction
375 if(fOpt==0 || fOpt==1)
376 fECorrection = 0; // For electron correction
377 else fECorrection = 1; // For electron correction
379 fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
380 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
383 //____________________________________________________________________________
384 void AliJetESDReader::InitUnitArray()
387 // Create and Initialise unit arrays
390 Int_t nElements = fTpcGrid->GetNEntries();
391 Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
392 if(fArrayInitialised) fUnitArray->Delete();
394 if(fTpcGrid->GetGridType()==0)
395 { // Fill the following quantities :
396 // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
397 // detector flag, in/out jet, pt cut, mass, cluster ID)
398 for(Int_t nBin = 1; nBin < nElements+1; nBin++)
400 fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
401 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
402 Deta = fTpcGrid->GetDeta();
403 Dphi = fTpcGrid->GetDphi();
404 new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
408 if(fTpcGrid->GetGridType()==1)
411 Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
415 // Define a grid of cell for the gaps between SM
416 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
417 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
418 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
419 fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
420 fGrid0->SetGridType(0);
421 fGrid0->SetMatrixIndexes();
422 fGrid0->SetIndexIJ();
423 n0 = fGrid0->GetNEntries();
424 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
425 fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
426 fGrid1->SetGridType(0);
427 fGrid1->SetMatrixIndexes();
428 fGrid1->SetIndexIJ();
429 n1 = fGrid1->GetNEntries();
430 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
431 fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
432 fGrid2->SetGridType(0);
433 fGrid2->SetMatrixIndexes();
434 fGrid2->SetIndexIJ();
435 n2 = fGrid2->GetNEntries();
436 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
437 fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
438 fGrid3->SetGridType(0);
439 fGrid3->SetMatrixIndexes();
440 fGrid3->SetIndexIJ();
441 n3 = fGrid3->GetNEntries();
442 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
443 fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
444 fGrid4->SetGridType(0);
445 fGrid4->SetMatrixIndexes();
446 fGrid4->SetIndexIJ();
447 n4 = fGrid4->GetNEntries();
451 cout << "n0 cells: " << n0 << "phimin0: " << phimin0 << ", phimax0: " << phimax0 << endl;
452 cout << "n1 cells: " << n1 << "phimin1: " << phimin1 << ", phimax1: " << phimax1 << endl;
453 cout << "n2 cells: " << n2 << "phimin2: " << phimin2 << ", phimax2: " << phimax2 << endl;
454 cout << "n3 cells: " << n3 << "phimin3: " << phimin3 << ", phimax3: " << phimax3 << endl;
455 cout << "n4 cells: " << n4 << "phimin4: " << phimin4 << ", phimax4: " << phimax4 << endl;
458 nGaps = n0+n1+n2+n3+n4;
462 for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
466 fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
467 // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
468 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
469 Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
470 Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
471 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
474 if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
475 fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
476 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
477 Deta = fTpcGrid->GetDeta();
478 Dphi = fTpcGrid->GetDphi();
479 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
483 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
484 if(nBin<fNumUnits+nElements+n0)
486 Float_t phi0 = eta = 0.;
487 fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
488 Deta = fGrid0->GetDeta();
489 Dphi = fGrid0->GetDphi();
490 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
492 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
494 Float_t phi = eta = 0.;
495 fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
496 Deta = fGrid1->GetDeta();
497 Dphi = fGrid1->GetDphi();
498 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
500 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
502 Float_t phi0 = eta = 0.;
503 fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
504 Deta = fGrid2->GetDeta();
505 Dphi = fGrid2->GetDphi();
506 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
508 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
510 Float_t phi0 = eta = 0.;
511 fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
512 Deta = fGrid3->GetDeta();
513 Dphi = fGrid3->GetDphi();
514 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
516 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
518 Float_t phi0 = eta = 0.;
519 fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
520 Deta = fGrid4->GetDeta();
521 Dphi = fGrid4->GetDphi();
522 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
528 } // end loop on nBin
529 } // end grid type == 1
530 fArrayInitialised = 1;