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();
146 if(!fCluster) { // Keep all digit information
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();
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);
160 Float_t etaD=-10., phiD=-10.;
161 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
162 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
164 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
166 // Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
168 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
170 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
171 uArray->SetUnitTrackID(digitID);
173 Float_t unitEnergy = 0.;
175 unitEnergy = uArray->GetUnitEnergy();
179 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
182 fRefArray->Add(uArray);
187 // Hadronic Correction
189 if (fApplyFractionHadronicCorrection) {
191 TArrayS* matched = new TArrayS();
192 GetTracksPointingToCell(matched, etaD, phiD,0.02);
195 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
197 const short indexS = matched->At(itrack);
200 AliAODTrack *track = fAOD->GetTrack(indexS);
205 correction = ptot * fFractionHadronicCorrection;
207 if (ptot>0 && fDebug>1 ) {
208 printf("SUM of track momentum for this cell: p=%f \n", ptot);
209 printf("fractional correction=%f \n", fFractionHadronicCorrection);
210 printf("TOTAL correction=%f \n", correction );
215 }//end hadronic correction
217 double enerCorr = digitEn;
218 if (correction >= enerCorr) enerCorr = 0;
219 else enerCorr -= correction;
220 if (correction>0 && fDebug>1) {
221 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
224 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
228 uArray->SetUnitDetectorFlag(kAll);
229 else uArray->SetUnitDetectorFlag(kEmcal);
231 uArray->SetUnitEnergy(unitEnergy+digitEt);
233 uArray->SetUnitCutFlag(kPtHigher);
235 // To be modified !!!
236 uArray->SetUnitSignalFlag(kGood);
238 // This is for jet multiplicity
239 uArray->SetUnitClusterID(index);
241 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
243 } // End loop over cells
244 } // end if !fCluster
245 else { // Keep digit information from clusterization
247 // Loop over calo clusters
248 //------------------------------------------------------------------
250 // select EMCAL clusters only
251 TRefArray * caloClusters = new TRefArray();
252 fAOD->GetEMCALClusters(caloClusters);
254 // Total number of EMCAL cluster
255 Int_t nclus = caloClusters->GetEntries() ;
260 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
262 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
263 // Retrieve cluster from aod
264 fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
266 // Get the cluster info
268 fClus->GetPosition(pos) ;
269 TVector3 vpos(pos[0],pos[1],pos[2]) ;
271 Int_t digMult = fClus->GetNCells() ;
272 UShort_t *digID = fClus->GetCellsAbsId() ;
273 Int_t trackIndex = -1;
274 if(fClus->GetNTracksMatched()!=0)
275 trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
277 // Do double-counted electron correction
278 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
280 // The electron correction go there
281 // Under construction !!!!
282 } // End of Electron correction
284 // Get CaloCells of cluster and fill the unitArray
285 for(Int_t i = 0; i < digMult ; i++) {
286 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
287 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
289 // Calibration for an energy in GeV
290 Float_t digitEn = digitAmp*0.0153;
292 Float_t etaD=-10., phiD=-10.;
293 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
294 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
296 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
298 // Hadronic Correction
300 if (fApplyFractionHadronicCorrection) {
302 TArrayS* matched = new TArrayS();
303 GetTracksPointingToCell(matched, etaD, phiD,0.02);
306 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
308 const short indexS = matched->At(itrack);
311 AliAODTrack *track = fAOD->GetTrack(indexS);
316 correction = ptot * fFractionHadronicCorrection;
318 if (ptot>0 && fDebug>1 ) {
319 printf("SUM of track momentum for this cell: p=%f \n", ptot);
320 printf("fractional correction=%f \n", fFractionHadronicCorrection);
321 printf("TOTAL correction=%f \n", correction );
326 }//end hadronic correction
328 double enerCorr = digitEn;
329 if (correction >= enerCorr) enerCorr = 0;
330 else enerCorr -= correction;
331 if (correction>0 && fDebug>1) {
332 printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
335 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
337 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
338 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
339 uArray->SetUnitTrackID(digitID);
341 Float_t unitEnergy = 0.;
343 unitEnergy = uArray->GetUnitEnergy();
347 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
350 fRefArray->Add(uArray);
357 uArray->SetUnitDetectorFlag(kAll);
358 else uArray->SetUnitDetectorFlag(kEmcal);
360 uArray->SetUnitEnergy(unitEnergy+digitEt);
362 uArray->SetUnitCutFlag(kPtHigher);
364 // Signal or background
365 // To be modified !!!
366 uArray->SetUnitSignalFlag(kGood);
368 // This is for jet multiplicity
369 uArray->SetUnitClusterID(index);
371 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
373 }// End loop over cells
374 } // End loop over clusters
381 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
382 printf("goodDigit : %d\n", goodDigit);
386 //____________________________________________________________________________
387 void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
392 for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop
394 AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
395 AliAODPid* pid = (AliAODPid*) track->GetDetPid();
399 pid->GetEMCALPosition(emcpos);
400 TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
402 Double_t deta = tpos.Eta() - eta;
403 Double_t dphi = tpos.Phi() - phi;
405 if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
406 if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
408 Double_t res = sqrt(dphi*dphi + deta*deta);
411 //add this track-index
414 array->AddAt( itrk, (size-1) );
415 if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
424 //____________________________________________________________________________
425 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
427 // Set calibration parameters:
428 // if calibration database exists, they are read from database,
429 // otherwise, they are taken from digitizer.
431 // It is a user responsilibity to open CDB before reconstruction,
433 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
435 //Check if calibration is stored in data base
437 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
439 AliCDBEntry *entry = (AliCDBEntry*)
440 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
441 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
445 printf("************* Calibration parameters not found in CDB! ****************");
446 // AliFatal("Calibration parameters not found in CDB!");
451 //____________________________________________________________________________
452 Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
455 // Convert digitized amplitude into energy.
456 // Calibration parameters are taken from calibration data base for raw data,
457 // or from digitizer parameters for simulated data.
462 printf("************* Did not get geometry from EMCALLoader ***************");
463 // AliFatal("Did not get geometry from EMCALLoader") ;
472 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
474 // fGeom->PrintGeometry();
475 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
479 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
481 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
482 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
484 return -fADCpedestalECA + amp * fADCchannelECA ;
487 else //Return energy with default parameters if calibration is not available
488 return -fADCpedestalECA + amp * fADCchannelECA ;