Adding includes for EMCAL_Utils and OADB PATH (A. Shabetai)
[u/mrichter/AliRoot.git] / JETAN / AliJetAODFillUnitArrayEMCalDigits.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 AOD reader for jet analysis
20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 //======================================================================
22
23 class TSystem;
24 class TLorentzVector;
25 class TGeoManager;
26 class TArrayS;
27
28 // --- AliRoot header files ---
29 #include "AliJetUnitArray.h"
30 #include "AliJetAODFillUnitArrayEMCalDigits.h"
31 class AliJetFinder;
32
33 // #include "AliEMCALCalibData.h"
34 // #include "AliCDBManager.h"
35 // class AliCDBStorage;
36 // #include "AliCDBEntry.h"
37
38
39 ClassImp(AliJetAODFillUnitArrayEMCalDigits)
40
41 //_____________________________________________________________________________
42 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits()
43   : AliJetFillUnitArray(),
44     fAOD(0),
45     fNIn(0),
46     fCluster(0),
47     fNCEMCAL(0),
48     fNCPHOS(0),
49     fNCCalo(0),
50     fApplyElectronCorrection(kFALSE),
51     fApplyFractionHadronicCorrection(kFALSE),
52     fFractionHadronicCorrection(0.3),
53     fClus(0x0),
54     fNDigitEmcal(0),
55     fNDigitEmcalCut(0)
56 {
57   // constructor
58 }
59
60 //_____________________________________________________________________________
61 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(AliAODEvent *aod)
62   : AliJetFillUnitArray(),
63     fAOD(aod),
64     fNIn(0),
65     fCluster(0),
66     fNCEMCAL(0),
67     fNCPHOS(0),
68     fNCCalo(0),
69     fApplyElectronCorrection(kFALSE),
70     fApplyFractionHadronicCorrection(kFALSE),
71     fFractionHadronicCorrection(0.3),
72     fClus(0x0),
73     fNDigitEmcal(0),
74     fNDigitEmcalCut(0)
75 {
76   // constructor
77 }
78
79 //_____________________________________________________________________________
80 AliJetAODFillUnitArrayEMCalDigits::AliJetAODFillUnitArrayEMCalDigits(const AliJetAODFillUnitArrayEMCalDigits &det)
81   : AliJetFillUnitArray(det),
82     fAOD(det.fAOD),
83     fNIn(det.fNIn),
84     fCluster(det.fCluster),
85     fNCEMCAL(det.fNCEMCAL),
86     fNCPHOS(det.fNCPHOS),
87     fNCCalo(det.fNCCalo),
88     fApplyElectronCorrection(det.fApplyElectronCorrection),
89     fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
90     fFractionHadronicCorrection(det.fFractionHadronicCorrection),
91     fClus(det.fClus),
92     fNDigitEmcal(det.fNDigitEmcal),
93     fNDigitEmcalCut(det.fNDigitEmcalCut)
94 {
95   // Copy constructor
96 }
97
98 //_____________________________________________________________________________
99 AliJetAODFillUnitArrayEMCalDigits& AliJetAODFillUnitArrayEMCalDigits::operator=(const AliJetAODFillUnitArrayEMCalDigits& other)
100 {
101   // Assignment
102
103   fAOD = other.fAOD;
104   fNIn = other.fNIn;
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;
112   fClus = other.fClus;
113   fNDigitEmcal = other.fNDigitEmcal;
114   fNDigitEmcalCut = other.fNDigitEmcalCut;
115
116   return (*this);
117
118 }
119
120 //_____________________________________________________________________________
121 AliJetAODFillUnitArrayEMCalDigits::~AliJetAODFillUnitArrayEMCalDigits()
122 {
123   // destructor
124 }
125
126
127
128 //_____________________________________________________________________________
129 void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
130 //void AliJetAODFillUnitArrayEMCalDigits::Exec(Option_t* option)
131 {
132   //
133   // Main method.
134   // Fill the unit array with the neutral particle information from the EMCal cells in AOD
135
136   fDebug = fReaderHeader->GetDebug();
137   fOpt = fReaderHeader->GetDetector();
138   fCluster = fReaderHeader->GetCluster();
139
140   // Init parameters
141   //  InitParameters();
142
143   Int_t   goodDigit      = 0;
144   Int_t   index          = 0;
145
146   if(!fCluster) { // Keep all digit information
147
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();
153
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);
159
160       Float_t etaD=-10., phiD=-10.;
161       fGeom->EtaPhiFromIndex(digitID,etaD,phiD); 
162       //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
163
164       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
165       
166       //      Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
167
168       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
169
170       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
171       uArray->SetUnitTrackID(digitID);
172       
173       Float_t unitEnergy = 0.;
174       Bool_t ok = kFALSE;
175       unitEnergy = uArray->GetUnitEnergy();
176
177       if(unitEnergy==0){
178         if(!fProcId){
179           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
180           fProcId = kTRUE;
181         }
182         fRefArray->Add(uArray);
183         fNDigitEmcal++;
184         ok = kTRUE;
185       }
186
187       // Hadronic Correction
188       double correction=0;
189       if (fApplyFractionHadronicCorrection) {
190       
191         TArrayS* matched = new TArrayS();
192         GetTracksPointingToCell(matched, etaD, phiD,0.02);
193
194         double ptot = 0;
195         for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
196           {
197             const short indexS = matched->At(itrack);
198             if (indexS>-1)
199               { 
200                 AliAODTrack *track = fAOD->GetTrack(indexS);
201                 ptot += track->P();
202               }
203           }
204         
205         correction = ptot * fFractionHadronicCorrection;
206
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 );
211         }
212         
213         delete matched;
214         
215       }//end hadronic correction
216
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);
222       }
223
224       Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
225
226       // Detector flag
227       if(unitEnergy>0.)
228         uArray->SetUnitDetectorFlag(kAll);
229       else uArray->SetUnitDetectorFlag(kEmcal);
230       
231       uArray->SetUnitEnergy(unitEnergy+digitEt);
232
233       uArray->SetUnitCutFlag(kPtHigher);
234
235       // To be modified !!!
236       uArray->SetUnitSignalFlag(kGood);
237         
238       // This is for jet multiplicity
239       uArray->SetUnitClusterID(index);
240         
241       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
242       
243     } // End loop over cells
244   } // end if !fCluster
245   else { // Keep digit information from clusterization
246
247     // Loop over calo clusters
248     //------------------------------------------------------------------
249
250     // select EMCAL clusters only
251     TRefArray * caloClusters  = new TRefArray();
252     fAOD->GetEMCALClusters(caloClusters);
253
254     // Total number of EMCAL cluster
255     Int_t nclus = caloClusters->GetEntries() ;
256     Int_t beg   = 0;
257     Float_t pos[3] ;
258
259     // Get CaloCells
260     AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
261
262     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
263       // Retrieve cluster from aod
264       fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
265
266       // Get the cluster info
267
268       fClus->GetPosition(pos) ;
269       TVector3 vpos(pos[0],pos[1],pos[2]) ;
270
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();
276
277       // Do double-counted electron correction 
278       if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
279         {
280           // The electron correction go there
281           // Under construction !!!!
282         }  // End of Electron correction
283
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) ;
288
289         // Calibration for an energy in GeV
290         Float_t digitEn = digitAmp*0.0153;
291
292         Float_t etaD=-10., phiD=-10.;
293         fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
294         //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
295
296         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
297
298         // Hadronic Correction
299         double correction=0;
300         if (fApplyFractionHadronicCorrection) {
301       
302           TArrayS* matched = new TArrayS();
303           GetTracksPointingToCell(matched, etaD, phiD,0.02);
304           
305           double ptot = 0;
306           for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
307             {
308               const short indexS = matched->At(itrack);
309               if (indexS>-1)
310                 {       
311                   AliAODTrack *track = fAOD->GetTrack(indexS);
312                   ptot += track->P();
313                 }
314             }
315         
316           correction = ptot * fFractionHadronicCorrection;
317
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 );
322           }
323         
324           delete matched;
325         
326         }//end hadronic correction
327         
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);
333         }
334
335         Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
336
337         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
338         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
339         uArray->SetUnitTrackID(digitID);
340
341         Float_t unitEnergy = 0.;
342         Bool_t ok = kFALSE;
343         unitEnergy = uArray->GetUnitEnergy();
344
345         if(unitEnergy==0){
346           if(!fProcId){
347             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
348             fProcId = kTRUE;
349           }
350           fRefArray->Add(uArray);
351           fNDigitEmcal++;
352           ok = kTRUE;
353         }
354
355         // Detector flag
356         if(unitEnergy>0.)
357           uArray->SetUnitDetectorFlag(kAll);
358         else uArray->SetUnitDetectorFlag(kEmcal);
359       
360         uArray->SetUnitEnergy(unitEnergy+digitEt);
361
362         uArray->SetUnitCutFlag(kPtHigher);
363
364         // Signal or background
365         // To be modified !!!
366         uArray->SetUnitSignalFlag(kGood);
367
368         // This is for jet multiplicity
369         uArray->SetUnitClusterID(index);
370         
371         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
372
373       }// End loop over cells
374     } // End loop over clusters
375   } // end else
376
377   fNIn += goodDigit;
378
379   if(fDebug>1)
380     {
381       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
382       printf("goodDigit : %d\n", goodDigit);
383     }
384 }
385
386 //____________________________________________________________________________
387 void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
388 {
389 // Get all tracks pointing to cell 
390   int size=0;
391   
392   for (Int_t itrk =  0; itrk <  fAOD->GetNumberOfTracks() ; itrk++) { //track loop
393   
394     AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;  
395     AliAODPid*    pid   = (AliAODPid*) track->GetDetPid();
396     
397     if(pid) {
398       Double_t emcpos[3];
399       pid->GetEMCALPosition(emcpos);      
400       TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
401
402       Double_t deta = tpos.Eta() - eta;
403       Double_t dphi = tpos.Phi() - phi;
404
405       if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
406       if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
407
408       Double_t res = sqrt(dphi*dphi + deta*deta);
409       
410       if (res< cut) {
411         //add this track-index
412         size++;
413         array->Set( size );
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);
416                       
417       }
418     }
419   }
420
421 }
422
423 /*
424 //____________________________________________________________________________
425 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
426 {
427   // Set calibration parameters:
428   // if calibration database exists, they are read from database,
429   // otherwise, they are taken from digitizer.
430   //
431   // It is a user responsilibity to open CDB before reconstruction,
432   // for example:
433   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
434
435   //Check if calibration is stored in data base
436
437   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
438     {
439       AliCDBEntry *entry = (AliCDBEntry*)
440         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
441       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
442     }
443
444   if(!fCalibData)
445     printf("************* Calibration parameters not found in CDB! ****************");
446 //    AliFatal("Calibration parameters not found in CDB!");
447
448
449 }
450
451 //____________________________________________________________________________
452 Float_t  AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
453 {
454
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.
458
459   if(fCalibData){
460
461     if (fGeom==0)
462       printf("************* Did not get geometry from EMCALLoader ***************");
463 //      AliFatal("Did not get geometry from EMCALLoader") ;
464
465     Int_t iSupMod = -1;
466     Int_t nModule  = -1;
467     Int_t nIphi   = -1;
468     Int_t nIeta   = -1;
469     Int_t iphi    = -1;
470     Int_t ieta    = -1;
471
472     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
473     if(!bCell) {
474       // fGeom->PrintGeometry();
475       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
476       assert(0);
477     }
478
479     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
480
481     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
482     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
483
484    return -fADCpedestalECA + amp * fADCchannelECA ;
485
486   }
487   else //Return energy with default parameters if calibration is not available
488     return -fADCpedestalECA + amp * fADCchannelECA ;
489
490 }
491 */
492
493