]>
Commit | Line | Data |
---|---|---|
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 --- | |
38 | class TSystem; | |
39 | class TGeoManager; | |
40 | class TArrayS; | |
41 | ||
42 | // --- AliRoot header files --- | |
43 | class AliJetFinder; | |
44 | class AliJetReader; | |
45 | class AliJetESDReader; | |
46 | // class AliEMCALCalibData; | |
47 | // class AliCDBManager; | |
48 | // class AliCDBStorage; | |
49 | // class AliCDBEntry; | |
50 | ||
3a7af7bd | 51 | using std::cout; |
52 | using std::endl; | |
be6e5811 | 53 | ClassImp(AliJetESDFillUnitArrayEMCalDigits) |
54 | ||
55 | //_____________________________________________________________________________ | |
56 | AliJetESDFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
77 | AliJetESDFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
98 | AliJetESDFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
119 | AliJetESDFillUnitArrayEMCalDigits& 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 | //_____________________________________________________________________________ | |
143 | AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits() | |
144 | { | |
145 | // destructor | |
146 | } | |
147 | ||
148 | //_____________________________________________________________________________ | |
149 | void 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 | //____________________________________________________________________________ | |
350 | void 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 | //____________________________________________________________________________ | |
377 | Float_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 |