199d64d20d605a6d2595b7b0c08c9cf6ec35ef18
[u/mrichter/AliRoot.git] / JETAN / AliJetESDFillUnitArrayEMCalDigits.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
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  **************************************************************************/
16  
17 //---------------------------------------------------------------------
18 // Fill Unit Array 
19 // Called by ESD reader for jet analysis
20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 //---------------------------------------------------------------------
22
23 // --- ROOT system ---
24 #include <TLorentzVector.h>
25 #include <TMath.h>
26
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"
36
37 // --- ROOT system ---
38 class TSystem;
39 class TGeoManager;
40 class TArrayS;
41
42 // --- AliRoot header files ---
43 class AliJetFinder;
44 class AliJetReader;
45 class AliJetESDReader;
46 // class AliEMCALCalibData;
47 // class AliCDBManager;
48 // class AliCDBStorage;
49 // class AliCDBEntry;
50
51 ClassImp(AliJetESDFillUnitArrayEMCalDigits)
52
53 //_____________________________________________________________________________
54 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
55   : AliJetFillUnitArray(),
56     fESD(0),
57     fNIn(0),
58     fOpt(0),
59     fCluster(0),
60     fDebug(0),
61     fNCEMCAL(0),
62     fNCPHOS(0),
63     fNCCalo(0),
64     fApplyElectronCorrection(kFALSE),
65     fApplyFractionHadronicCorrection(kFALSE),
66     fFractionHadronicCorrection(0.3),
67     fClus(0x0),
68     fNDigitEmcal(0),
69     fNDigitEmcalCut(0)
70 {
71   // constructor
72 }
73
74 //_____________________________________________________________________________
75 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
76   : AliJetFillUnitArray(),
77     fESD(esd),
78     fNIn(0),
79     fOpt(0),
80     fCluster(0),
81     fDebug(0),
82     fNCEMCAL(0),
83     fNCPHOS(0),
84     fNCCalo(0),
85     fApplyElectronCorrection(kFALSE),
86     fApplyFractionHadronicCorrection(kFALSE),
87     fFractionHadronicCorrection(0.3),
88     fClus(0x0),
89     fNDigitEmcal(0),
90     fNDigitEmcalCut(0)
91 {
92   // constructor
93 }
94
95 //_____________________________________________________________________________
96 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
97   : AliJetFillUnitArray(det),
98     fESD(det.fESD),
99     fNIn(det.fNIn),
100     fOpt(det.fOpt),
101     fCluster(det.fCluster),
102     fDebug(det.fDebug),
103     fNCEMCAL(det.fNCEMCAL),
104     fNCPHOS(det.fNCPHOS),
105     fNCCalo(det.fNCCalo),
106     fApplyElectronCorrection(det.fApplyElectronCorrection),
107     fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
108     fFractionHadronicCorrection(det.fFractionHadronicCorrection),
109     fClus(det.fClus),
110     fNDigitEmcal(det.fNDigitEmcal),
111     fNDigitEmcalCut(det.fNDigitEmcalCut)
112 {
113   // Copy constructor
114 }
115
116 //_____________________________________________________________________________
117 AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
118 {
119   // Assignment
120
121   fESD = other.fESD;
122   fNIn = other.fNIn;
123   fOpt = other.fOpt;
124   fCluster = other.fCluster;
125   fDebug = other.fDebug;
126   fNCEMCAL = other.fNCEMCAL;
127   fNCPHOS = other.fNCPHOS;
128   fNCCalo = other.fNCCalo;
129   fApplyElectronCorrection = other.fApplyElectronCorrection;
130   fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
131   fFractionHadronicCorrection = other.fFractionHadronicCorrection;
132   fClus = other.fClus;
133   fNDigitEmcal = other.fNDigitEmcal;
134   fNDigitEmcalCut = other.fNDigitEmcalCut;
135
136   return (*this);
137
138 }
139
140 //_____________________________________________________________________________
141 AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
142 {
143   // destructor
144 }
145
146 //_____________________________________________________________________________
147 void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
148 {
149   //
150   // Main method.
151   // Fill the unit array with the neutral particle information from the EMCal cells in ESD
152   //
153
154   fDebug = fReaderHeader->GetDebug();
155   fOpt = fReaderHeader->GetDetector();
156   fCluster = fReaderHeader->GetCluster();
157
158   // Init parameters
159   InitParameters();
160
161 //(not used ?)  Int_t   nDigitTot      = 0;
162   Int_t   goodDigit      = 0;
163 //(not used ?)  Int_t   beg            = 0;
164 //(not used ?)  Int_t   end            = 0; 
165   Int_t   index          = 0;
166 //(not used ?)  Int_t   count          = 0;
167
168 //(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
169
170   if(!fCluster) { // Keep all digit information
171     // Loop over all cell information
172     //------------------------------------------------------------------
173     AliESDCaloCells &cells= *(fESD->GetEMCALCells());
174     Int_t ncell = cells.GetNumberOfCells() ;
175 //(not used ?)    Int_t type = cells.GetType();
176
177     for (Int_t icell=  0; icell <  ncell; icell++) {
178       Int_t   digitID   = cells.GetCellNumber(icell);
179       Float_t digitAmp  = cells.GetAmplitude(icell);
180 //(not used ?)      Float_t digitTime = cells.GetTime(icell);
181       Float_t digitEn   = digitAmp*0.0153; // Last correct
182       //  Float_t digitEn = Calibrate(digitAmp,digitID);
183
184       Float_t etaD=-10., phiD=-10.;
185       fGeom->EtaPhiFromIndex(digitID,etaD,phiD); 
186       //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
187
188       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
189       
190       Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
191
192       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
193
194       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
195       uArray->SetUnitTrackID(digitID);
196
197       Float_t unitEnergy = 0.;
198       Bool_t ok = kFALSE;
199       unitEnergy = uArray->GetUnitEnergy();
200
201       if(unitEnergy==0){
202         if(!fProcId){
203           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
204           fProcId = kTRUE;
205         }
206         fRefArray->Add(uArray);
207         fNDigitEmcal++;
208         ok = kTRUE;
209       }
210
211       // Detector flag
212       if(unitEnergy>0.)
213         uArray->SetUnitDetectorFlag(kAll);
214       else uArray->SetUnitDetectorFlag(kEmcal);
215       
216       uArray->SetUnitEnergy(unitEnergy+digitEt);
217
218       uArray->SetUnitCutFlag(kPtHigher);
219
220       // To be modified !!!
221       uArray->SetUnitSignalFlag(kGood);
222
223       // This is for jet multiplicity
224       uArray->SetUnitClusterID(index);
225         
226       if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
227
228     } // End loop over cells
229   } // end if !fCluster
230   else { // Keep digit information from clusterization
231
232     // Loop over calo clusters
233     //------------------------------------------------------------------
234
235     //select EMCAL clusters only
236     TRefArray * caloClusters  = new TRefArray();
237     fESD->GetEMCALClusters(caloClusters);
238
239     // Total number of EMCAL cluster
240     Int_t nclus = caloClusters->GetEntries() ;
241     Int_t beg   = 0;
242     Float_t pos[3] ;
243
244     // Get reconstructed vertex position
245     Double_t vertexPosition[3] ;
246     fESD->GetVertex()->GetXYZ(vertexPosition) ;
247
248     // Get CaloCells
249     AliESDCaloCells &cells= *(fESD->GetEMCALCells());
250 //(not used ?)    Int_t ncell = cells.GetNumberOfCells() ;
251 //(not used ?)    Int_t type = cells.GetType();
252
253     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
254       // Retrieve cluster from esd
255       AliESDCaloCluster *fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
256
257       // Get the cluster info
258 //(not used ?)      Float_t energy         = fClus->E() ;
259 //(not used ?)      Int_t   iprim          = fClus->GetLabel();
260
261       fClus->GetPosition(pos) ;
262       TVector3 vpos(pos[0],pos[1],pos[2]) ;
263       TLorentzVector p ;
264       fClus->GetMomentum(p,vertexPosition);
265
266       Int_t     digMult = fClus->GetNCells() ;
267       UShort_t *digID   = fClus->GetCellsAbsId() ;
268 //(not used ?)      Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
269       Int_t     trackIndex = fClus->GetTrackMatched();
270
271       // Do double-counted electron correction 
272       if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
273         {
274           // The electron correction go there
275           // Under construction !!!!
276         }  // End of Electron correction
277
278       // Get CaloCells of cluster and fill the unitArray
279       for(Int_t i = 0; i < digMult ; i++){
280         Int_t    digitID   = digID[i]; // or clus->GetCellNumber(i) ;
281 //(not used ?)        Double_t digitAmpFrac = digAmpFrac[i];
282         Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
283 //(not used ?)        Float_t  digitTime    = cells.GetCellTime(digitID);
284
285         // Calibration for an energy in GeV
286         Float_t digitEn = digitAmp*0.0153;
287
288         Float_t etaD=-10., phiD=-10.;
289         fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
290         //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
291
292         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
293
294         Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
295
296         cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
297
298         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
299         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
300         uArray->SetUnitTrackID(digitID);
301
302         Float_t unitEnergy = 0.;
303         Bool_t ok = kFALSE;
304         unitEnergy = uArray->GetUnitEnergy();
305
306         if(unitEnergy==0){
307           if(!fProcId){
308             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
309             fProcId = kTRUE;
310           }
311           fRefArray->Add(uArray);
312           fNDigitEmcal++;
313           ok = kTRUE;
314         }
315
316         // Detector flag
317         if(unitEnergy>0.)
318           uArray->SetUnitDetectorFlag(kAll);
319         else uArray->SetUnitDetectorFlag(kEmcal);
320       
321         uArray->SetUnitEnergy(unitEnergy+digitEt);
322
323         uArray->SetUnitCutFlag(kPtHigher);
324
325         // To be modified !!!
326         uArray->SetUnitSignalFlag(kGood);
327
328         // This is for jet multiplicity
329         uArray->SetUnitClusterID(index);
330         
331         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
332       
333       } // End loop over cells
334     } // End loop over clusters
335   } // end else
336
337   fNIn += goodDigit;
338
339   if(fDebug>1)
340     {
341       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
342       printf("goodDigit : %d\n", goodDigit);
343     }
344 }
345
346 /*
347 //____________________________________________________________________________
348 void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
349 {
350   // Set calibration parameters:
351   // if calibration database exists, they are read from database,
352   // otherwise, they are taken from digitizer.
353   //
354   // It is a user responsilibity to open CDB before reconstruction,
355   // for example:
356   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
357
358   //Check if calibration is stored in data base
359
360   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
361     {
362       AliCDBEntry *entry = (AliCDBEntry*)
363         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
364       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
365     }
366
367   if(!fCalibData)
368     printf("************* Calibration parameters not found in CDB! ****************");
369 //    AliFatal("Calibration parameters not found in CDB!");
370
371
372 }
373
374 //____________________________________________________________________________
375 Float_t  AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
376 {
377
378   // Convert digitized amplitude into energy.
379   // Calibration parameters are taken from calibration data base for raw data,
380   // or from digitizer parameters for simulated data.
381
382   if(fCalibData){
383
384     if (fGeom==0)
385       printf("************* Did not get geometry from EMCALLoader ***************");
386 //      AliFatal("Did not get geometry from EMCALLoader") ;
387
388     Int_t iSupMod = -1;
389     Int_t nModule  = -1;
390     Int_t nIphi   = -1;
391     Int_t nIeta   = -1;
392     Int_t iphi    = -1;
393     Int_t ieta    = -1;
394
395     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
396     if(!bCell) {
397       // fGeom->PrintGeometry();
398       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
399       assert(0);
400     }
401
402     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
403
404     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
405     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
406
407    return -fADCpedestalECA + amp * fADCchannelECA ;
408
409   }
410   else //Return energy with default parameters if calibration is not available
411     return -fADCpedestalECA + amp * fADCchannelECA ;
412
413 }
414 */
415