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