Fixes to track selection for embedding in cluster task (Bastian), Do no select events...
[u/mrichter/AliRoot.git] / JETAN / AliJetAODFillUnitArrayEMCalDigits.cxx
CommitLineData
be6e5811 1
2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17//======================================================================
18// Fill Unit Array
19// Called by AOD reader for jet analysis
20// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21//======================================================================
22
23class TSystem;
24class TLorentzVector;
25class TGeoManager;
26class TArrayS;
27
28// --- AliRoot header files ---
29#include "AliJetUnitArray.h"
30#include "AliJetAODFillUnitArrayEMCalDigits.h"
31class AliJetFinder;
32
33// #include "AliEMCALCalibData.h"
34// #include "AliCDBManager.h"
35// class AliCDBStorage;
36// #include "AliCDBEntry.h"
37
38
39ClassImp(AliJetAODFillUnitArrayEMCalDigits)
40
41//_____________________________________________________________________________
42AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits()
43 : AliJetFillUnitArray(),
44 fAOD(0),
45 fNIn(0),
46 fCluster(0),
47 fNCEMCAL(0),
48 fNCPHOS(0),
49 fNCCalo(0),
50 fApplyElectronCorrection(kFALSE),
51 fApplyFractionHadronicCorrection(kFALSE),
52 fFractionHadronicCorrection(0.3),
53 fClus(0x0),
54 fNDigitEmcal(0),
55 fNDigitEmcalCut(0)
56{
57 // constructor
58}
59
60//_____________________________________________________________________________
61AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod)
62 : AliJetFillUnitArray(),
63 fAOD(aod),
64 fNIn(0),
65 fCluster(0),
66 fNCEMCAL(0),
67 fNCPHOS(0),
68 fNCCalo(0),
69 fApplyElectronCorrection(kFALSE),
70 fApplyFractionHadronicCorrection(kFALSE),
71 fFractionHadronicCorrection(0.3),
72 fClus(0x0),
73 fNDigitEmcal(0),
74 fNDigitEmcalCut(0)
75{
76 // constructor
77}
78
79//_____________________________________________________________________________
80AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det)
81 : AliJetFillUnitArray(det),
82 fAOD(det.fAOD),
83 fNIn(det.fNIn),
84 fCluster(det.fCluster),
85 fNCEMCAL(det.fNCEMCAL),
86 fNCPHOS(det.fNCPHOS),
87 fNCCalo(det.fNCCalo),
88 fApplyElectronCorrection(det.fApplyElectronCorrection),
89 fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
90 fFractionHadronicCorrection(det.fFractionHadronicCorrection),
91 fClus(det.fClus),
92 fNDigitEmcal(det.fNDigitEmcal),
93 fNDigitEmcalCut(det.fNDigitEmcalCut)
94{
95 // Copy constructor
96}
97
98//_____________________________________________________________________________
99AliJetAODFillUnitArrayEMCalDigits& AliJetAODFillUnitArrayEMCalDigits::operator=(const AliJetAODFillUnitArrayEMCalDigits& other)
100{
101 // Assignment
102
103 fAOD = other.fAOD;
104 fNIn = other.fNIn;
105 fCluster = other.fCluster;
106 fNCEMCAL = other.fNCEMCAL;
107 fNCPHOS = other.fNCPHOS;
108 fNCCalo = other.fNCCalo;
109 fApplyElectronCorrection = other.fApplyElectronCorrection;
110 fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
111 fFractionHadronicCorrection = other.fFractionHadronicCorrection;
112 fClus = other.fClus;
113 fNDigitEmcal = other.fNDigitEmcal;
114 fNDigitEmcalCut = other.fNDigitEmcalCut;
115
116 return (*this);
117
118}
119
120//_____________________________________________________________________________
121AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits()
122{
123 // destructor
124}
125
126
127
128//_____________________________________________________________________________
129void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
130//void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* option)
131{
132 //
133 // Main method.
134 // Fill the unit array with the neutral particle information from the EMCal cells in AOD
135
136 fDebug = fReaderHeader->GetDebug();
137 fOpt = fReaderHeader->GetDetector();
138 fCluster = fReaderHeader->GetCluster();
139
140 // Init parameters
e36a3f22 141 // InitParameters();
be6e5811 142
be6e5811 143 Int_t goodDigit = 0;
be6e5811 144 Int_t index = 0;
be6e5811 145
146 if(!fCluster) { // Keep all digit information
e36a3f22 147
be6e5811 148 // Loop over all cell information
149 //------------------------------------------------------------------
150 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
151 Int_t ncell = cells.GetNumberOfCells() ;
152//(not used ?) Int_t type = cells.GetType();
153
154 for (Int_t icell= 0; icell < ncell; icell++) {
155 Int_t digitID = cells.GetCellNumber(icell);
156 Float_t digitAmp = cells.GetAmplitude(icell);
157 Float_t digitEn = digitAmp*0.0153; // Last correct
158 // Float_t digitEn = Calibrate(digitAmp,digitID);
159
160 Float_t etaD=-10., phiD=-10.;
161 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
162 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
163
164 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
165
e36a3f22 166 // Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
be6e5811 167
168 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
169
170 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
171 uArray->SetUnitTrackID(digitID);
172
173 Float_t unitEnergy = 0.;
174 Bool_t ok = kFALSE;
175 unitEnergy = uArray->GetUnitEnergy();
176
177 if(unitEnergy==0){
178 if(!fProcId){
179 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
180 fProcId = kTRUE;
181 }
182 fRefArray->Add(uArray);
183 fNDigitEmcal++;
184 ok = kTRUE;
185 }
186
e36a3f22 187 // Hadronic Correction
188 double correction=0;
189 if (fApplyFractionHadronicCorrection) {
190
191 TArrayS* matched = new TArrayS();
192 GetTracksPointingToCell(matched, etaD, phiD,0.02);
193
194 double ptot = 0;
195 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
196 {
08c28025 197 const short indexS = matched->At(itrack);
198 if (indexS>-1)
e36a3f22 199 {
08c28025 200 AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22 201 ptot += track->P();
e36a3f22 202 }
203 }
204
205 correction = ptot * fFractionHadronicCorrection;
206
207 if (ptot>0 && fDebug>1 ) {
208 printf("SUM of track momentum for this cell: p=%f \n", ptot);
209 printf("fractional correction=%f \n", fFractionHadronicCorrection);
210 printf("TOTAL correction=%f \n", correction );
211 }
212
213 delete matched;
214
215 }//end hadronic correction
216
217 double enerCorr = digitEn;
218 if (correction >= enerCorr) enerCorr = 0;
219 else enerCorr -= correction;
220 if (correction>0 && fDebug>1) {
221 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
222 }
223
224 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
225
be6e5811 226 // Detector flag
227 if(unitEnergy>0.)
228 uArray->SetUnitDetectorFlag(kAll);
229 else uArray->SetUnitDetectorFlag(kEmcal);
230
231 uArray->SetUnitEnergy(unitEnergy+digitEt);
232
233 uArray->SetUnitCutFlag(kPtHigher);
234
235 // To be modified !!!
236 uArray->SetUnitSignalFlag(kGood);
237
238 // This is for jet multiplicity
239 uArray->SetUnitClusterID(index);
240
241 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
242
243 } // End loop over cells
244 } // end if !fCluster
245 else { // Keep digit information from clusterization
246
247 // Loop over calo clusters
248 //------------------------------------------------------------------
249
250 // select EMCAL clusters only
251 TRefArray * caloClusters = new TRefArray();
252 fAOD->GetEMCALClusters(caloClusters);
253
254 // Total number of EMCAL cluster
255 Int_t nclus = caloClusters->GetEntries() ;
256 Int_t beg = 0;
257 Float_t pos[3] ;
258
259 // Get CaloCells
260 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
be6e5811 261
262 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
263 // Retrieve cluster from aod
08c28025 264 fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
be6e5811 265
266 // Get the cluster info
be6e5811 267
268 fClus->GetPosition(pos) ;
269 TVector3 vpos(pos[0],pos[1],pos[2]) ;
270
271 Int_t digMult = fClus->GetNCells() ;
272 UShort_t *digID = fClus->GetCellsAbsId() ;
be6e5811 273 Int_t trackIndex = -1;
274 if(fClus->GetNTracksMatched()!=0)
275 trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
e36a3f22 276
be6e5811 277 // Do double-counted electron correction
278 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
279 {
280 // The electron correction go there
281 // Under construction !!!!
282 } // End of Electron correction
283
284 // Get CaloCells of cluster and fill the unitArray
285 for(Int_t i = 0; i < digMult ; i++) {
286 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
be6e5811 287 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
288
289 // Calibration for an energy in GeV
290 Float_t digitEn = digitAmp*0.0153;
291
292 Float_t etaD=-10., phiD=-10.;
293 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
294 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
295
296 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
297
e36a3f22 298 // Hadronic Correction
299 double correction=0;
300 if (fApplyFractionHadronicCorrection) {
301
302 TArrayS* matched = new TArrayS();
303 GetTracksPointingToCell(matched, etaD, phiD,0.02);
304
305 double ptot = 0;
306 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
307 {
08c28025 308 const short indexS = matched->At(itrack);
309 if (indexS>-1)
e36a3f22 310 {
08c28025 311 AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22 312 ptot += track->P();
e36a3f22 313 }
314 }
315
316 correction = ptot * fFractionHadronicCorrection;
317
318 if (ptot>0 && fDebug>1 ) {
319 printf("SUM of track momentum for this cell: p=%f \n", ptot);
320 printf("fractional correction=%f \n", fFractionHadronicCorrection);
321 printf("TOTAL correction=%f \n", correction );
322 }
323
324 delete matched;
325
326 }//end hadronic correction
327
328 double enerCorr = digitEn;
329 if (correction >= enerCorr) enerCorr = 0;
330 else enerCorr -= correction;
331 if (correction>0 && fDebug>1) {
332 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
333 }
334
335 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
be6e5811 336
337 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
338 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
339 uArray->SetUnitTrackID(digitID);
340
341 Float_t unitEnergy = 0.;
342 Bool_t ok = kFALSE;
343 unitEnergy = uArray->GetUnitEnergy();
344
345 if(unitEnergy==0){
346 if(!fProcId){
347 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
348 fProcId = kTRUE;
349 }
350 fRefArray->Add(uArray);
351 fNDigitEmcal++;
352 ok = kTRUE;
353 }
354
355 // Detector flag
356 if(unitEnergy>0.)
357 uArray->SetUnitDetectorFlag(kAll);
358 else uArray->SetUnitDetectorFlag(kEmcal);
359
360 uArray->SetUnitEnergy(unitEnergy+digitEt);
361
362 uArray->SetUnitCutFlag(kPtHigher);
363
364 // Signal or background
365 // To be modified !!!
366 uArray->SetUnitSignalFlag(kGood);
367
368 // This is for jet multiplicity
369 uArray->SetUnitClusterID(index);
370
371 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
372
373 }// End loop over cells
374 } // End loop over clusters
375 } // end else
376
377 fNIn += goodDigit;
378
379 if(fDebug>1)
380 {
381 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
382 printf("goodDigit : %d\n", goodDigit);
383 }
384}
385
e36a3f22 386//____________________________________________________________________________
387void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
388{
1240edf5 389// Get all tracks pointing to cell
e36a3f22 390 int size=0;
391
392 for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop
393
394 AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
395 AliAODPid* pid = (AliAODPid*) track->GetDetPid();
396
397 if(pid) {
398 Double_t emcpos[3];
399 pid->GetEMCALPosition(emcpos);
400 TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
401
402 Double_t deta = tpos.Eta() - eta;
403 Double_t dphi = tpos.Phi() - phi;
404
405 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
406 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
407
408 Double_t res = sqrt(dphi*dphi + deta*deta);
409
410 if (res< cut) {
411 //add this track-index
412 size++;
413 array->Set( size );
414 array->AddAt( itrk, (size-1) );
415 if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
416
417 }
418 }
419 }
420
421}
422
be6e5811 423/*
424//____________________________________________________________________________
425void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
426{
427 // Set calibration parameters:
428 // if calibration database exists, they are read from database,
429 // otherwise, they are taken from digitizer.
430 //
431 // It is a user responsilibity to open CDB before reconstruction,
432 // for example:
433 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
434
435 //Check if calibration is stored in data base
436
437 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
438 {
439 AliCDBEntry *entry = (AliCDBEntry*)
440 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
441 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
442 }
443
444 if(!fCalibData)
445 printf("************* Calibration parameters not found in CDB! ****************");
446// AliFatal("Calibration parameters not found in CDB!");
447
448
449}
450
451//____________________________________________________________________________
452Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
453{
454
455 // Convert digitized amplitude into energy.
456 // Calibration parameters are taken from calibration data base for raw data,
457 // or from digitizer parameters for simulated data.
458
459 if(fCalibData){
460
461 if (fGeom==0)
462 printf("************* Did not get geometry from EMCALLoader ***************");
463// AliFatal("Did not get geometry from EMCALLoader") ;
464
465 Int_t iSupMod = -1;
466 Int_t nModule = -1;
467 Int_t nIphi = -1;
468 Int_t nIeta = -1;
469 Int_t iphi = -1;
470 Int_t ieta = -1;
471
472 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
473 if(!bCell) {
474 // fGeom->PrintGeometry();
475 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
476 assert(0);
477 }
478
479 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
480
481 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
482 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
483
484 return -fADCpedestalECA + amp * fADCchannelECA ;
485
486 }
487 else //Return energy with default parameters if calibration is not available
488 return -fADCpedestalECA + amp * fADCchannelECA ;
489
490}
491*/
492
493