]>
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 AOD reader for jet analysis | |
20 | // Author: Magali Estienne (magali.estienne@ires.in2p3.fr) | |
21 | //====================================================================== | |
22 | ||
23 | class TSystem; | |
24 | class TLorentzVector; | |
25 | class TGeoManager; | |
26 | class TArrayS; | |
27 | ||
28 | // --- AliRoot header files --- | |
29 | #include "AliJetUnitArray.h" | |
30 | #include "AliJetAODFillUnitArrayEMCalDigits.h" | |
31 | class AliJetFinder; | |
32 | ||
33 | // #include "AliEMCALCalibData.h" | |
34 | // #include "AliCDBManager.h" | |
35 | // class AliCDBStorage; | |
36 | // #include "AliCDBEntry.h" | |
37 | ||
38 | ||
39 | ClassImp(AliJetAODFillUnitArrayEMCalDigits) | |
40 | ||
41 | //_____________________________________________________________________________ | |
42 | AliJetAODFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
61 | AliJetAODFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
80 | AliJetAODFillUnitArrayEMCalDigits::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 | //_____________________________________________________________________________ | |
99 | AliJetAODFillUnitArrayEMCalDigits& 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 | //_____________________________________________________________________________ | |
121 | AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits() | |
122 | { | |
123 | // destructor | |
124 | } | |
125 | ||
126 | ||
127 | ||
128 | //_____________________________________________________________________________ | |
129 | void 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 | //____________________________________________________________________________ | |
427 | void 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 | //____________________________________________________________________________ | |
454 | Float_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 |