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