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 **************************************************************************/
18 // Called by ESD reader for jet analysis
19 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 // --- Standard library ---
22 #include <Riostream.h>
25 // --- ROOT system ---
27 #include <TLorentzVector.h>
28 #include <TRefArray.h>
30 #include <TGeoManager.h>
32 #include <TClonesArray.h>
35 // --- AliRoot header files ---
36 #include "AliJetFinder.h"
37 #include "AliJetReaderHeader.h"
38 #include "AliJetReader.h"
39 #include "AliJetESDReader.h"
40 #include "AliJetESDReaderHeader.h"
41 #include "AliESDEvent.h"
42 #include "AliESDVertex.h"
43 #include "AliJetDummyGeo.h"
44 #include "AliESDCaloCluster.h"
45 #include "AliESDCaloCells.h"
46 #include "AliJetUnitArray.h"
47 #include "AliJetFillUnitArrayEMCalDigits.h"
48 // Remove CDB dependence under construction
49 //#include "AliEMCALCalibData.h"
50 //#include "AliCDBManager.h"
52 //class AliCDBStorage;
53 //#include "AliCDBEntry.h"
56 ClassImp(AliJetFillUnitArrayEMCalDigits)
58 //_____________________________________________________________________________
59 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
60 : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
93 //_____________________________________________________________________________
94 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent *esd)
95 : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
128 //_____________________________________________________________________________
129 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det)
130 : TTask(det),//"AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
132 fNumUnits(det.fNumUnits),
133 fEtaMinCal(det.fEtaMinCal),
134 fEtaMaxCal(det.fEtaMaxCal),
135 fPhiMinCal(det.fPhiMinCal),
136 fPhiMaxCal(det.fPhiMaxCal),
139 fCluster(det.fCluster),
141 fNCEMCAL(det.fNCEMCAL),
142 fNCPHOS(det.fNCPHOS),
143 fNCCalo(det.fNCCalo),
144 fTPCGrid(det.fTPCGrid),
145 fEMCalGrid(det.fEMCalGrid),
146 fECorrection(det.fECorrection),
147 fReaderHeader(det.fReaderHeader),
148 fMomentumArray(det.fMomentumArray),
149 fUnitArray(det.fUnitArray),
150 fRefArray(det.fRefArray),
151 fProcId(det.fProcId),
154 fNDigitEmcal(det.fNDigitEmcal),
155 fNDigitEmcalCut(det.fNDigitEmcalCut),
156 fCalibData(det.fCalibData),
157 fADCchannelECA(det.fADCchannelECA),
158 fADCpedestalECA(det.fADCpedestalECA)
163 //_____________________________________________________________________________
164 AliJetFillUnitArrayEMCalDigits& AliJetFillUnitArrayEMCalDigits::operator=(const AliJetFillUnitArrayEMCalDigits& other)
169 fNumUnits = other.fNumUnits;
170 fEtaMinCal = other.fEtaMinCal;
171 fEtaMaxCal = other.fEtaMaxCal;
172 fPhiMinCal = other.fPhiMinCal;
173 fPhiMaxCal = other.fPhiMaxCal;
176 fCluster = other.fCluster;
177 fDebug = other.fDebug;
178 fNCEMCAL = other.fNCEMCAL;
179 fNCPHOS = other.fNCPHOS;
180 fNCCalo = other.fNCCalo;
181 fTPCGrid = other.fTPCGrid;
182 fEMCalGrid = other.fEMCalGrid;
183 fECorrection = other.fECorrection;
184 fReaderHeader = other.fReaderHeader;
185 fMomentumArray = other.fMomentumArray;
186 fUnitArray = other.fUnitArray;
187 fRefArray = other.fRefArray;
188 fProcId = other.fProcId;
191 fNDigitEmcal = other.fNDigitEmcal;
192 fNDigitEmcalCut = other.fNDigitEmcalCut;
193 fCalibData = other.fCalibData;
194 fADCchannelECA = other.fADCchannelECA;
195 fADCpedestalECA = other.fADCpedestalECA;
201 //____________________________________________________________________________
202 void AliJetFillUnitArrayEMCalDigits::InitParameters()
204 fNumUnits = fGeom->GetNCells(); // Number of towers in EMCAL
206 fEtaMinCal = fGeom->GetArm1EtaMin();
207 fEtaMaxCal = fGeom->GetArm1EtaMax();
208 fPhiMinCal = fGeom->GetArm1PhiMin();
209 fPhiMaxCal = fGeom->GetArm1PhiMax();
212 // Get calibration parameters from file or digitizer default values.
213 // Under construction
214 // GetCalibrationParameters() ;
216 if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
220 //_____________________________________________________________________________
221 AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
226 //_____________________________________________________________________________
227 void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
233 fDebug = fReaderHeader->GetDebug();
234 fOpt = fReaderHeader->GetDetector();
235 fCluster = fReaderHeader->GetCluster();
243 if(!fCluster) { // Keep all digit information
244 // Loop over all cell information
245 //------------------------------------------------------------------
246 AliESDCaloCells &cells = *(fESD->GetEMCALCells());
247 Int_t ncell = cells.GetNumberOfCells() ;
249 for (Int_t icell= 0; icell < ncell; icell++) {
250 Int_t digitID = cells.GetCellNumber(icell);
251 Double_t digitAmp = cells.GetAmplitude(icell);
252 Float_t digitEn = digitAmp*0.0153; // Last correct
253 // Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
255 Float_t etaD=-10., phiD=-10.;
256 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
257 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
259 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
261 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
263 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
265 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
266 uArray->SetUnitTrackID(digitID);
268 Float_t unitEnergy = 0.;
270 unitEnergy = uArray->GetUnitEnergy();
274 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
277 fRefArray->Add(uArray);
284 uArray->SetUnitDetectorFlag(kAll);
285 else uArray->SetUnitDetectorFlag(kEmcal);
287 uArray->SetUnitEnergy(unitEnergy+digitEt);
289 uArray->SetUnitCutFlag(kPtHigher);
291 // To be modified !!!
292 uArray->SetUnitSignalFlag(kGood);
294 // This is for jet multiplicity
295 uArray->SetUnitClusterID(index);
297 if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
299 } // End loop over cells
301 } // end if !fCluster
302 else { // Keep digit information from clusterization
304 // Loop over calo clusters
305 //------------------------------------------------------------------
307 //select EMCAL clusters only
308 TRefArray * caloClusters = new TRefArray();
309 fESD->GetEMCALClusters(caloClusters);
311 // Total number of EMCAL cluster
312 Int_t nclus = caloClusters->GetEntries() ;
316 // Get reconstructed vertex position
317 Double_t vertex_position[3] ;
318 fESD->GetVertex()->GetXYZ(vertex_position) ;
321 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
323 for(Int_t j = beg; j < nclus; j++) { // loop over clusters
324 // Retrieve cluster from esd
325 fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
327 // Get the cluster info
328 Float_t energy = fClus->E() ;
329 Int_t iprim = fClus->GetLabel();
330 Int_t trackIndex = fClus->GetTrackMatched();
332 fClus->GetPosition(pos) ;
333 TVector3 vpos(pos[0],pos[1],pos[2]) ;
335 fClus->GetMomentum(p,vertex_position);
337 Int_t digMult = fClus->GetNCells() ;
338 UShort_t *digID = fClus->GetCellsAbsId() ;
340 if(fDebug>2) cout<<"Cluster "<< j <<"; digits mult "<<digMult<<"; type "<<(Int_t )fClus->GetClusterType()
341 <<"; Energy "<<energy<< "; transverse energy:" << energy*TMath::Abs(TMath::Sin(EtaToTheta(vpos.Eta())))
342 <<"; Phi "<<vpos.Phi()<<"; Eta "<<vpos.Eta() <<"; label "<<iprim<<endl;
344 // Do double-counted electron correction
345 if (fECorrection != 0 && trackIndex !=-1 )
347 // The electron correction go there
348 // Under construction !!!!
350 } // End of Electron correction
352 // Get CaloCells of cluster and fill the unitArray
353 for(Int_t i = 0; i < digMult ; i++)
355 Int_t digitID = digID[i]; // or clus->GetCellNumber(i) ;
356 Double_t digitAmp = cells.GetCellAmplitude(digitID) ;
358 // Calibration for an energy in GeV
359 Float_t digitEn = digitAmp*0.0153;
360 // Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
362 Float_t etaD=-10., phiD=-10.;
363 fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
364 // fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
366 phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
368 Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
370 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
371 if(uArray->GetUnitEnergy() == 0.) goodDigit++;
372 uArray->SetUnitTrackID(digitID);
374 Float_t unitEnergy = 0.;
376 unitEnergy = uArray->GetUnitEnergy();
380 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
383 fRefArray->Add(uArray);
390 uArray->SetUnitDetectorFlag(kAll);
391 else uArray->SetUnitDetectorFlag(kEmcal);
393 uArray->SetUnitEnergy(unitEnergy+digitEt);
394 uArray->SetUnitCutFlag(kPtHigher);
396 // To be modified !!!
397 uArray->SetUnitSignalFlag(kGood);
399 // This is for jet multiplicity
400 uArray->SetUnitClusterID(index);
402 } // End loop over cells
403 } // End loop over clusters
410 printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
411 printf("goodDigit : %d\n", goodDigit);
416 // //____________________________________________________________________________
417 // void AliJetFillUnitArrayEMCalDigits::GetCalibrationParameters()
419 // // Set calibration parameters:
420 // // if calibration database exists, they are read from database,
421 // // otherwise, they are taken from digitizer.
423 // // It is a user responsilibity to open CDB before reconstruction,
425 // // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
427 // //Check if calibration is stored in data base
429 // if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
431 // AliCDBEntry *entry = (AliCDBEntry*)
432 // AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
433 // if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
437 // printf("************* Calibration parameters not found in CDB! ****************");
438 // // AliFatal("Calibration parameters not found in CDB!");
443 // //____________________________________________________________________________
444 // Float_t AliJetFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
447 // // Convert digitized amplitude into energy.
448 // // Calibration parameters are taken from calibration data base for raw data,
449 // // or from digitizer parameters for simulated data.
454 // printf("************* Did not get geometry from EMCALLoader ***************");
456 // Int_t iSupMod = -1;
457 // Int_t nModule = -1;
463 // Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
465 // // fGeom->PrintGeometry();
466 // Error("Calibrate()"," Wrong cell id number : %i", AbsId);
470 // fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
472 // fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
473 // fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
475 // return -fADCpedestalECA + amp * fADCchannelECA ;
478 // else //Return energy with default parameters if calibration is not available
479 // return -fADCpedestalECA + amp * fADCchannelECA ;
484 //_____________________________________________________________________________
485 Float_t AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
487 // return (180./TMath::Pi())*2.*atan(exp(-arg));
488 return 2.*atan(exp(-arg));