removing unused variable fPtCut to avoid confusion, same is uses in the ReaderHeader
[u/mrichter/AliRoot.git] / JETAN / AliJetAODReader.cxx
CommitLineData
586f2bc3 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
16//-------------------------------------------------------------------------
17// Jet AOD Reader
18// AOD reader for jet analysis
e53baffe 19// This is the reader which must be used if the jet analysis task
20// is executed after the ESD filter task, in order to read its output
21//
586f2bc3 22// Author: Davide Perrino <davide.perrino@cern.ch>
23//-------------------------------------------------------------------------
24
25
26#include <Riostream.h>
27#include <TSystem.h>
28#include <TLorentzVector.h>
29#include <TVector3.h>
ee7de0dd 30#include <TChain.h>
be6e5811 31#include <TFile.h>
32#include <TTask.h>
33#include <TGeoManager.h>
ee7de0dd 34
586f2bc3 35#include "AliJetAODReader.h"
36#include "AliJetAODReaderHeader.h"
37#include "AliAODEvent.h"
38#include "AliAODTrack.h"
0e4398c8 39#include "AliAODMCParticle.h"
be6e5811 40#include "AliJetDummyGeo.h"
41#include "AliJetAODFillUnitArrayTracks.h"
42#include "AliJetAODFillUnitArrayEMCalDigits.h"
e36a3f22 43#include "AliJetHadronCorrectionv1.h"
be6e5811 44#include "AliJetUnitArray.h"
586f2bc3 45
46ClassImp(AliJetAODReader)
47
48AliJetAODReader::AliJetAODReader():
49 AliJetReader(),
586f2bc3 50 fAOD(0x0),
51 fRef(new TRefArray),
52 fDebug(0),
be6e5811 53 fOpt(0),
54 fGeom(0),
55 fHadCorr(0x0),
56 fTpcGrid(0x0),
57 fEmcalGrid(0x0),
58 fGrid0(0),
59 fGrid1(0),
60 fGrid2(0),
61 fGrid3(0),
62 fGrid4(0),
be6e5811 63 fApplyElectronCorrection(kFALSE),
be6e5811 64 fApplyMIPCorrection(kTRUE),
65 fApplyFractionHadronicCorrection(kFALSE),
66 fFractionHadronicCorrection(0.3),
67 fNumUnits(0),
68 fMass(0),
69 fSign(0),
70 fNIn(0),
71 fDZ(0),
72 fNeta(0),
73 fNphi(0),
be6e5811 74 fRefArray(0x0),
75 fProcId(kFALSE)
586f2bc3 76{
77 // Constructor
78}
79
80//____________________________________________________________________________
81
82AliJetAODReader::~AliJetAODReader()
83{
84 // Destructor
586f2bc3 85 delete fAOD;
86 delete fRef;
be6e5811 87 delete fTpcGrid;
88 delete fEmcalGrid;
89 if(fDZ)
90 {
91 delete fGrid0;
92 delete fGrid1;
93 delete fGrid2;
94 delete fGrid3;
95 delete fGrid4;
96 }
97
586f2bc3 98}
99
100//____________________________________________________________________________
101
102void AliJetAODReader::OpenInputFiles()
103{
104 // Open the necessary input files
105 // chain for the AODs
106 fChain = new TChain("aodTree");
107
108 // get directory and pattern name from the header
109 const char* dirName=fReaderHeader->GetDirectory();
110 const char* pattern=fReaderHeader->GetPattern();
111
112// // Add files matching patters to the chain
113
114 void *dir = gSystem->OpenDirectory(dirName);
115 const char *name = 0x0;
116 int naod = ((AliJetAODReaderHeader*) fReaderHeader)->GetNaod();
117 int a = 0;
118 while ((name = gSystem->GetDirEntry(dir))){
119 if (a>=naod) continue;
120
121 if (strstr(name,pattern)){
122 char path[256];
123 sprintf(path,"%s/%s/aod.root",dirName,name);
124 fChain->AddFile(path);
125 a++;
126 }
127 }
128
129 gSystem->FreeDirectory(dir);
130
586f2bc3 131 fAOD = 0;
132 fChain->SetBranchAddress("AOD",&fAOD);
133
134 int nMax = fChain->GetEntries();
135
be6e5811 136 printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
586f2bc3 137
138 // set number of events in header
139 if (fReaderHeader->GetLastEvent() == -1)
140 fReaderHeader->SetLastEvent(nMax);
141 else {
142 Int_t nUsr = fReaderHeader->GetLastEvent();
143 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
144 }
145}
146
147//____________________________________________________________________________
148
149void AliJetAODReader::ConnectTree(TTree* tree, TObject* /*data*/) {
150 // Connect the tree
151 // For AOD reader it's needed only to set the number of events
152 fChain = (TChain*) tree;
153
154 Int_t nMax = fChain->GetEntries();
be6e5811 155 printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
586f2bc3 156 // set number of events in header
157 if (fReaderHeader->GetLastEvent() == -1)
158 fReaderHeader->SetLastEvent(nMax);
159 else {
160 Int_t nUsr = fReaderHeader->GetLastEvent();
161 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
162 }
163}
164
0e4398c8 165
166Bool_t AliJetAODReader::AcceptAODMCParticle(AliAODMCParticle *mcP,Short_t flag){
167
168 //
169 // filter for charge and for charged and neutral, no detector
170 // response yet
171 // physical priamries already selected
172
173 Int_t pdg = TMath::Abs(mcP->GetPdgCode());
174
175 // exclude neutrinos anyway
176 if((pdg == 12 || pdg == 14 || pdg == 16)) return kFALSE;
177
178 if(flag==AliJetAODReaderHeader::kAllMC)return kTRUE;
179 if(flag==AliJetAODReaderHeader::kChargedMC){
180 if(mcP->Charge()==0)return kFALSE;
181 return kTRUE;
182 }
183
184 return kTRUE;
185
186}
187
188
189Bool_t AliJetAODReader::FillMomentumArrayMC(){
190
191 //
192 // This routine fetches the MC particles from the AOD
193 // Depending on the flag all particles except neurinos are use
194 // or only charged particles
195 //
196
197 TClonesArray *mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
198 if(!mcArray){
199 Printf("%s:%d No MC particle branch found",(char*)__FILE__,__LINE__);
200 return kFALSE;
201 }
202
203 Int_t nMC = mcArray->GetEntriesFast();
204
205 // get number of tracks in event (for the loop)
206 if(fDebug>0)printf("AOD MC tracks: %5d \t", nMC);
207
208 // temporary storage of signal and pt cut flag
209 Int_t* sflag = new Int_t[nMC];
210 Int_t* cflag = new Int_t[nMC];
211
212 // get cuts set by user
213 Float_t ptMin = fReaderHeader->GetPtCut();
214 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
215 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
216 Int_t mcTrack = 0;
217 Float_t pt, eta;
218 TVector3 p3;
219
220 Short_t readerFlag = ((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC();
221
222
223 for (Int_t it = 0; it < nMC; it++) {
224 AliAODMCParticle *track = (AliAODMCParticle*)mcArray->At(it);
225 if(!track->IsPhysicalPrimary())continue;
226
227 p3.SetXYZ(track->Px(),track->Py(),track->Pz());
228 eta = p3.Eta();
229 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
230 if(!AcceptAODMCParticle(track,readerFlag))continue;
231 pt = p3.Pt();
232
233
234
235
f7ce54e3 236 if (mcTrack == 0){
237 fRef->Delete(); // make sure to delete before placement new...
238 new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
239 }
0e4398c8 240 new ((*fMomentumArray)[mcTrack]) TLorentzVector(p3,p3.Mag());
241 sflag[mcTrack] = 1;
242 cflag[mcTrack] = ( pt > ptMin ) ? 1: 0;
243 fRef->Add(track);
244 mcTrack++;
245 }
246 if(fDebug>0)printf("Used MC tracks: %5d \n", mcTrack);
247 // set the signal flags
248 fSignalFlag.Set(mcTrack,sflag);
249 fCutFlag.Set(mcTrack,cflag);
250
251 delete [] sflag;
252 delete [] cflag;
253
254 return kTRUE;
255}
256
586f2bc3 257//____________________________________________________________________________
258
9e4cc50d 259Bool_t AliJetAODReader::FillMomentumArray()
586f2bc3 260{
261 // Clear momentum array
262 ClearArray();
e53baffe 263 fRef->Clear();
586f2bc3 264 fDebug = fReaderHeader->GetDebug();
0e4398c8 265
586f2bc3 266 if (!fAOD) {
267 return kFALSE;
268 }
f7ce54e3 269
270 if(((AliJetAODReaderHeader*)fReaderHeader)->GetReadAODMC()){
271 return FillMomentumArrayMC();
272 }
273
586f2bc3 274 // get number of tracks in event (for the loop)
275 Int_t nt = fAOD->GetNTracks();
0e4398c8 276 if(fDebug>0)printf("AOD tracks: %5d \t", nt);
586f2bc3 277
278 // temporary storage of signal and pt cut flag
279 Int_t* sflag = new Int_t[nt];
280 Int_t* cflag = new Int_t[nt];
281
282 // get cuts set by user
283 Float_t ptMin = fReaderHeader->GetPtCut();
284 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
285 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
d61f057d 286 UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
e53baffe 287
586f2bc3 288 //loop over tracks
289 Int_t aodTrack = 0;
290 Float_t pt, eta;
291 TVector3 p3;
d61f057d 292
586f2bc3 293 for (Int_t it = 0; it < nt; it++) {
d61f057d 294 AliAODTrack *track = fAOD->GetTrack(it);
e53baffe 295 UInt_t status = track->GetStatus();
d61f057d 296
297 Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
298 p3.SetXYZ(mom[0],mom[1],mom[2]);
299 pt = p3.Pt();
300 eta = p3.Eta();
e53baffe 301 if (status == 0) continue;
d61f057d 302 if((filterMask>0)&&!(track->TestFilterBit(filterMask)))continue;
d61f057d 303 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
e53baffe 304
6473bac6 305 if (aodTrack == 0){
306 fRef->Delete(); // make sure to delete before placement new...
307 new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
308 }
e53baffe 309 new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
310 sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
311 cflag[aodTrack] = ( pt > ptMin ) ? 1: 0;
312 aodTrack++;
d61f057d 313 fRef->Add(track);
586f2bc3 314 }
0e4398c8 315 if(fDebug>0)printf("Used AOD tracks: %5d \n", aodTrack);
586f2bc3 316 // set the signal flags
317 fSignalFlag.Set(aodTrack,sflag);
318 fCutFlag.Set(aodTrack,cflag);
d61f057d 319
680ed75f 320 delete [] sflag;
321 delete [] cflag;
d61f057d 322
586f2bc3 323 return kTRUE;
324}
be6e5811 325
326//__________________________________________________________
327void AliJetAODReader::SetApplyMIPCorrection(Bool_t val)
328{
329 //
330 // Set flag to apply MIP correction fApplyMIPCorrection
331 // - exclusive with fApplyFractionHadronicCorrection
332 //
333
334 fApplyMIPCorrection = val;
335 if(fApplyMIPCorrection == kTRUE)
336 {
337 SetApplyFractionHadronicCorrection(kFALSE);
338 printf("Enabling MIP Correction \n");
339 }
340 else
341 {
342 printf("Disabling MIP Correction \n");
343 }
344}
345
346//__________________________________________________________
347void AliJetAODReader::SetApplyFractionHadronicCorrection(Bool_t val)
348{
349 //
350 // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
351 // - exclusive with fApplyMIPCorrection
352 //
353
354 fApplyFractionHadronicCorrection = val;
355 if(fApplyFractionHadronicCorrection == kTRUE)
356 {
357 SetApplyMIPCorrection(kFALSE);
358 printf("Enabling Fraction Hadronic Correction \n");
359 }
360 else
361 {
362 printf("Disabling Fraction Hadronic Correction \n");
363 }
364}
365
366//__________________________________________________________
367void AliJetAODReader::SetFractionHadronicCorrection(Double_t val)
368{
369 //
370 // Set value to fFractionHadronicCorrection (default is 0.3)
371 // apply EMC hadronic correction fApplyFractionHadronicCorrection
372 // - exclusive with fApplyMIPCorrection
373 //
374
375 fFractionHadronicCorrection = val;
376 if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
377 {
378 SetApplyFractionHadronicCorrection(kTRUE);
379 printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
380 }
381 else
382 {
383 SetApplyFractionHadronicCorrection(kFALSE);
384 }
385}
386
387//____________________________________________________________________________
388void AliJetAODReader::CreateTasks(TChain* tree)
389{
4751efb5 390 //
391 // For reader task initialization
392 //
393
be6e5811 394 fDebug = fReaderHeader->GetDebug();
395 fDZ = fReaderHeader->GetDZ();
396 fTree = tree;
397
398 // Init EMCAL geometry and create UnitArray object
399 SetEMCALGeometry();
400 // cout << "In create task" << endl;
401 InitParameters();
402 InitUnitArray();
403
404 fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
405 fFillUAFromTracks = new AliJetAODFillUnitArrayTracks();
406 fFillUAFromTracks->SetReaderHeader(fReaderHeader);
407 fFillUAFromTracks->SetGeom(fGeom);
408 fFillUAFromTracks->SetTPCGrid(fTpcGrid);
409 fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
410
411 if(fDZ)
412 {
413 fFillUAFromTracks->SetGrid0(fGrid0);
414 fFillUAFromTracks->SetGrid1(fGrid1);
415 fFillUAFromTracks->SetGrid2(fGrid2);
416 fFillUAFromTracks->SetGrid3(fGrid3);
417 fFillUAFromTracks->SetGrid4(fGrid4);
418 }
419 fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
420 fFillUAFromTracks->SetHadCorrector(fHadCorr);
421 fFillUAFromEMCalDigits = new AliJetAODFillUnitArrayEMCalDigits();
422 fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
423 fFillUAFromEMCalDigits->SetGeom(fGeom);
424 fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
425 fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
426 fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
427 fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
428 fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
e36a3f22 429 // Add the task to global task
be6e5811 430 fFillUnitArray->Add(fFillUAFromTracks);
431 fFillUnitArray->Add(fFillUAFromEMCalDigits);
432 fFillUAFromTracks->SetActive(kFALSE);
433 fFillUAFromEMCalDigits->SetActive(kFALSE);
434
435 cout << "Tasks instantiated at that stage ! " << endl;
436 cout << "You can loop over events now ! " << endl;
437
438}
439
440//____________________________________________________________________________
441Bool_t AliJetAODReader::ExecTasks(Bool_t procid, TRefArray* refArray)
442{
443 //
444 // Main function
445 // Fill the reader part
446 //
447
448 fProcId = procid;
449 fRefArray = refArray;
450//(not used ?) Int_t nEntRef = fRefArray->GetEntries();
451//(not used ?) Int_t nEntUnit = fUnitArray->GetEntries();
452
453 // clear momentum array
454 ClearArray();
e36a3f22 455 fRef->Clear();
be6e5811 456
457 fDebug = fReaderHeader->GetDebug();
458 fOpt = fReaderHeader->GetDetector();
459
460 if(!fAOD) {
461 return kFALSE;
462 }
463
464 // TPC only or Digits+TPC or Clusters+TPC
465 if(fOpt%2==!0 && fOpt!=0){
466 fFillUAFromTracks->SetAOD(fAOD);
467 fFillUAFromTracks->SetActive(kTRUE);
468 fFillUAFromTracks->SetUnitArray(fUnitArray);
469 fFillUAFromTracks->SetRefArray(fRefArray);
e36a3f22 470 fFillUAFromTracks->SetReferences(fRef);
471 fFillUAFromTracks->SetSignalFlag(fSignalFlag);
472 fFillUAFromTracks->SetCutFlag(fCutFlag);
be6e5811 473 fFillUAFromTracks->SetProcId(fProcId);
474 // fFillUAFromTracks->ExecuteTask("tpc"); // => Temporarily changed
475 fFillUAFromTracks->Exec("tpc");
476 if(fOpt==1){
477 fNumCandidate = fFillUAFromTracks->GetMult();
478 fNumCandidateCut = fFillUAFromTracks->GetMultCut();
479 }
e36a3f22 480 fSignalFlag = fFillUAFromTracks->GetSignalFlag();
481 fCutFlag = fFillUAFromTracks->GetCutFlag();
be6e5811 482 }
483
484 // Digits only or Digits+TPC
485 if(fOpt>=2 && fOpt<=3){
486 fFillUAFromEMCalDigits->SetAOD(fAOD);
487 fFillUAFromEMCalDigits->SetActive(kTRUE);
488 fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
489 fFillUAFromEMCalDigits->SetRefArray(fRefArray);
490 fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
491 fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
492 fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
493 fFillUAFromEMCalDigits->Exec("digits"); // => Temporarily added
494 fNumCandidate = fFillUAFromEMCalDigits->GetMult();
495 fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
496 }
497
498 // fFillUnitArray->ExecuteTask(); // => Temporarily commented
499
500 return kTRUE;
501}
502
503//____________________________________________________________________________
504Bool_t AliJetAODReader::SetEMCALGeometry()
505{
506 //
507 // Set the EMCal Geometry
508 //
509
510 if (!fTree->GetFile())
511 return kFALSE;
512
513 TString geomFile(fTree->GetFile()->GetName());
514 geomFile.ReplaceAll("AliESDs", "geometry");
515
516 // temporary workaround for PROOF bug #18505
517 geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
518 if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
519
520 // Define EMCAL geometry to be able to read ESDs
521 fGeom = AliJetDummyGeo::GetInstance();
522 if (fGeom == 0)
523 fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
524
525 // To be setted to run some AliEMCALGeometry functions
526 TGeoManager::Import(geomFile);
527 fGeom->GetTransformationForSM();
528 printf("\n EMCal Geometry set ! \n");
529
530 return kTRUE;
531
532}
533
534//____________________________________________________________________________
535void AliJetAODReader::InitParameters()
536{
537 // Initialise parameters
538 fOpt = fReaderHeader->GetDetector();
be6e5811 539 fNumUnits = fGeom->GetNCells(); // Number of cells in EMCAL
540 if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
541}
542
543//____________________________________________________________________________
544void AliJetAODReader::InitUnitArray()
545{
546 //Initialises unit arrays
547 Int_t nElements = fTpcGrid->GetNEntries();
f23f7908 548 Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
be6e5811 549 if(fArrayInitialised) fUnitArray->Delete();
550
551 if(fTpcGrid->GetGridType()==0)
552 { // Fill the following quantities :
f23f7908 553 // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
be6e5811 554 // detector flag, in/out jet, pt cut, mass, cluster ID)
555 for(Int_t nBin = 1; nBin < nElements+1; nBin++)
556 {
557 // fTpcGrid->GetEtaPhiFromIndex2(nBin,eta,phi);
558 fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
559 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908 560 deltaEta = fTpcGrid->GetDeta();
561 deltaPhi = fTpcGrid->GetDphi();
562 new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 563 }
564 }
565
566 if(fTpcGrid->GetGridType()==1)
567 {
568 Int_t nGaps = 0;
569 Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
570
571 if(fDZ)
572 {
573 // Define a grid of cell for the gaps between SM
574 Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
575 Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
576 fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
577 fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
578 fGrid0->SetGridType(0);
579 fGrid0->SetMatrixIndexes();
580 fGrid0->SetIndexIJ();
581 n0 = fGrid0->GetNEntries();
582 fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
583 fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
584 fGrid1->SetGridType(0);
585 fGrid1->SetMatrixIndexes();
586 fGrid1->SetIndexIJ();
587 n1 = fGrid1->GetNEntries();
588 fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
589 fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
590 fGrid2->SetGridType(0);
591 fGrid2->SetMatrixIndexes();
592 fGrid2->SetIndexIJ();
593 n2 = fGrid2->GetNEntries();
594 fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
595 fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
596 fGrid3->SetGridType(0);
597 fGrid3->SetMatrixIndexes();
598 fGrid3->SetIndexIJ();
599 n3 = fGrid3->GetNEntries();
600 fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
601 fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
602 fGrid4->SetGridType(0);
603 fGrid4->SetMatrixIndexes();
604 fGrid4->SetIndexIJ();
605 n4 = fGrid4->GetNEntries();
606
607 nGaps = n0+n1+n2+n3+n4;
608
609 }
610
611 for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
612 {
613 if(nBin<fNumUnits)
614 {
615 fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
616 // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
617 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908 618 deltaEta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
619 deltaPhi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
620 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 621 }
622 else {
623 if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
624 fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
625 phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
f23f7908 626 deltaEta = fTpcGrid->GetDeta();
627 deltaPhi = fTpcGrid->GetDphi();
628 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 629 }
630 else {
631 if(fDZ) {
632 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
633 if(nBin<fNumUnits+nElements+n0)
634 {
f23f7908 635 phi = eta = 0.;
be6e5811 636 fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
f23f7908 637 deltaEta = fGrid0->GetDeta();
638 deltaPhi = fGrid0->GetDphi();
639 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 640 }
641 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
642 {
f23f7908 643 phi = eta = 0.;
be6e5811 644 fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
f23f7908 645 deltaEta = fGrid1->GetDeta();
646 deltaPhi = fGrid1->GetDphi();
647 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 648 }
649 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
650 {
f23f7908 651 phi = eta = 0.;
be6e5811 652 fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
f23f7908 653 deltaEta = fGrid2->GetDeta();
654 deltaPhi = fGrid2->GetDphi();
655 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 656 }
657 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
658 {
f23f7908 659 phi = eta = 0.;
be6e5811 660 fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
f23f7908 661 deltaEta = fGrid3->GetDeta();
662 deltaPhi = fGrid3->GetDphi();
663 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 664 }
665 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
666 {
f23f7908 667 phi = eta = 0.;
be6e5811 668 fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
f23f7908 669 deltaEta = fGrid4->GetDeta();
670 deltaPhi = fGrid4->GetDphi();
671 new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,deltaEta,deltaPhi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
be6e5811 672 }
673 }
674 } // end if(fDZ)
675 } // end else 2
676 } // end else 1
677 } // end loop on nBin
678 } // end grid type == 1
679 fArrayInitialised = 1;
680}
0e4398c8 681