2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //---------------------------------------------------------------------
19 // Called by ESD reader for jet analysis
20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 //---------------------------------------------------------------------
23 // --- ROOT system ---
24 #include <TLorentzVector.h>
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"
37 // --- ROOT system ---
42 // --- AliRoot header files ---
45 class AliJetESDReader;
46 // class AliEMCALCalibData;
47 // class AliCDBManager;
48 // class AliCDBStorage;
53 ClassImp(AliJetESDFillUnitArrayEMCalDigits)
55 //_____________________________________________________________________________
56 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
57 : AliJetFillUnitArray(),
66 fApplyElectronCorrection(kFALSE),
67 fApplyFractionHadronicCorrection(kFALSE),
68 fFractionHadronicCorrection(0.3),
76 //_____________________________________________________________________________
77 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
78 : AliJetFillUnitArray(),
87 fApplyElectronCorrection(kFALSE),
88 fApplyFractionHadronicCorrection(kFALSE),
89 fFractionHadronicCorrection(0.3),
97 //_____________________________________________________________________________
98 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
99 : AliJetFillUnitArray(det),
103 fCluster(det.fCluster),
105 fNCEMCAL(det.fNCEMCAL),
106 fNCPHOS(det.fNCPHOS),
107 fNCCalo(det.fNCCalo),
108 fApplyElectronCorrection(det.fApplyElectronCorrection),
109 fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
110 fFractionHadronicCorrection(det.fFractionHadronicCorrection),
112 fNDigitEmcal(det.fNDigitEmcal),
113 fNDigitEmcalCut(det.fNDigitEmcalCut)
118 //_____________________________________________________________________________
119 AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
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;
135 fNDigitEmcal = other.fNDigitEmcal;
136 fNDigitEmcalCut = other.fNDigitEmcalCut;
142 //_____________________________________________________________________________
143 AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
148 //_____________________________________________________________________________
149 void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
153 // Fill the unit array with the neutral particle information from the EMCal cells in ESD
156 fDebug = fReaderHeader->GetDebug();
157 fOpt = fReaderHeader->GetDetector();
158 fCluster = fReaderHeader->GetCluster();
163 //(not used ?) Int_t nDigitTot = 0;
165 //(not used ?) Int_t beg = 0;
166 //(not used ?) Int_t end = 0;
168 //(not used ?) Int_t count = 0;
170 //(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
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();
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);
186 Float_t etaD=-10., phiD=-10.;
187 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
188 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
190 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
192 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
194 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
196 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
197 uArray->SetUnitTrackID(digitID);
199 Float_t unitEnergy = 0.;
201 unitEnergy = uArray->GetUnitEnergy();
205 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
208 fRefArray->Add(uArray);
215 uArray->SetUnitDetectorFlag(kAll);
216 else uArray->SetUnitDetectorFlag(kEmcal);
218 uArray->SetUnitEnergy(unitEnergy+digitEt);
220 uArray->SetUnitCutFlag(kPtHigher);
222 // To be modified !!!
223 uArray->SetUnitSignalFlag(kGood);
225 // This is for jet multiplicity
226 uArray->SetUnitClusterID(index);
228 if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
230 } // End loop over cells
231 } // end if !fCluster
232 else { // Keep digit information from clusterization
234 // Loop over calo clusters
235 //------------------------------------------------------------------
237 //select EMCAL clusters only
238 TRefArray * caloClusters = new TRefArray();
239 fESD->GetEMCALClusters(caloClusters);
241 // Total number of EMCAL cluster
242 Int_t nclus = caloClusters->GetEntries() ;
246 // Get reconstructed vertex position
247 Double_t vertexPosition[3] ;
248 fESD->GetVertex()->GetXYZ(vertexPosition) ;
251 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
252 //(not used ?) Int_t ncell = cells.GetNumberOfCells() ;
253 //(not used ?) Int_t type = cells.GetType();
255 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
256 // Retrieve cluster from esd
257 fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
259 // Get the cluster info
260 //(not used ?) Float_t energy = fClus->E() ;
261 //(not used ?) Int_t iprim = fClus->GetLabel();
263 fClus->GetPosition(pos) ;
264 TVector3 vpos(pos[0],pos[1],pos[2]) ;
266 fClus->GetMomentum(p,vertexPosition);
268 Int_t digMult = fClus->GetNCells() ;
269 UShort_t *digID = fClus->GetCellsAbsId() ;
270 //(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
271 Int_t trackIndex = fClus->GetTrackMatchedIndex();
273 // Do double-counted electron correction
274 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
276 // The electron correction go there
277 // Under construction !!!!
278 } // End of Electron correction
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);
287 // Calibration for an energy in GeV
288 Float_t digitEn = digitAmp*0.0153;
290 Float_t etaD=-10., phiD=-10.;
291 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
292 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
294 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
296 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
298 cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
300 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
301 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
302 uArray->SetUnitTrackID(digitID);
304 Float_t unitEnergy = 0.;
306 unitEnergy = uArray->GetUnitEnergy();
310 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
313 fRefArray->Add(uArray);
320 uArray->SetUnitDetectorFlag(kAll);
321 else uArray->SetUnitDetectorFlag(kEmcal);
323 uArray->SetUnitEnergy(unitEnergy+digitEt);
325 uArray->SetUnitCutFlag(kPtHigher);
327 // To be modified !!!
328 uArray->SetUnitSignalFlag(kGood);
330 // This is for jet multiplicity
331 uArray->SetUnitClusterID(index);
333 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
335 } // End loop over cells
336 } // End loop over clusters
343 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
344 printf("goodDigit : %d\n", goodDigit);
349 //____________________________________________________________________________
350 void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
352 // Set calibration parameters:
353 // if calibration database exists, they are read from database,
354 // otherwise, they are taken from digitizer.
356 // It is a user responsilibity to open CDB before reconstruction,
358 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
360 //Check if calibration is stored in data base
362 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
364 AliCDBEntry *entry = (AliCDBEntry*)
365 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
366 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
370 printf("************* Calibration parameters not found in CDB! ****************");
371 // AliFatal("Calibration parameters not found in CDB!");
376 //____________________________________________________________________________
377 Float_t AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
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.
387 printf("************* Did not get geometry from EMCALLoader ***************");
388 // AliFatal("Did not get geometry from EMCALLoader") ;
397 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
399 // fGeom->PrintGeometry();
400 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
404 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
406 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
407 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
409 return -fADCpedestalECA + amp * fADCchannelECA ;
412 else //Return energy with default parameters if calibration is not available
413 return -fADCpedestalECA + amp * fADCchannelECA ;