f3f157fb408b47acc859e15094fd21394e165963
[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   fCentV0MvsFMD(0),
80   fCentTKLvsV0M(0),
81   fCentZEMvsZDC(0),
82   fHtempV0M(0),
83   fHtempFMD(0),
84   fHtempTRK(0),
85   fHtempTKL(0),
86   fHtempCL0(0),
87   fHtempV0MvsFMD(0),
88   fHtempTKLvsV0M(0),
89   fHtempZEMvsZDC(0)
90 {   
91   // Default constructor
92 }   
93
94 //________________________________________________________________________
95 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
96   AliAnalysisTaskSE(name),
97   fDebug(0),
98   fAnalysisInput("ESD"),
99   fIsMCInput(kFALSE),
100   fFile(0),
101   fFile2(0),
102   fCentfilename(0),
103   fCentfilename2(0),
104   fCentV0M(0),
105   fCentFMD(0),
106   fCentTRK(0),
107   fCentTKL(0),
108   fCentCL0(0),
109   fCentV0MvsFMD(0),
110   fCentTKLvsV0M(0),
111   fCentZEMvsZDC(0),
112   fHtempV0M(0),
113   fHtempFMD(0),
114   fHtempTRK(0),
115   fHtempTKL(0),
116   fHtempCL0(0),
117   fHtempV0MvsFMD(0),
118   fHtempTKLvsV0M(0),
119   fHtempZEMvsZDC(0)
120 {
121   // Default constructor
122 }
123
124 //________________________________________________________________________
125 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
126 {
127   // Assignment operator
128   if (this!=&c) {
129     AliAnalysisTaskSE::operator=(c);
130   }
131   return *this;
132 }
133
134 //________________________________________________________________________
135 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
136   AliAnalysisTaskSE(ana),
137   fDebug(ana.fDebug),     
138   fAnalysisInput(ana.fDebug),
139   fIsMCInput(ana.fIsMCInput),
140   fFile(ana.fFile),
141   fFile2(ana.fFile2),
142   fCentfilename(ana.fCentfilename),
143   fCentfilename2(ana.fCentfilename2),
144   fCentV0M(ana.fCentV0M),
145   fCentFMD(ana.fCentFMD),
146   fCentTRK(ana.fCentTRK),
147   fCentTKL(ana.fCentTKL),
148   fCentCL0(ana.fCentCL0),
149   fCentV0MvsFMD(ana.fCentV0MvsFMD),
150   fCentTKLvsV0M(ana.fCentTKLvsV0M),
151   fCentZEMvsZDC(ana.fCentZEMvsZDC),
152   fHtempV0M(ana.fHtempV0M),
153   fHtempFMD(ana.fHtempFMD),
154   fHtempTRK(ana.fHtempTRK),
155   fHtempTKL(ana.fHtempTKL),
156   fHtempCL0(ana.fHtempCL0),
157   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
158   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
159   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
160 {
161   // Copy Constructor   
162 }
163  
164 //________________________________________________________________________
165  AliCentralitySelectionTask::~AliCentralitySelectionTask()
166  {
167    // Destructor
168  }  
169
170 //________________________________________________________________________
171 void AliCentralitySelectionTask::UserCreateOutputObjects()
172 {  
173   // Create the output containers
174   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
175 }
176
177 //________________________________________________________________________
178 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
179
180   // Execute analysis for current event:
181   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
182   
183   Float_t  zncEnergy;               //  ZNC Energy
184   Float_t  zpcEnergy;               //  ZPC Energy
185   Float_t  znaEnergy;               //  ZNA Energy
186   Float_t  zpaEnergy;               //  ZPA Energy
187   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
188   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
189   
190   Int_t    nTracks    = 0;          //  no. tracks
191   Int_t    nTracklets = 0;          //  no. tracklets
192   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
193   Int_t    nChips[2];               //  no. chips on 2 SPD layers
194
195   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
196   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
197   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
198   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
199
200   AliESDCentrality *esdCent = 0;
201
202   if(fAnalysisInput.CompareTo("ESD")==0){
203
204     AliVEvent* event = InputEvent();
205     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
206
207     esdCent = esd->GetCentrality();
208
209     // ***** V0 info    
210     AliESDVZERO* esdV0 = esd->GetVZEROData();
211     multV0A=esdV0->GetMTotV0A();
212     multV0C=esdV0->GetMTotV0C();
213     
214     // ***** CB info (tracklets, clusters, chips)
215     nTracks    = event->GetNumberOfTracks();     
216
217     const AliMultiplicity *mult = esd->GetMultiplicity();
218
219     nTracklets = mult->GetNumberOfTracklets();
220
221     for(Int_t ilay=0; ilay<6; ilay++){
222       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
223     }
224
225     for(Int_t ilay=0; ilay<2; ilay++){
226       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
227     }
228     
229
230     // ***** FMD info
231     AliESDFMD *fmd = esd->GetFMDData();
232     Float_t totalMultA = 0;
233     Float_t totalMultC = 0;
234     const Float_t fFMDLowCut = 0.4;
235     
236     for(UShort_t det=1;det<=3;det++) {
237       Int_t nRings = (det==1 ? 1 : 2);
238       for (UShort_t ir = 0; ir < nRings; ir++) {          
239         Char_t   ring = (ir == 0 ? 'I' : 'O');
240         UShort_t nsec = (ir == 0 ? 20  : 40);
241         UShort_t nstr = (ir == 0 ? 512 : 256);
242         for(UShort_t sec =0; sec < nsec;  sec++)  {
243           for(UShort_t strip = 0; strip < nstr; strip++) {
244             
245             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
246             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
247             
248             Float_t nParticles=0;
249             
250             if(FMDmult > fFMDLowCut) {
251               nParticles = 1.;
252             }
253             
254             if (det<3) totalMultA = totalMultA + nParticles;
255             else totalMultC = totalMultC + nParticles;
256             
257           }
258         }
259       }
260     }
261     multFMDA = totalMultA;
262     multFMDC = totalMultC;
263     
264     // ***** ZDC info
265     AliESDZDC *esdZDC = esd->GetESDZDC();
266     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
267     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
268     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
269     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
270     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
271     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
272     
273   }   
274   else if(fAnalysisInput.CompareTo("AOD")==0){
275     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
276     // to be implemented
277     printf("  AOD analysis not yet implemented!!!\n\n");
278     return;
279   }
280
281   // ***** Centrality Selection
282   fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
283   fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
284   fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
285   fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
286   fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
287   
288   fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
289   fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
290   fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
291   
292   esdCent->SetCentralityV0M(fCentV0M);
293   esdCent->SetCentralityFMD(fCentFMD);
294   esdCent->SetCentralityTRK(fCentTRK);
295   esdCent->SetCentralityTKL(fCentTKL);
296   esdCent->SetCentralityCL0(fCentCL0);
297   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
298   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
299   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
300 }
301
302 //________________________________________________________________________
303 void AliCentralitySelectionTask::ReadCentralityHistos() 
304 {
305   fFile  = new TFile(fCentfilename);
306   fHtempV0M  = (TH1D*) (fFile->Get("hmultV0_percentile")); 
307   fHtempFMD  = (TH1D*) (fFile->Get("hmultFMD_percentile")); 
308   fHtempTRK  = (TH1D*) (fFile->Get("hNtracks_percentile")); 
309   fHtempTKL  = (TH1D*) (fFile->Get("hNtracklets_percentile")); 
310   fHtempCL0  = (TH1D*) (fFile->Get("hNclusters0_percentile")); 
311 }  
312   
313 void AliCentralitySelectionTask::ReadCentralityHistos2() 
314 {
315   fFile2  = new TFile(fCentfilename2);
316   fHtempV0MvsFMD  = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile")); 
317   fHtempTKLvsV0M  = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile")); 
318   fHtempZEMvsZDC  = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));  
319 }
320
321 void AliCentralitySelectionTask::SetPercentileFile(TString filename) 
322 {
323   fCentfilename = filename;
324   ReadCentralityHistos();
325 }
326
327 void AliCentralitySelectionTask::SetPercentileFile2(TString filename) 
328 {
329   fCentfilename2 = filename;
330   ReadCentralityHistos2();
331 }
332
333 //________________________________________________________________________
334 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
335 {
336   // Terminate analysis
337   fFile->Close();  
338   fFile2->Close();  
339 }
340 //________________________________________________________________________
341