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