Adding rho's depenence on trigger track (M. Verweij)
[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 // Temporarily commented - will be updated removing the AliAODpid dependence
190 /*
191       if (fApplyFractionHadronicCorrection) {
192         TArrayS* matched = new TArrayS();
193         GetTracksPointingToCell(matched, etaD, phiD,0.02);
194
195         double ptot = 0;
196         for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
197           {
198             const short indexS = matched->At(itrack);
199             if (indexS>-1)
200               { 
201                 AliAODTrack *track = fAOD->GetTrack(indexS);
202                 ptot += track->P();
203               }
204           }
205         
206         correction = ptot * fFractionHadronicCorrection;
207
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 );
212         }
213         
214         delete matched;
215         
216       }//end hadronic correction
217 */
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);
223       }
224
225       Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
226
227       // Detector flag
228       if(unitEnergy>0.)
229         uArray->SetUnitDetectorFlag(kAll);
230       else uArray->SetUnitDetectorFlag(kEmcal);
231       
232       uArray->SetUnitEnergy(unitEnergy+digitEt);
233
234       uArray->SetUnitCutFlag(kPtHigher);
235
236       // To be modified !!!
237       uArray->SetUnitSignalFlag(kGood);
238         
239       // This is for jet multiplicity
240       uArray->SetUnitClusterID(index);
241         
242       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
243       
244     } // End loop over cells
245   } // end if !fCluster
246   else { // Keep digit information from clusterization
247
248     // Loop over calo clusters
249     //------------------------------------------------------------------
250
251     // select EMCAL clusters only
252     TRefArray * caloClusters  = new TRefArray();
253     fAOD->GetEMCALClusters(caloClusters);
254
255     // Total number of EMCAL cluster
256     Int_t nclus = caloClusters->GetEntries() ;
257     Int_t beg   = 0;
258     Float_t pos[3] ;
259
260     // Get CaloCells
261     AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
262
263     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
264       // Retrieve cluster from aod
265       fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
266
267       // Get the cluster info
268
269       fClus->GetPosition(pos) ;
270       TVector3 vpos(pos[0],pos[1],pos[2]) ;
271
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();
277
278       // Do double-counted electron correction 
279       if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
280         {
281           // The electron correction go there
282           // Under construction !!!!
283         }  // End of Electron correction
284
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) ;
289
290         // Calibration for an energy in GeV
291         Float_t digitEn = digitAmp*0.0153;
292
293         Float_t etaD=-10., phiD=-10.;
294         fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
295         //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
296
297         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
298
299         // Hadronic Correction
300         double correction=0;
301 // Temporarily commented - will be updated removing the AliAODpid dependence
302 /*
303         if (fApplyFractionHadronicCorrection) {
304           TArrayS* matched = new TArrayS();
305           GetTracksPointingToCell(matched, etaD, phiD,0.02);
306           
307           double ptot = 0;
308           for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
309             {
310               const short indexS = matched->At(itrack);
311               if (indexS>-1)
312                 {       
313                   AliAODTrack *track = fAOD->GetTrack(indexS);
314                   ptot += track->P();
315                 }
316             }
317         
318           correction = ptot * fFractionHadronicCorrection;
319
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 );
324           }
325         
326           delete matched;
327         
328         }//end hadronic correction
329 */      
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);
335         }
336
337         Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
338
339         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
340         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
341         uArray->SetUnitTrackID(digitID);
342
343         Float_t unitEnergy = 0.;
344         Bool_t ok = kFALSE;
345         unitEnergy = uArray->GetUnitEnergy();
346
347         if(unitEnergy==0){
348           if(!fProcId){
349             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
350             fProcId = kTRUE;
351           }
352           fRefArray->Add(uArray);
353           fNDigitEmcal++;
354           ok = kTRUE;
355         }
356
357         // Detector flag
358         if(unitEnergy>0.)
359           uArray->SetUnitDetectorFlag(kAll);
360         else uArray->SetUnitDetectorFlag(kEmcal);
361       
362         uArray->SetUnitEnergy(unitEnergy+digitEt);
363
364         uArray->SetUnitCutFlag(kPtHigher);
365
366         // Signal or background
367         // To be modified !!!
368         uArray->SetUnitSignalFlag(kGood);
369
370         // This is for jet multiplicity
371         uArray->SetUnitClusterID(index);
372         
373         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
374
375       }// End loop over cells
376     } // End loop over clusters
377   } // end else
378
379   fNIn += goodDigit;
380
381   if(fDebug>1)
382     {
383       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
384       printf("goodDigit : %d\n", goodDigit);
385     }
386 }
387
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 
392 //  int size=0;
393 //  
394 //  for (Int_t itrk =  0; itrk <  fAOD->GetNumberOfTracks() ; itrk++) { //track loop
395 //  
396 //    AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;  
397 //    AliAODPid*    pid   = (AliAODPid*) track->GetDetPid();
398 //    
399 //    if(pid) {
400 //      Double_t emcpos[3];
401 //      pid->GetEMCALPosition(emcpos);      
402 //      TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
403 //
404 //      Double_t deta = tpos.Eta() - eta;
405 //      Double_t dphi = tpos.Phi() - phi;
406 //
407 //      if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
408 //      if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
409 //
410 //      Double_t res = sqrt(dphi*dphi + deta*deta);
411 //      
412 //      if (res< cut) {
413 //        //add this track-index
414 //        size++;
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);
418 //                    
419 //      }
420 //    }
421 //  }
422 //
423 //}
424
425 /*
426 //____________________________________________________________________________
427 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
428 {
429   // Set calibration parameters:
430   // if calibration database exists, they are read from database,
431   // otherwise, they are taken from digitizer.
432   //
433   // It is a user responsilibity to open CDB before reconstruction,
434   // for example:
435   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
436
437   //Check if calibration is stored in data base
438
439   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
440     {
441       AliCDBEntry *entry = (AliCDBEntry*)
442         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
443       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
444     }
445
446   if(!fCalibData)
447     printf("************* Calibration parameters not found in CDB! ****************");
448 //    AliFatal("Calibration parameters not found in CDB!");
449
450
451 }
452
453 //____________________________________________________________________________
454 Float_t  AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
455 {
456
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.
460
461   if(fCalibData){
462
463     if (fGeom==0)
464       printf("************* Did not get geometry from EMCALLoader ***************");
465 //      AliFatal("Did not get geometry from EMCALLoader") ;
466
467     Int_t iSupMod = -1;
468     Int_t nModule  = -1;
469     Int_t nIphi   = -1;
470     Int_t nIeta   = -1;
471     Int_t iphi    = -1;
472     Int_t ieta    = -1;
473
474     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
475     if(!bCell) {
476       // fGeom->PrintGeometry();
477       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
478       assert(0);
479     }
480
481     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
482
483     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
484     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
485
486    return -fADCpedestalECA + amp * fADCchannelECA ;
487
488   }
489   else //Return energy with default parameters if calibration is not available
490     return -fADCpedestalECA + amp * fADCchannelECA ;
491
492 }
493 */
494
495