Changes in the ACORDE libraries to compile on Windows/Cygwin
[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"
ee7de0dd 48
49ClassImp(AliJetFillUnitArrayEMCalDigits)
50
51//_____________________________________________________________________________
52AliJetFillUnitArrayEMCalDigits::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//_____________________________________________________________________________
81AliJetFillUnitArrayEMCalDigits::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//____________________________________________________________________________
111void 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//_____________________________________________________________________________
126AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
127{
128 // destructor
129}
130
131//_____________________________________________________________________________
132void 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
8ada0ffe 169 if (type != AliESDCaloCluster::kEMCALClusterv1) continue;
ee7de0dd 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
90dc70f0 178 TArrayS *digID = fClus->GetDigitIndex();
179 TArrayS *digEnergy = fClus->GetDigitAmplitude();
180 Float_t *digitEnergy = new Float_t[nD];
ee7de0dd 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
90dc70f0 187 Int_t idF = (Int_t)digID->At(k);
ee7de0dd 188 // Calibration for an energy in GeV
90dc70f0 189 digitEnergy[k] = (Float_t)digEnergy->At(k)/500.;
ee7de0dd 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//_____________________________________________________________________________
239Float_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}