Add a protection to avoid crash in QA (Julian)
[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;
aa862876 189// Temporarily commented - will be updated removing the AliAODpid dependence
190/*
e36a3f22 191 if (fApplyFractionHadronicCorrection) {
e36a3f22 192 TArrayS* matched = new TArrayS();
193 GetTracksPointingToCell(matched, etaD, phiD,0.02);
194
195 double ptot = 0;
196 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
197 {
08c28025 198 const short indexS = matched->At(itrack);
199 if (indexS>-1)
e36a3f22 200 {
08c28025 201 AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22 202 ptot += track->P();
e36a3f22 203 }
204 }
205
206 correction = ptot * fFractionHadronicCorrection;
207
208 if (ptot>0 && fDebug>1 ) {
209 printf("SUM of track momentum for this cell: p=%f \n", ptot);
210 printf("fractional correction=%f \n", fFractionHadronicCorrection);
211 printf("TOTAL correction=%f \n", correction );
212 }
213
214 delete matched;
215
216 }//end hadronic correction
aa862876 217*/
e36a3f22 218 double enerCorr = digitEn;
219 if (correction >= enerCorr) enerCorr = 0;
220 else enerCorr -= correction;
221 if (correction>0 && fDebug>1) {
222 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
223 }
224
225 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
226
be6e5811 227 // Detector flag
228 if(unitEnergy>0.)
229 uArray->SetUnitDetectorFlag(kAll);
230 else uArray->SetUnitDetectorFlag(kEmcal);
231
232 uArray->SetUnitEnergy(unitEnergy+digitEt);
233
234 uArray->SetUnitCutFlag(kPtHigher);
235
236 // To be modified !!!
237 uArray->SetUnitSignalFlag(kGood);
238
239 // This is for jet multiplicity
240 uArray->SetUnitClusterID(index);
241
242 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
243
244 } // End loop over cells
245 } // end if !fCluster
246 else { // Keep digit information from clusterization
247
248 // Loop over calo clusters
249 //------------------------------------------------------------------
250
251 // select EMCAL clusters only
252 TRefArray * caloClusters = new TRefArray();
253 fAOD->GetEMCALClusters(caloClusters);
254
255 // Total number of EMCAL cluster
256 Int_t nclus = caloClusters->GetEntries() ;
257 Int_t beg = 0;
258 Float_t pos[3] ;
259
260 // Get CaloCells
261 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
be6e5811 262
263 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
264 // Retrieve cluster from aod
08c28025 265 fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
be6e5811 266
267 // Get the cluster info
be6e5811 268
269 fClus->GetPosition(pos) ;
270 TVector3 vpos(pos[0],pos[1],pos[2]) ;
271
272 Int_t digMult = fClus->GetNCells() ;
273 UShort_t *digID = fClus->GetCellsAbsId() ;
be6e5811 274 Int_t trackIndex = -1;
275 if(fClus->GetNTracksMatched()!=0)
276 trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
e36a3f22 277
be6e5811 278 // Do double-counted electron correction
279 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
280 {
281 // The electron correction go there
282 // Under construction !!!!
283 } // End of Electron correction
284
285 // Get CaloCells of cluster and fill the unitArray
286 for(Int_t i = 0; i < digMult ; i++) {
287 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
be6e5811 288 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
289
290 // Calibration for an energy in GeV
291 Float_t digitEn = digitAmp*0.0153;
292
293 Float_t etaD=-10., phiD=-10.;
294 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
295 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
296
297 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
298
e36a3f22 299 // Hadronic Correction
300 double correction=0;
aa862876 301// Temporarily commented - will be updated removing the AliAODpid dependence
302/*
e36a3f22 303 if (fApplyFractionHadronicCorrection) {
e36a3f22 304 TArrayS* matched = new TArrayS();
305 GetTracksPointingToCell(matched, etaD, phiD,0.02);
306
307 double ptot = 0;
308 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
309 {
08c28025 310 const short indexS = matched->At(itrack);
311 if (indexS>-1)
e36a3f22 312 {
08c28025 313 AliAODTrack *track = fAOD->GetTrack(indexS);
e36a3f22 314 ptot += track->P();
e36a3f22 315 }
316 }
317
318 correction = ptot * fFractionHadronicCorrection;
319
320 if (ptot>0 && fDebug>1 ) {
321 printf("SUM of track momentum for this cell: p=%f \n", ptot);
322 printf("fractional correction=%f \n", fFractionHadronicCorrection);
323 printf("TOTAL correction=%f \n", correction );
324 }
325
326 delete matched;
327
328 }//end hadronic correction
aa862876 329*/
e36a3f22 330 double enerCorr = digitEn;
331 if (correction >= enerCorr) enerCorr = 0;
332 else enerCorr -= correction;
333 if (correction>0 && fDebug>1) {
334 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
335 }
336
337 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
be6e5811 338
339 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
340 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
341 uArray->SetUnitTrackID(digitID);
342
343 Float_t unitEnergy = 0.;
344 Bool_t ok = kFALSE;
345 unitEnergy = uArray->GetUnitEnergy();
346
347 if(unitEnergy==0){
348 if(!fProcId){
349 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
350 fProcId = kTRUE;
351 }
352 fRefArray->Add(uArray);
353 fNDigitEmcal++;
354 ok = kTRUE;
355 }
356
357 // Detector flag
358 if(unitEnergy>0.)
359 uArray->SetUnitDetectorFlag(kAll);
360 else uArray->SetUnitDetectorFlag(kEmcal);
361
362 uArray->SetUnitEnergy(unitEnergy+digitEt);
363
364 uArray->SetUnitCutFlag(kPtHigher);
365
366 // Signal or background
367 // To be modified !!!
368 uArray->SetUnitSignalFlag(kGood);
369
370 // This is for jet multiplicity
371 uArray->SetUnitClusterID(index);
372
373 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
374
375 }// End loop over cells
376 } // End loop over clusters
377 } // end else
378
379 fNIn += goodDigit;
380
381 if(fDebug>1)
382 {
383 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
384 printf("goodDigit : %d\n", goodDigit);
385 }
386}
387
aa862876 388////____________________________________________________________________________
389//void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
390//{ // Temporarily commented -> will be corrected removing the dependence to AliAODPid
391//// Get all tracks pointing to cell
392// int size=0;
393//
394// for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop
395//
396// AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
397// AliAODPid* pid = (AliAODPid*) track->GetDetPid();
398//
399// if(pid) {
400// Double_t emcpos[3];
401// pid->GetEMCALPosition(emcpos);
402// TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
403//
404// Double_t deta = tpos.Eta() - eta;
405// Double_t dphi = tpos.Phi() - phi;
406//
407// if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
408// if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
409//
410// Double_t res = sqrt(dphi*dphi + deta*deta);
411//
412// if (res< cut) {
413// //add this track-index
414// size++;
415// array->Set( size );
416// array->AddAt( itrk, (size-1) );
417// if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
418//
419// }
420// }
421// }
422//
423//}
e36a3f22 424
be6e5811 425/*
426//____________________________________________________________________________
427void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
428{
429 // Set calibration parameters:
430 // if calibration database exists, they are read from database,
431 // otherwise, they are taken from digitizer.
432 //
433 // It is a user responsilibity to open CDB before reconstruction,
434 // for example:
435 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
436
437 //Check if calibration is stored in data base
438
439 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
440 {
441 AliCDBEntry *entry = (AliCDBEntry*)
442 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
443 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
444 }
445
446 if(!fCalibData)
447 printf("************* Calibration parameters not found in CDB! ****************");
448// AliFatal("Calibration parameters not found in CDB!");
449
450
451}
452
453//____________________________________________________________________________
454Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
455{
456
457 // Convert digitized amplitude into energy.
458 // Calibration parameters are taken from calibration data base for raw data,
459 // or from digitizer parameters for simulated data.
460
461 if(fCalibData){
462
463 if (fGeom==0)
464 printf("************* Did not get geometry from EMCALLoader ***************");
465// AliFatal("Did not get geometry from EMCALLoader") ;
466
467 Int_t iSupMod = -1;
468 Int_t nModule = -1;
469 Int_t nIphi = -1;
470 Int_t nIeta = -1;
471 Int_t iphi = -1;
472 Int_t ieta = -1;
473
474 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
475 if(!bCell) {
476 // fGeom->PrintGeometry();
477 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
478 assert(0);
479 }
480
481 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
482
483 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
484 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
485
486 return -fADCpedestalECA + amp * fADCchannelECA ;
487
488 }
489 else //Return energy with default parameters if calibration is not available
490 return -fADCpedestalECA + amp * fADCchannelECA ;
491
492}
493*/
494
495