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 AOD reader for jet analysis
20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 //======================================================================
28 // --- AliRoot header files ---
29 #include "AliJetUnitArray.h"
30 #include "AliJetAODFillUnitArrayEMCalDigits.h"
33 // #include "AliEMCALCalibData.h"
34 // #include "AliCDBManager.h"
35 // class AliCDBStorage;
36 // #include "AliCDBEntry.h"
39 ClassImp(AliJetAODFillUnitArrayEMCalDigits)
41 //_____________________________________________________________________________
42 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits()
43 : AliJetFillUnitArray(),
50 fApplyElectronCorrection(kFALSE),
51 fApplyFractionHadronicCorrection(kFALSE),
52 fFractionHadronicCorrection(0.3),
60 //_____________________________________________________________________________
61 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod)
62 : AliJetFillUnitArray(),
69 fApplyElectronCorrection(kFALSE),
70 fApplyFractionHadronicCorrection(kFALSE),
71 fFractionHadronicCorrection(0.3),
79 //_____________________________________________________________________________
80 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det)
81 : AliJetFillUnitArray(det),
84 fCluster(det.fCluster),
85 fNCEMCAL(det.fNCEMCAL),
88 fApplyElectronCorrection(det.fApplyElectronCorrection),
89 fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
90 fFractionHadronicCorrection(det.fFractionHadronicCorrection),
92 fNDigitEmcal(det.fNDigitEmcal),
93 fNDigitEmcalCut(det.fNDigitEmcalCut)
98 //_____________________________________________________________________________
99 AliJetAODFillUnitArrayEMCalDigits& AliJetAODFillUnitArrayEMCalDigits::operator=(const AliJetAODFillUnitArrayEMCalDigits& other)
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;
113 fNDigitEmcal = other.fNDigitEmcal;
114 fNDigitEmcalCut = other.fNDigitEmcalCut;
120 //_____________________________________________________________________________
121 AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits()
128 //_____________________________________________________________________________
129 void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
130 //void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* option)
134 // Fill the unit array with the neutral particle information from the EMCal cells in AOD
136 fDebug = fReaderHeader->GetDebug();
137 fOpt = fReaderHeader->GetDetector();
138 fCluster = fReaderHeader->GetCluster();
143 //(not used ?) Int_t nDigitTot = 0;
145 //(not used ?) Int_t beg = 0;
146 //(not used ?) Int_t end = 0;
148 //(not used ?) Int_t count = 0;
150 //(not used ?) Int_t nRefEnt = fRefArray->GetEntries();
152 if(!fCluster) { // Keep all digit information
153 // Loop over all cell information
154 //------------------------------------------------------------------
155 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
156 Int_t ncell = cells.GetNumberOfCells() ;
157 //(not used ?) Int_t type = cells.GetType();
159 for (Int_t icell= 0; icell < ncell; icell++) {
160 Int_t digitID = cells.GetCellNumber(icell);
161 Float_t digitAmp = cells.GetAmplitude(icell);
162 Float_t digitEn = digitAmp*0.0153; // Last correct
163 // Float_t digitEn = Calibrate(digitAmp,digitID);
165 Float_t etaD=-10., phiD=-10.;
166 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
167 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
169 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
171 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
173 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
175 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
176 uArray->SetUnitTrackID(digitID);
178 Float_t unitEnergy = 0.;
180 unitEnergy = uArray->GetUnitEnergy();
184 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
187 fRefArray->Add(uArray);
194 uArray->SetUnitDetectorFlag(kAll);
195 else uArray->SetUnitDetectorFlag(kEmcal);
197 uArray->SetUnitEnergy(unitEnergy+digitEt);
199 uArray->SetUnitCutFlag(kPtHigher);
201 // To be modified !!!
202 uArray->SetUnitSignalFlag(kGood);
204 // This is for jet multiplicity
205 uArray->SetUnitClusterID(index);
207 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
209 } // End loop over cells
210 } // end if !fCluster
211 else { // Keep digit information from clusterization
213 // Loop over calo clusters
214 //------------------------------------------------------------------
216 // select EMCAL clusters only
217 TRefArray * caloClusters = new TRefArray();
218 fAOD->GetEMCALClusters(caloClusters);
220 // Total number of EMCAL cluster
221 Int_t nclus = caloClusters->GetEntries() ;
226 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
227 //(not used ?) Int_t ncell = cells.GetNumberOfCells() ;
228 //(not used ?) Int_t type = cells.GetType();
230 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
231 // Retrieve cluster from aod
232 AliAODCaloCluster *fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
234 // Get the cluster info
235 //(not used ?) Float_t energy = fClus->E() ;
237 fClus->GetPosition(pos) ;
238 TVector3 vpos(pos[0],pos[1],pos[2]) ;
240 Int_t digMult = fClus->GetNCells() ;
241 UShort_t *digID = fClus->GetCellsAbsId() ;
242 //(not used ?) Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
243 Int_t trackIndex = -1;
244 if(fClus->GetNTracksMatched()!=0)
245 trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
246 //(not used ?) Int_t cLabel = fClus->GetLabel(0);
247 // Do double-counted electron correction
248 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
250 // The electron correction go there
251 // Under construction !!!!
252 } // End of Electron correction
254 // Get CaloCells of cluster and fill the unitArray
255 for(Int_t i = 0; i < digMult ; i++) {
256 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
257 //(not used ?) Double_t digitAmpFrac = digAmpFrac[i];
258 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
260 // Calibration for an energy in GeV
261 Float_t digitEn = digitAmp*0.0153;
263 Float_t etaD=-10., phiD=-10.;
264 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
265 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
267 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
269 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
271 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
272 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
273 uArray->SetUnitTrackID(digitID);
275 Float_t unitEnergy = 0.;
277 unitEnergy = uArray->GetUnitEnergy();
281 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
284 fRefArray->Add(uArray);
291 uArray->SetUnitDetectorFlag(kAll);
292 else uArray->SetUnitDetectorFlag(kEmcal);
294 uArray->SetUnitEnergy(unitEnergy+digitEt);
296 uArray->SetUnitCutFlag(kPtHigher);
298 // Signal or background
299 // To be modified !!!
300 uArray->SetUnitSignalFlag(kGood);
302 // This is for jet multiplicity
303 uArray->SetUnitClusterID(index);
305 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
307 }// End loop over cells
308 } // End loop over clusters
315 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
316 printf("goodDigit : %d\n", goodDigit);
321 //____________________________________________________________________________
322 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
324 // Set calibration parameters:
325 // if calibration database exists, they are read from database,
326 // otherwise, they are taken from digitizer.
328 // It is a user responsilibity to open CDB before reconstruction,
330 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
332 //Check if calibration is stored in data base
334 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
336 AliCDBEntry *entry = (AliCDBEntry*)
337 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
338 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
342 printf("************* Calibration parameters not found in CDB! ****************");
343 // AliFatal("Calibration parameters not found in CDB!");
348 //____________________________________________________________________________
349 Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
352 // Convert digitized amplitude into energy.
353 // Calibration parameters are taken from calibration data base for raw data,
354 // or from digitizer parameters for simulated data.
359 printf("************* Did not get geometry from EMCALLoader ***************");
360 // AliFatal("Did not get geometry from EMCALLoader") ;
369 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
371 // fGeom->PrintGeometry();
372 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
376 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
378 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
379 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
381 return -fADCpedestalECA + amp * fADCchannelECA ;
384 else //Return energy with default parameters if calibration is not available
385 return -fADCpedestalECA + amp * fADCchannelECA ;