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 // Temporarily commented - will be updated removing the AliAODpid dependence
191 if (fApplyFractionHadronicCorrection) {
192 TArrayS* matched = new TArrayS();
193 GetTracksPointingToCell(matched, etaD, phiD,0.02);
196 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
198 const short indexS = matched->At(itrack);
201 AliAODTrack *track = fAOD->GetTrack(indexS);
206 correction = ptot * fFractionHadronicCorrection;
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 );
216 }//end hadronic correction
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);
225 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
229 uArray->SetUnitDetectorFlag(kAll);
230 else uArray->SetUnitDetectorFlag(kEmcal);
232 uArray->SetUnitEnergy(unitEnergy+digitEt);
234 uArray->SetUnitCutFlag(kPtHigher);
236 // To be modified !!!
237 uArray->SetUnitSignalFlag(kGood);
239 // This is for jet multiplicity
240 uArray->SetUnitClusterID(index);
242 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
244 } // End loop over cells
245 } // end if !fCluster
246 else { // Keep digit information from clusterization
248 // Loop over calo clusters
249 //------------------------------------------------------------------
251 // select EMCAL clusters only
252 TRefArray * caloClusters = new TRefArray();
253 fAOD->GetEMCALClusters(caloClusters);
255 // Total number of EMCAL cluster
256 Int_t nclus = caloClusters->GetEntries() ;
261 AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
263 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
264 // Retrieve cluster from aod
265 fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
267 // Get the cluster info
269 fClus->GetPosition(pos) ;
270 TVector3 vpos(pos[0],pos[1],pos[2]) ;
272 Int_t digMult = fClus->GetNCells() ;
273 UShort_t *digID = fClus->GetCellsAbsId() ;
274 Int_t trackIndex = -1;
275 if(fClus->GetNTracksMatched()!=0)
276 trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
278 // Do double-counted electron correction
279 if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
281 // The electron correction go there
282 // Under construction !!!!
283 } // End of Electron correction
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) ;
288 Float_t digitAmp = cells.GetCellAmplitude(digitID) ;
290 // Calibration for an energy in GeV
291 Float_t digitEn = digitAmp*0.0153;
293 Float_t etaD=-10., phiD=-10.;
294 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
295 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
297 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
299 // Hadronic Correction
301 // Temporarily commented - will be updated removing the AliAODpid dependence
303 if (fApplyFractionHadronicCorrection) {
304 TArrayS* matched = new TArrayS();
305 GetTracksPointingToCell(matched, etaD, phiD,0.02);
308 for(int itrack = 0; itrack < matched->GetSize(); itrack++)
310 const short indexS = matched->At(itrack);
313 AliAODTrack *track = fAOD->GetTrack(indexS);
318 correction = ptot * fFractionHadronicCorrection;
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 );
328 }//end hadronic correction
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);
337 Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
339 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
340 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
341 uArray->SetUnitTrackID(digitID);
343 Float_t unitEnergy = 0.;
345 unitEnergy = uArray->GetUnitEnergy();
349 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
352 fRefArray->Add(uArray);
359 uArray->SetUnitDetectorFlag(kAll);
360 else uArray->SetUnitDetectorFlag(kEmcal);
362 uArray->SetUnitEnergy(unitEnergy+digitEt);
364 uArray->SetUnitCutFlag(kPtHigher);
366 // Signal or background
367 // To be modified !!!
368 uArray->SetUnitSignalFlag(kGood);
370 // This is for jet multiplicity
371 uArray->SetUnitClusterID(index);
373 if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
375 }// End loop over cells
376 } // End loop over clusters
383 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
384 printf("goodDigit : %d\n", goodDigit);
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
394 // for (Int_t itrk = 0; itrk < fAOD->GetNumberOfTracks() ; itrk++) { //track loop
396 // AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
397 // AliAODPid* pid = (AliAODPid*) track->GetDetPid();
400 // Double_t emcpos[3];
401 // pid->GetEMCALPosition(emcpos);
402 // TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
404 // Double_t deta = tpos.Eta() - eta;
405 // Double_t dphi = tpos.Phi() - phi;
407 // if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
408 // if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
410 // Double_t res = sqrt(dphi*dphi + deta*deta);
413 // //add this track-index
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);
426 //____________________________________________________________________________
427 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
429 // Set calibration parameters:
430 // if calibration database exists, they are read from database,
431 // otherwise, they are taken from digitizer.
433 // It is a user responsilibity to open CDB before reconstruction,
435 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
437 //Check if calibration is stored in data base
439 if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
441 AliCDBEntry *entry = (AliCDBEntry*)
442 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
443 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
447 printf("************* Calibration parameters not found in CDB! ****************");
448 // AliFatal("Calibration parameters not found in CDB!");
453 //____________________________________________________________________________
454 Float_t AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
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.
464 printf("************* Did not get geometry from EMCALLoader ***************");
465 // AliFatal("Did not get geometry from EMCALLoader") ;
474 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
476 // fGeom->PrintGeometry();
477 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
481 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
483 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
484 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
486 return -fADCpedestalECA + amp * fADCchannelECA ;
489 else //Return energy with default parameters if calibration is not available
490 return -fADCpedestalECA + amp * fADCchannelECA ;