Added methods to set/get centrailty from clusters in outer layer, which is less noisy...
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
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
61 ClassImp(AliCentralitySelectionTask)
62
63
64 //________________________________________________________________________
65 AliCentralitySelectionTask::AliCentralitySelectionTask():
66 AliAnalysisTaskSE(),
67   fDebug(0),
68   fAnalysisInput("ESD"),
69   fIsMCInput(kFALSE),
70   fFile(0),
71   fFile2(0),
72   fCentfilename(0),
73   fCentfilename2(0),
74   fCentV0M(0),
75   fCentFMD(0),
76   fCentTRK(0),
77   fCentTKL(0),
78   fCentCL0(0),
79   fCentCL1(0),
80   fCentV0MvsFMD(0),
81   fCentTKLvsV0M(0),
82   fCentZEMvsZDC(0),
83   fHtempV0M(0),
84   fHtempFMD(0),
85   fHtempTRK(0),
86   fHtempTKL(0),
87   fHtempCL0(0),
88   fHtempCL1(0),
89   fHtempV0MvsFMD(0),
90   fHtempTKLvsV0M(0),
91   fHtempZEMvsZDC(0)
92 {   
93   // Default constructor
94 }   
95
96 //________________________________________________________________________
97 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
98   AliAnalysisTaskSE(name),
99   fDebug(0),
100   fAnalysisInput("ESD"),
101   fIsMCInput(kFALSE),
102   fFile(0),
103   fFile2(0),
104   fCentfilename(0),
105   fCentfilename2(0),
106   fCentV0M(0),
107   fCentFMD(0),
108   fCentTRK(0),
109   fCentTKL(0),
110   fCentCL0(0),
111   fCentCL1(0),
112   fCentV0MvsFMD(0),
113   fCentTKLvsV0M(0),
114   fCentZEMvsZDC(0),
115   fHtempV0M(0),
116   fHtempFMD(0),
117   fHtempTRK(0),
118   fHtempTKL(0),
119   fHtempCL0(0),
120   fHtempCL1(0),
121   fHtempV0MvsFMD(0),
122   fHtempTKLvsV0M(0),
123   fHtempZEMvsZDC(0)
124 {
125   // Default constructor
126 }
127
128 //________________________________________________________________________
129 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
130 {
131   // Assignment operator
132   if (this!=&c) {
133     AliAnalysisTaskSE::operator=(c);
134   }
135   return *this;
136 }
137
138 //________________________________________________________________________
139 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
140   AliAnalysisTaskSE(ana),
141   fDebug(ana.fDebug),     
142   fAnalysisInput(ana.fDebug),
143   fIsMCInput(ana.fIsMCInput),
144   fFile(ana.fFile),
145   fFile2(ana.fFile2),
146   fCentfilename(ana.fCentfilename),
147   fCentfilename2(ana.fCentfilename2),
148   fCentV0M(ana.fCentV0M),
149   fCentFMD(ana.fCentFMD),
150   fCentTRK(ana.fCentTRK),
151   fCentTKL(ana.fCentTKL),
152   fCentCL0(ana.fCentCL0),
153   fCentCL1(ana.fCentCL1),
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),
162   fHtempCL1(ana.fHtempCL1),
163   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
164   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
165   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
166 {
167   // Copy Constructor   
168 }
169  
170 //________________________________________________________________________
171  AliCentralitySelectionTask::~AliCentralitySelectionTask()
172  {
173    // Destructor
174  }  
175
176 //________________________________________________________________________
177 void AliCentralitySelectionTask::UserCreateOutputObjects()
178 {  
179   // Create the output containers
180   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
181 }
182
183 //________________________________________________________________________
184 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
185
186   // Execute analysis for current event:
187   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
188   
189   Float_t  zncEnergy;          //  ZNC Energy
190   Float_t  zpcEnergy;          //  ZPC Energy
191   Float_t  znaEnergy;          //  ZNA Energy
192   Float_t  zpaEnergy;          //  ZPA Energy
193   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
194   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
195   
196   Int_t    nTracks = 0;             //  no. tracks
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
200
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
205
206   AliESDCentrality *esdCent = 0;
207
208   if(fAnalysisInput.CompareTo("ESD")==0){
209
210     AliVEvent* event = InputEvent();
211     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
212
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();
272     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
273     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
274     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
275     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
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
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]));
293   fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
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.));
298   
299   esdCent->SetCentralityV0M(fCentV0M);
300   esdCent->SetCentralityFMD(fCentFMD);
301   esdCent->SetCentralityTRK(fCentTRK);
302   esdCent->SetCentralityTKL(fCentTKL);
303   esdCent->SetCentralityCL0(fCentCL0);
304   esdCent->SetCentralityCL1(fCentCL1);
305   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
306   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
307   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
308 }
309
310 //________________________________________________________________________
311 void AliCentralitySelectionTask::ReadCentralityHistos() 
312 {
313   fFile  = new TFile(fCentfilename);
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")); 
319   fHtempCL1  = (TH1D*) (fFile->Get("hNclusters1_percentile")); 
320 }  
321   
322 void 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
330 void AliCentralitySelectionTask::SetPercentileFile(TString filename) 
331 {
332   fCentfilename = filename;
333   ReadCentralityHistos();
334 }
335
336 void AliCentralitySelectionTask::SetPercentileFile2(TString filename) 
337 {
338   fCentfilename2 = filename;
339   ReadCentralityHistos2();
340 }
341
342 //________________________________________________________________________
343 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
344 {
345   // Terminate analysis
346   fFile->Close();  
347   fFile2->Close();  
348 }
349 //________________________________________________________________________
350