Protections for coverity: DIVIDE_BY_ZERO
[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>
f45e6f72 23//
24// **February 2011
f445ed00 25// implemented standard geometry (AliEMCALGeometry) instead of dummy one (AliJetDummyGeo)
f45e6f72 26// moved geometry definition in AliJetReader
27// marco.bregant@subatech.in2p3.fr
586f2bc3 28//-------------------------------------------------------------------------
29
30
31#include <Riostream.h>
32#include <TSystem.h>
33#include <TLorentzVector.h>
34#include <TVector3.h>
ee7de0dd 35#include <TChain.h>
be6e5811 36#include <TFile.h>
37#include <TTask.h>
38#include <TGeoManager.h>
f45e6f72 39#include <TGeoMatrix.h>
ee7de0dd 40
586f2bc3 41#include "AliJetAODReader.h"
42#include "AliJetAODReaderHeader.h"
43#include "AliAODEvent.h"
44#include "AliAODTrack.h"
0e4398c8 45#include "AliAODMCParticle.h"
f445ed00 46#include "AliEMCALGeometry.h"
be6e5811 47#include "AliJetAODFillUnitArrayTracks.h"
48#include "AliJetAODFillUnitArrayEMCalDigits.h"
e36a3f22 49#include "AliJetHadronCorrectionv1.h"
be6e5811 50#include "AliJetUnitArray.h"
f45e6f72 51#include "AliOADBContainer.h"
586f2bc3 52
53ClassImp(AliJetAODReader)
54
55AliJetAODReader::AliJetAODReader():
56 AliJetReader(),
586f2bc3 57 fAOD(0x0),
58 fRef(new TRefArray),
59 fDebug(0),
be6e5811 60 fOpt(0),
be6e5811 61 fHadCorr(0x0),
62 fTpcGrid(0x0),
63 fEmcalGrid(0x0),
64 fGrid0(0),
65 fGrid1(0),
66 fGrid2(0),
67 fGrid3(0),
68 fGrid4(0),
be6e5811 69 fApplyElectronCorrection(kFALSE),
be6e5811 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),
be6e5811 80 fRefArray(0x0),
81 fProcId(kFALSE)
586f2bc3 82{
83 // Constructor
84}
85
86//____________________________________________________________________________
87
88AliJetAODReader::~AliJetAODReader()
89{
90 // Destructor
586f2bc3 91 delete fAOD;
92 delete fRef;
be6e5811 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
586f2bc3 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)){
225f0094 128 fChain->AddFile(Form("%s/%s/aod.root",dirName,name));
586f2bc3 129 a++;
130 }
131 }
132
133 gSystem->FreeDirectory(dir);
134
586f2bc3 135 fAOD = 0;
136 fChain->SetBranchAddress("AOD",&fAOD);
137
138 int nMax = fChain->GetEntries();
139
be6e5811 140 printf("\n AliJetAODReader: Total number of events in chain= %d \n",nMax);
586f2bc3 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();
be6e5811 159 printf("\n AliJetAODReader: Total number of events in chain= %5d \n", nMax);
586f2bc3 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
0e4398c8 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
f7ce54e3 240 if (mcTrack == 0){
241 fRef->Delete(); // make sure to delete before placement new...
242 new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
243 }
0e4398c8 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
586f2bc3 261//____________________________________________________________________________
262
9e4cc50d 263Bool_t AliJetAODReader::FillMomentumArray()
586f2bc3 264{
265 // Clear momentum array
266 ClearArray();
e53baffe 267 fRef->Clear();
586f2bc3 268 fDebug = fReaderHeader->GetDebug();
0e4398c8 269
586f2bc3 270 if (!fAOD) {
271 return kFALSE;
272 }
f7ce54e3 273
586f2bc3 274 // get cuts set by user
275 Float_t ptMin = fReaderHeader->GetPtCut();
276 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
277 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
d61f057d 278 UInt_t filterMask = ((AliJetAODReaderHeader*)fReaderHeader)->GetTestFilterMask();
e53baffe 279
d73efa93 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
586f2bc3 314 Int_t aodTrack = 0;
315 Float_t pt, eta;
316 TVector3 p3;
d61f057d 317
d73efa93 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
d61f057d 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();
d73efa93 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
e53baffe 386
6473bac6 387 if (aodTrack == 0){
388 fRef->Delete(); // make sure to delete before placement new...
389 new(fRef) TRefArray(TProcessID::GetProcessWithUID(track));
390 }
e53baffe 391 new ((*fMomentumArray)[aodTrack]) TLorentzVector(p3,p3.Mag());
392 sflag[aodTrack] = (TMath::Abs(track->GetLabel()) < 10000) ? 1 : 0;
d73efa93 393 cflag[aodTrack] = ( pt > ptMin ) ? 1 : 0;
e53baffe 394 aodTrack++;
d61f057d 395 fRef->Add(track);
d73efa93 396 // printf("added non-standard track\n");
586f2bc3 397 }
d73efa93 398
586f2bc3 399 // set the signal flags
400 fSignalFlag.Set(aodTrack,sflag);
401 fCutFlag.Set(aodTrack,cflag);
d61f057d 402
680ed75f 403 delete [] sflag;
404 delete [] cflag;
d73efa93 405
586f2bc3 406 return kTRUE;
407}
be6e5811 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{
4751efb5 473 //
474 // For reader task initialization
475 //
476
be6e5811 477 fDebug = fReaderHeader->GetDebug();
478 fDZ = fReaderHeader->GetDZ();
479 fTree = tree;
480
f45e6f72 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
be6e5811 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);
e36a3f22 514 // Add the task to global task
be6e5811 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();
e36a3f22 540 fRef->Clear();
be6e5811 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);
e36a3f22 555 fFillUAFromTracks->SetReferences(fRef);
556 fFillUAFromTracks->SetSignalFlag(fSignalFlag);
557 fFillUAFromTracks->SetCutFlag(fCutFlag);
be6e5811 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 }
e36a3f22 565 fSignalFlag = fFillUAFromTracks->GetSignalFlag();
566 fCutFlag = fFillUAFromTracks->GetCutFlag();
be6e5811 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
be6e5811 588
be6e5811 589
590//____________________________________________________________________________
591void AliJetAODReader::InitParameters()
592{
593 // Initialise parameters
594 fOpt = fReaderHeader->GetDetector();
f45e6f72 595 fNumUnits = fGeom->GetEMCGeometry()->GetNCells(); // Number of cells in EMCAL
be6e5811 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();
f23f7908 604 Float_t eta = 0., phi = 0., deltaEta = 0., deltaPhi = 0.;
be6e5811 605 if(fArrayInitialised) fUnitArray->Delete();
606
607 if(fTpcGrid->GetGridType()==0)
608 { // Fill the following quantities :
f23f7908 609 // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, deltaEta, deltaPhi,
be6e5811 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);
f23f7908 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);
be6e5811 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.;
f45e6f72 632 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
be6e5811 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();
f45e6f72 638 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
be6e5811 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();
f45e6f72 644 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
be6e5811 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();
f45e6f72 650 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
be6e5811 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();
f45e6f72 656 fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
be6e5811 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);
f23f7908 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);
be6e5811 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);
f23f7908 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);
be6e5811 685 }
686 else {
687 if(fDZ) {
688 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
689 if(nBin<fNumUnits+nElements+n0)
690 {
f23f7908 691 phi = eta = 0.;
be6e5811 692 fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi,eta);
f23f7908 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);
be6e5811 696 }
697 else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
698 {
f23f7908 699 phi = eta = 0.;
be6e5811 700 fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi,eta);
f23f7908 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);
be6e5811 704 }
705 else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
706 {
f23f7908 707 phi = eta = 0.;
be6e5811 708 fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi,eta);
f23f7908 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);
be6e5811 712 }
713 else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
714 {
f23f7908 715 phi = eta = 0.;
be6e5811 716 fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi,eta);
f23f7908 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);
be6e5811 720 }
721 else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
722 {
f23f7908 723 phi = eta = 0.;
be6e5811 724 fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi,eta);
f23f7908 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);
be6e5811 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}
0e4398c8 737