Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / AliJetESDFillUnitArrayEMCalDigits.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 ESD reader for jet analysis
20// Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21//---------------------------------------------------------------------
22
23// --- ROOT system ---
24#include <TLorentzVector.h>
25#include <TMath.h>
26
27// --- AliRoot header files ---
28#include "AliESDCaloCluster.h"
29#include "AliESDCaloCells.h"
30#include "AliJetUnitArray.h"
31#include "AliJetESDFillUnitArrayEMCalDigits.h"
32// #include "AliEMCALCalibData.h"
33// #include "AliCDBManager.h"
34// // class AliCDBStorage;
35// #include "AliCDBEntry.h"
36
37// --- ROOT system ---
38class TSystem;
39class TGeoManager;
40class TArrayS;
41
42// --- AliRoot header files ---
43class AliJetFinder;
44class AliJetReader;
45class AliJetESDReader;
46// class AliEMCALCalibData;
47// class AliCDBManager;
48// class AliCDBStorage;
49// class AliCDBEntry;
50
3a7af7bd 51using std::cout;
52using std::endl;
be6e5811 53ClassImp(AliJetESDFillUnitArrayEMCalDigits)
54
55//_____________________________________________________________________________
56AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
57 : AliJetFillUnitArray(),
58 fESD(0),
59 fNIn(0),
60 fOpt(0),
61 fCluster(0),
62 fDebug(0),
63 fNCEMCAL(0),
64 fNCPHOS(0),
65 fNCCalo(0),
66 fApplyElectronCorrection(kFALSE),
67 fApplyFractionHadronicCorrection(kFALSE),
68 fFractionHadronicCorrection(0.3),
69 fClus(0x0),
70 fNDigitEmcal(0),
71 fNDigitEmcalCut(0)
72{
73 // constructor
74}
75
76//_____________________________________________________________________________
77AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
78 : AliJetFillUnitArray(),
79 fESD(esd),
80 fNIn(0),
81 fOpt(0),
82 fCluster(0),
83 fDebug(0),
84 fNCEMCAL(0),
85 fNCPHOS(0),
86 fNCCalo(0),
87 fApplyElectronCorrection(kFALSE),
88 fApplyFractionHadronicCorrection(kFALSE),
89 fFractionHadronicCorrection(0.3),
90 fClus(0x0),
91 fNDigitEmcal(0),
92 fNDigitEmcalCut(0)
93{
94 // constructor
95}
96
97//_____________________________________________________________________________
98AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
99 : AliJetFillUnitArray(det),
100 fESD(det.fESD),
101 fNIn(det.fNIn),
102 fOpt(det.fOpt),
103 fCluster(det.fCluster),
104 fDebug(det.fDebug),
105 fNCEMCAL(det.fNCEMCAL),
106 fNCPHOS(det.fNCPHOS),
107 fNCCalo(det.fNCCalo),
108 fApplyElectronCorrection(det.fApplyElectronCorrection),
109 fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
110 fFractionHadronicCorrection(det.fFractionHadronicCorrection),
111 fClus(det.fClus),
112 fNDigitEmcal(det.fNDigitEmcal),
113 fNDigitEmcalCut(det.fNDigitEmcalCut)
114{
115 // Copy constructor
116}
117
118//_____________________________________________________________________________
119AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
120{
121 // Assignment
122
123 fESD = other.fESD;
124 fNIn = other.fNIn;
125 fOpt = other.fOpt;
126 fCluster = other.fCluster;
127 fDebug = other.fDebug;
128 fNCEMCAL = other.fNCEMCAL;
129 fNCPHOS = other.fNCPHOS;
130 fNCCalo = other.fNCCalo;
131 fApplyElectronCorrection = other.fApplyElectronCorrection;
132 fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
133 fFractionHadronicCorrection = other.fFractionHadronicCorrection;
134 fClus = other.fClus;
135 fNDigitEmcal = other.fNDigitEmcal;
136 fNDigitEmcalCut = other.fNDigitEmcalCut;
137
138 return (*this);
139
140}
141
142//_____________________________________________________________________________
143AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
144{
145 // destructor
146}
147
148//_____________________________________________________________________________
149void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
150{
151 //
152 // Main method.
153 // Fill the unit array with the neutral particle information from the EMCal cells in ESD
154 //
155
156 fDebug = fReaderHeader->GetDebug();
157 fOpt = fReaderHeader->GetDetector();
158 fCluster = fReaderHeader->GetCluster();
159
160 // Init parameters
161 InitParameters();
162
163//(not used ?) Int_t nDigitTot = 0;
164 Int_t goodDigit = 0;
165//(not used ?) Int_t beg = 0;
166//(not used ?) Int_t end = 0;
167 Int_t index = 0;
168//(not used ?) Int_t count = 0;
169
170//(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
171
172 if(!fCluster) { // Keep all digit information
173 // Loop over all cell information
174 //------------------------------------------------------------------
175 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
176 Int_t ncell = cells.GetNumberOfCells() ;
177//(not used ?) Int_t type = cells.GetType();
178
179 for (Int_t icell= 0; icell < ncell; icell++) {
180 Int_t digitID = cells.GetCellNumber(icell);
181 Float_t digitAmp = cells.GetAmplitude(icell);
182//(not used ?) Float_t digitTime = cells.GetTime(icell);
183 Float_t digitEn = digitAmp*0.0153; // Last correct
184 // Float_t digitEn = Calibrate(digitAmp,digitID);
185
186 Float_t etaD=-10., phiD=-10.;
187 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
188 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
189
190 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
191
192 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
193
194 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
195
196 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
197 uArray->SetUnitTrackID(digitID);
198
199 Float_t unitEnergy = 0.;
200 Bool_t ok = kFALSE;
201 unitEnergy = uArray->GetUnitEnergy();
202
203 if(unitEnergy==0){
204 if(!fProcId){
205 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
206 fProcId = kTRUE;
207 }
208 fRefArray->Add(uArray);
209 fNDigitEmcal++;
210 ok = kTRUE;
211 }
212
213 // Detector flag
214 if(unitEnergy>0.)
215 uArray->SetUnitDetectorFlag(kAll);
216 else uArray->SetUnitDetectorFlag(kEmcal);
217
218 uArray->SetUnitEnergy(unitEnergy+digitEt);
219
220 uArray->SetUnitCutFlag(kPtHigher);
221
222 // To be modified !!!
223 uArray->SetUnitSignalFlag(kGood);
224
225 // This is for jet multiplicity
226 uArray->SetUnitClusterID(index);
227
228 if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
229
230 } // End loop over cells
231 } // end if !fCluster
232 else { // Keep digit information from clusterization
233
234 // Loop over calo clusters
235 //------------------------------------------------------------------
236
237 //select EMCAL clusters only
238 TRefArray * caloClusters = new TRefArray();
239 fESD->GetEMCALClusters(caloClusters);
240
241 // Total number of EMCAL cluster
242 Int_t nclus = caloClusters->GetEntries() ;
243 Int_t beg = 0;
244 Float_t pos[3] ;
245
246 // Get reconstructed vertex position
247 Double_t vertexPosition[3] ;
248 fESD->GetVertex()->GetXYZ(vertexPosition) ;
249
250 // Get CaloCells
251 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
252//(not used ?) Int_t ncell = cells.GetNumberOfCells() ;
253//(not used ?) Int_t type = cells.GetType();
254
255 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
256 // Retrieve cluster from esd
08c28025 257 fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
be6e5811 258
259 // Get the cluster info
260//(not used ?) Float_t energy = fClus->E() ;
261//(not used ?) Int_t iprim = fClus->GetLabel();
262
263 fClus->GetPosition(pos) ;
264 TVector3 vpos(pos[0],pos[1],pos[2]) ;
265 TLorentzVector p ;
266 fClus->GetMomentum(p,vertexPosition);
267
268 Int_t digMult = fClus->GetNCells() ;
269 UShort_t *digID = fClus->GetCellsAbsId() ;
270//(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
c8fe2783 271 Int_t trackIndex = fClus->GetTrackMatchedIndex();
be6e5811 272
273 // Do double-counted electron correction
274 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
275 {
276 // The electron correction go there
277 // Under construction !!!!
278 } // End of Electron correction
279
280 // Get CaloCells of cluster and fill the unitArray
281 for(Int_t i = 0; i < digMult ; i++){
282 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
283//(not used ?) Double_t digitAmpFrac = digAmpFrac[i];
284 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
285//(not used ?) Float_t digitTime = cells.GetCellTime(digitID);
286
287 // Calibration for an energy in GeV
288 Float_t digitEn = digitAmp*0.0153;
289
290 Float_t etaD=-10., phiD=-10.;
291 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
292 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
293
294 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
295
296 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
297
298 cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
299
300 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
301 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
302 uArray->SetUnitTrackID(digitID);
303
304 Float_t unitEnergy = 0.;
305 Bool_t ok = kFALSE;
306 unitEnergy = uArray->GetUnitEnergy();
307
308 if(unitEnergy==0){
309 if(!fProcId){
310 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
311 fProcId = kTRUE;
312 }
313 fRefArray->Add(uArray);
314 fNDigitEmcal++;
315 ok = kTRUE;
316 }
317
318 // Detector flag
319 if(unitEnergy>0.)
320 uArray->SetUnitDetectorFlag(kAll);
321 else uArray->SetUnitDetectorFlag(kEmcal);
322
323 uArray->SetUnitEnergy(unitEnergy+digitEt);
324
325 uArray->SetUnitCutFlag(kPtHigher);
326
327 // To be modified !!!
328 uArray->SetUnitSignalFlag(kGood);
329
330 // This is for jet multiplicity
331 uArray->SetUnitClusterID(index);
332
333 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
334
335 } // End loop over cells
336 } // End loop over clusters
337 } // end else
338
339 fNIn += goodDigit;
340
341 if(fDebug>1)
342 {
343 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
344 printf("goodDigit : %d\n", goodDigit);
345 }
346}
347
348/*
349//____________________________________________________________________________
350void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
351{
352 // Set calibration parameters:
353 // if calibration database exists, they are read from database,
354 // otherwise, they are taken from digitizer.
355 //
356 // It is a user responsilibity to open CDB before reconstruction,
357 // for example:
358 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
359
360 //Check if calibration is stored in data base
361
362 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
363 {
364 AliCDBEntry *entry = (AliCDBEntry*)
365 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
366 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
367 }
368
369 if(!fCalibData)
370 printf("************* Calibration parameters not found in CDB! ****************");
371// AliFatal("Calibration parameters not found in CDB!");
372
373
374}
375
376//____________________________________________________________________________
377Float_t AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
378{
379
380 // Convert digitized amplitude into energy.
381 // Calibration parameters are taken from calibration data base for raw data,
382 // or from digitizer parameters for simulated data.
383
384 if(fCalibData){
385
386 if (fGeom==0)
387 printf("************* Did not get geometry from EMCALLoader ***************");
388// AliFatal("Did not get geometry from EMCALLoader") ;
389
390 Int_t iSupMod = -1;
391 Int_t nModule = -1;
392 Int_t nIphi = -1;
393 Int_t nIeta = -1;
394 Int_t iphi = -1;
395 Int_t ieta = -1;
396
397 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
398 if(!bCell) {
399 // fGeom->PrintGeometry();
400 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
401 assert(0);
402 }
403
404 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
405
406 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
407 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
408
409 return -fADCpedestalECA + amp * fADCchannelECA ;
410
411 }
412 else //Return energy with default parameters if calibration is not available
413 return -fADCpedestalECA + amp * fADCchannelECA ;
414
415}
416*/
417