Clean-up in includes.
[u/mrichter/AliRoot.git] / JETAN / AliJetFillUnitArrayEMCalDigits.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 with EMCal information 
19 // Called by ESD reader for jet analysis
20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
21 //------------------------------------------------------------- 
22
23 // --- Standard library ---
24 #include <Riostream.h>
25
26 // --- ROOT system ---
27 #include <TSystem.h>
28 #include <TLorentzVector.h>
29 #include <TRefArray.h> 
30 #include "TTask.h"
31 #include <TGeoManager.h>
32 #include <TMath.h>
33 #include <TClonesArray.h>
34 #include <TArrayS.h>
35
36 // --- AliRoot header files ---
37 #include "AliJetFinder.h"
38 #include "AliJetReaderHeader.h"
39 #include "AliJetReader.h"
40 #include "AliJetESDReader.h"
41 #include "AliJetESDReaderHeader.h"
42 //#include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliJetDummyGeo.h"
45 #include "AliESDCaloCluster.h"
46 #include "AliJetUnitArray.h"
47 #include "AliJetFillUnitArrayEMCalDigits.h"
48
49 ClassImp(AliJetFillUnitArrayEMCalDigits)
50
51 //_____________________________________________________________________________
52 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
53   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
54     fESD(0),    
55     fNumUnits(0), 
56     fEtaMinCal(0.),
57     fEtaMaxCal(0.),
58     fPhiMinCal(0.),
59     fPhiMaxCal(0.),
60     fNIn(0),
61     fOpt(0),
62     fDebug(0),
63     fNCEMCAL(0),
64     fNCPHOS(0),
65     fNCCalo(0),
66     fTPCGrid(0x0),
67     fEMCalGrid(0x0),
68     fReaderHeader(0x0),
69     fMomentumArray(0x0),
70     fUnitArray(0x0),
71     fRefArray(0x0),
72     fGeom(0x0),
73     fClus(0x0),
74     fNDigitEmcal(0),
75     fNDigitEmcalCut(0)
76 {
77   // constructor
78 }
79
80 //_____________________________________________________________________________
81 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/)
82   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
83     fESD(0),    
84     fNumUnits(0), 
85     fEtaMinCal(0.),
86     fEtaMaxCal(0.),
87     fPhiMinCal(0.),
88     fPhiMaxCal(0.),
89     fNIn(0),
90     fOpt(0),
91     fDebug(0),
92     fNCEMCAL(0),
93     fNCPHOS(0),
94     fNCCalo(0),
95     fTPCGrid(0x0),
96     fEMCalGrid(0x0),
97     fReaderHeader(0x0),
98     fMomentumArray(0x0),
99     fUnitArray(0x0),
100     fRefArray(0x0),
101     fGeom(0x0),
102     fClus(0x0),
103     fNDigitEmcal(0),
104     fNDigitEmcalCut(0)
105 {
106   // constructor
107 }
108
109
110 //____________________________________________________________________________
111 void AliJetFillUnitArrayEMCalDigits::InitParameters()
112 {
113   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
114
115   fEtaMinCal = fGeom->GetArm1EtaMin();
116   fEtaMaxCal = fGeom->GetArm1EtaMax();
117   fPhiMinCal = fGeom->GetArm1PhiMin();
118   fPhiMaxCal = fGeom->GetArm1PhiMax(); 
119   fClus      = 0;
120
121   if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
122
123 }
124
125 //_____________________________________________________________________________
126 AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
127 {
128   // destructor
129 }
130
131 //_____________________________________________________________________________
132 void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
133 {
134   //
135   // Main method.
136   // Explain
137
138   fDebug = fReaderHeader->GetDebug();
139   fOpt = fReaderHeader->GetDetector();
140
141   // Init parameters
142   InitParameters();
143
144   // Get number of clusters from EMCAL
145   
146   Int_t   nDigitTot      = 0;
147   Int_t   goodDigit      = 0;
148   Int_t   beg            = 0;
149   Int_t   end            = 0;
150   Float_t ptMin          = fReaderHeader->GetPtCut();
151
152   // Loop over calo clusters 
153   //------------------------------------------------------------------
154   Int_t type = 0; 
155   Int_t index = 0;
156
157   // Total number of EMCAL cluster
158   end =  fESD->GetNumberOfCaloClusters();
159
160   for(Int_t j = beg; j < end; j++) {
161       fClus = fESD->GetCaloCluster(j);
162       if(!fClus->IsEMCAL()) continue;
163
164       type = fClus->GetClusterType(); 
165       index = fClus->GetID();
166       nDigitTot = fClus->GetNumberOfDigits();
167       
168       // Keep clusters or pseudo clusters
169       if (type != AliESDCaloCluster::kEMCALClusterv1) continue;
170       //      if (type != AliESDCaloCluster::kPseudoCluster) continue;
171
172       // Get the digit index and the digit information
173       //============================================================
174           
175       // Get number of digits in a cluster
176       Int_t nD = fClus->GetNumberOfDigits();
177
178       TArrayS *digID = fClus->GetDigitIndex();      
179       TArrayS *digEnergy = fClus->GetDigitAmplitude();
180       Float_t *digitEnergy = new Float_t[nD];      
181       //      Float_t digitEn = 0.;
182           
183       // Loop over digits
184       for(Int_t k=0; k<nD; k++) {
185         
186         // Convert energy in GeV
187         Int_t idF = (Int_t)digID->At(k);
188         // Calibration for an energy in GeV
189         digitEnergy[k] = (Float_t)digEnergy->At(k)/500.; 
190
191         // Second method to extract eta, phi positions of a digit
192         //=================================================================
193
194         Float_t etaD=-10., phiD=-10.;
195         fGeom->EtaPhiFromIndex(idF,etaD,phiD); 
196         //          fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD);
197         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
198
199         Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
200
201         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idF);
202         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
203
204         Float_t unitEnergy = 0.;
205         Bool_t ok = kFALSE;
206         unitEnergy = uArray->GetUnitEnergy();
207         if(unitEnergy==0){
208           fRefArray->Add(uArray);
209           fNDigitEmcal++;
210           ok = kTRUE;
211         }
212         uArray->SetUnitEnergy(unitEnergy+etDigit);
213         // Put a cut flag
214         if(uArray->GetUnitEnergy()<ptMin)
215           uArray->SetUnitCutFlag(kPtSmaller);
216         else {
217           uArray->SetUnitCutFlag(kPtHigher);
218           if(ok) fNDigitEmcalCut++;
219         }
220         // Detector flag
221         if(unitEnergy>0.) 
222           uArray->SetUnitDetectorFlag(kAll);
223         else uArray->SetUnitDetectorFlag(kEmcal);
224         
225         // This is for jet multiplicity
226         uArray->SetUnitClusterID(index);
227         
228         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
229         
230       } // End loop over digits
231       
232   } // End loop over clusters
233   
234   fNIn += goodDigit;
235
236 }
237
238 //_____________________________________________________________________________
239 Float_t  AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
240 {
241   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
242   return 2.*atan(exp(-arg));
243
244
245 }