Added methods to set/get centrailty from clusters in outer layer, which is less noisy...
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
CommitLineData
e8472cd0 1/**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//*****************************************************
17// Class AliCentralitySelectionTask
18// Class to analyze determine centrality
19// author: Alberica Toia
20//*****************************************************
21
22#include <TTree.h>
23#include <TList.h>
24#include <TH1F.h>
25#include <TH2F.h>
26#include <TProfile.h>
27#include <TFile.h>
28#include <TString.h>
29#include <TCanvas.h>
30#include <iostream>
31
32#include "AliAnalysisManager.h"
33#include "AliVEvent.h"
34#include "AliESD.h"
35#include "AliESDEvent.h"
36#include "AliESDHeader.h"
37#include "AliESDInputHandler.h"
38#include "AliESDZDC.h"
39#include "AliESDFMD.h"
40#include "AliESDVZERO.h"
41#include "AliESDCentrality.h"
42#include "AliMultiplicity.h"
43#include "AliAODHandler.h"
44#include "AliAODEvent.h"
45#include "AliAODVertex.h"
46#include "AliAODMCHeader.h"
47#include "AliMCEvent.h"
48#include "AliMCEventHandler.h"
49#include "AliMCParticle.h"
50#include "AliStack.h"
51#include "AliHeader.h"
52#include "AliAODMCParticle.h"
53#include "AliAnalysisTaskSE.h"
54#include "AliGenEventHeader.h"
55#include "AliGenHijingEventHeader.h"
56#include "AliPhysicsSelectionTask.h"
57#include "AliPhysicsSelection.h"
58#include "AliBackgroundSelection.h"
59#include "AliCentralitySelectionTask.h"
60
61ClassImp(AliCentralitySelectionTask)
62
63
64//________________________________________________________________________
65AliCentralitySelectionTask::AliCentralitySelectionTask():
66AliAnalysisTaskSE(),
67 fDebug(0),
68 fAnalysisInput("ESD"),
69 fIsMCInput(kFALSE),
70 fFile(0),
d15bf53f 71 fFile2(0),
e8472cd0 72 fCentfilename(0),
d15bf53f 73 fCentfilename2(0),
74 fCentV0M(0),
75 fCentFMD(0),
76 fCentTRK(0),
77 fCentTKL(0),
78 fCentCL0(0),
be0d4e9b 79 fCentCL1(0),
d15bf53f 80 fCentV0MvsFMD(0),
81 fCentTKLvsV0M(0),
82 fCentZEMvsZDC(0),
83 fHtempV0M(0),
84 fHtempFMD(0),
85 fHtempTRK(0),
86 fHtempTKL(0),
87 fHtempCL0(0),
be0d4e9b 88 fHtempCL1(0),
d15bf53f 89 fHtempV0MvsFMD(0),
90 fHtempTKLvsV0M(0),
91 fHtempZEMvsZDC(0)
e8472cd0 92{
93 // Default constructor
94}
95
96//________________________________________________________________________
97AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
98 AliAnalysisTaskSE(name),
99 fDebug(0),
100 fAnalysisInput("ESD"),
101 fIsMCInput(kFALSE),
102 fFile(0),
d15bf53f 103 fFile2(0),
e8472cd0 104 fCentfilename(0),
d15bf53f 105 fCentfilename2(0),
106 fCentV0M(0),
107 fCentFMD(0),
108 fCentTRK(0),
109 fCentTKL(0),
110 fCentCL0(0),
be0d4e9b 111 fCentCL1(0),
d15bf53f 112 fCentV0MvsFMD(0),
113 fCentTKLvsV0M(0),
114 fCentZEMvsZDC(0),
115 fHtempV0M(0),
116 fHtempFMD(0),
117 fHtempTRK(0),
118 fHtempTKL(0),
119 fHtempCL0(0),
be0d4e9b 120 fHtempCL1(0),
d15bf53f 121 fHtempV0MvsFMD(0),
122 fHtempTKLvsV0M(0),
123 fHtempZEMvsZDC(0)
e8472cd0 124{
125 // Default constructor
126}
127
128//________________________________________________________________________
129AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
130{
131 // Assignment operator
132 if (this!=&c) {
133 AliAnalysisTaskSE::operator=(c);
134 }
135 return *this;
136}
137
138//________________________________________________________________________
139AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
140 AliAnalysisTaskSE(ana),
141 fDebug(ana.fDebug),
142 fAnalysisInput(ana.fDebug),
143 fIsMCInput(ana.fIsMCInput),
144 fFile(ana.fFile),
d15bf53f 145 fFile2(ana.fFile2),
e8472cd0 146 fCentfilename(ana.fCentfilename),
d15bf53f 147 fCentfilename2(ana.fCentfilename2),
148 fCentV0M(ana.fCentV0M),
149 fCentFMD(ana.fCentFMD),
150 fCentTRK(ana.fCentTRK),
151 fCentTKL(ana.fCentTKL),
152 fCentCL0(ana.fCentCL0),
be0d4e9b 153 fCentCL1(ana.fCentCL1),
d15bf53f 154 fCentV0MvsFMD(ana.fCentV0MvsFMD),
155 fCentTKLvsV0M(ana.fCentTKLvsV0M),
156 fCentZEMvsZDC(ana.fCentZEMvsZDC),
157 fHtempV0M(ana.fHtempV0M),
158 fHtempFMD(ana.fHtempFMD),
159 fHtempTRK(ana.fHtempTRK),
160 fHtempTKL(ana.fHtempTKL),
161 fHtempCL0(ana.fHtempCL0),
be0d4e9b 162 fHtempCL1(ana.fHtempCL1),
d15bf53f 163 fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
164 fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
165 fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
e8472cd0 166{
167 // Copy Constructor
168}
169
170//________________________________________________________________________
171 AliCentralitySelectionTask::~AliCentralitySelectionTask()
172 {
173 // Destructor
174 }
175
176//________________________________________________________________________
177void AliCentralitySelectionTask::UserCreateOutputObjects()
178{
179 // Create the output containers
180 if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
181}
182
183//________________________________________________________________________
184void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
185{
186 // Execute analysis for current event:
187 if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
188
be0d4e9b 189 Float_t zncEnergy; // ZNC Energy
190 Float_t zpcEnergy; // ZPC Energy
191 Float_t znaEnergy; // ZNA Energy
192 Float_t zpaEnergy; // ZPA Energy
d15bf53f 193 Float_t zem1Energy = 0.; // ZEM1 Energy
194 Float_t zem2Energy = 0.; // ZEM2 Energy
e8472cd0 195
be0d4e9b 196 Int_t nTracks = 0; // no. tracks
d15bf53f 197 Int_t nTracklets = 0; // no. tracklets
198 Int_t nClusters[6]; // no. clusters on 6 ITS layers
199 Int_t nChips[2]; // no. chips on 2 SPD layers
e8472cd0 200
d15bf53f 201 Float_t multV0A = 0; // multiplicity from V0 reco side A
202 Float_t multV0C = 0; // multiplicity from V0 reco side C
203 Float_t multFMDA = 0; // multiplicity from FMD on detector A
204 Float_t multFMDC = 0; // multiplicity from FMD on detector C
e8472cd0 205
206 AliESDCentrality *esdCent = 0;
207
208 if(fAnalysisInput.CompareTo("ESD")==0){
209
210 AliVEvent* event = InputEvent();
211 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
d15bf53f 212
e8472cd0 213 esdCent = esd->GetCentrality();
214
215 // ***** V0 info
216 AliESDVZERO* esdV0 = esd->GetVZEROData();
217 multV0A=esdV0->GetMTotV0A();
218 multV0C=esdV0->GetMTotV0C();
219
220 // ***** CB info (tracklets, clusters, chips)
221 nTracks = event->GetNumberOfTracks();
222
223 const AliMultiplicity *mult = esd->GetMultiplicity();
224
225 nTracklets = mult->GetNumberOfTracklets();
226
227 for(Int_t ilay=0; ilay<6; ilay++){
228 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
229 }
230
231 for(Int_t ilay=0; ilay<2; ilay++){
232 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
233 }
234
235
236 // ***** FMD info
237 AliESDFMD *fmd = esd->GetFMDData();
238 Float_t totalMultA = 0;
239 Float_t totalMultC = 0;
240 const Float_t fFMDLowCut = 0.4;
241
242 for(UShort_t det=1;det<=3;det++) {
243 Int_t nRings = (det==1 ? 1 : 2);
244 for (UShort_t ir = 0; ir < nRings; ir++) {
245 Char_t ring = (ir == 0 ? 'I' : 'O');
246 UShort_t nsec = (ir == 0 ? 20 : 40);
247 UShort_t nstr = (ir == 0 ? 512 : 256);
248 for(UShort_t sec =0; sec < nsec; sec++) {
249 for(UShort_t strip = 0; strip < nstr; strip++) {
250
251 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
252 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
253
254 Float_t nParticles=0;
255
256 if(FMDmult > fFMDLowCut) {
257 nParticles = 1.;
258 }
259
260 if (det<3) totalMultA = totalMultA + nParticles;
261 else totalMultC = totalMultC + nParticles;
262
263 }
264 }
265 }
266 }
267 multFMDA = totalMultA;
268 multFMDC = totalMultC;
269
270 // ***** ZDC info
271 AliESDZDC *esdZDC = esd->GetESDZDC();
d15bf53f 272 zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
273 zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
274 znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
275 zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
e8472cd0 276 zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
277 zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
278
279 }
280 else if(fAnalysisInput.CompareTo("AOD")==0){
281 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
282 // to be implemented
283 printf(" AOD analysis not yet implemented!!!\n\n");
284 return;
285 }
286
d15bf53f 287 // ***** Centrality Selection
288 fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
289 fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
290 fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
291 fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
292 fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
be0d4e9b 293 fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
d15bf53f 294
295 fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
296 fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
297 fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
e8472cd0 298
d15bf53f 299 esdCent->SetCentralityV0M(fCentV0M);
300 esdCent->SetCentralityFMD(fCentFMD);
301 esdCent->SetCentralityTRK(fCentTRK);
302 esdCent->SetCentralityTKL(fCentTKL);
303 esdCent->SetCentralityCL0(fCentCL0);
be0d4e9b 304 esdCent->SetCentralityCL1(fCentCL1);
d15bf53f 305 esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
306 esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
307 esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
e8472cd0 308}
309
310//________________________________________________________________________
d15bf53f 311void AliCentralitySelectionTask::ReadCentralityHistos()
e8472cd0 312{
e8472cd0 313 fFile = new TFile(fCentfilename);
d15bf53f 314 fHtempV0M = (TH1D*) (fFile->Get("hmultV0_percentile"));
315 fHtempFMD = (TH1D*) (fFile->Get("hmultFMD_percentile"));
316 fHtempTRK = (TH1D*) (fFile->Get("hNtracks_percentile"));
317 fHtempTKL = (TH1D*) (fFile->Get("hNtracklets_percentile"));
318 fHtempCL0 = (TH1D*) (fFile->Get("hNclusters0_percentile"));
be0d4e9b 319 fHtempCL1 = (TH1D*) (fFile->Get("hNclusters1_percentile"));
d15bf53f 320}
e8472cd0 321
d15bf53f 322void AliCentralitySelectionTask::ReadCentralityHistos2()
323{
324 fFile2 = new TFile(fCentfilename2);
325 fHtempV0MvsFMD = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
326 fHtempTKLvsV0M = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
327 fHtempZEMvsZDC = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
328}
329
330void AliCentralitySelectionTask::SetPercentileFile(TString filename)
331{
332 fCentfilename = filename;
333 ReadCentralityHistos();
334}
335
336void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
337{
338 fCentfilename2 = filename;
339 ReadCentralityHistos2();
e8472cd0 340}
341
342//________________________________________________________________________
343void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
344{
345 // Terminate analysis
346 fFile->Close();
d15bf53f 347 fFile2->Close();
e8472cd0 348}
349//________________________________________________________________________
350