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