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