8b8ccb8b57a6cbbb819d0752e69acb55e0ecacce
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraAODEventCuts.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-2009, 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 //         AliSpectraAODEventCuts class
19 //-----------------------------------------------------------------
20
21 #include "TChain.h"
22 #include "TTree.h"
23 #include "TLegend.h"
24 #include "TH1F.h"
25 #include "TH1I.h"
26 #include "TH2F.h"
27 #include "TCanvas.h"
28 #include "TProfile.h"
29 #include "TSpline.h"
30 #include "AliAnalysisTask.h"
31 #include "AliAnalysisManager.h"
32 #include "AliAODTrack.h"
33 #include "AliAODMCParticle.h"
34 #include "AliAODEvent.h"
35 #include "AliAODInputHandler.h"
36 #include "AliAnalysisTaskESDfilter.h"
37 #include "AliAnalysisDataContainer.h"
38 #include "AliESDVZERO.h"
39 #include "AliAODVZERO.h"
40 #include "AliSpectraAODEventCuts.h"
41 #include "AliSpectraAODTrackCuts.h"
42 #include <iostream>
43
44 using namespace std;
45
46 ClassImp(AliSpectraAODEventCuts)
47
48 AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) : 
49 TNamed(name, "AOD Event Cuts"),
50 fAOD(0),
51 fSelectBit(AliVEvent::kMB),
52 fCentralityMethod("V0M"),
53 fTrackBits(1),
54 fIsMC(0),
55 fIsLHC10h(1),
56 fTrackCuts(0),
57 fIsSelected(0),
58 fCentralityCutMin(0.),
59 fCentralityCutMax(999),
60 fQVectorCutMin(-999.),
61 fQVectorCutMax(999.),
62 fVertexCutMin(-10.),
63 fVertexCutMax(10.),
64 fMultiplicityCutMin(-999.),
65 fMultiplicityCutMax(99999.),
66 fqTPC(-999.),
67 fqV0C(-999.),
68 fqV0A(-999.),
69 fqV0Cx(-999.),
70 fqV0Ax(-999.),
71 fqV0Cy(-999.),
72 fqV0Ay(-999.),
73 fPsiV0C(-999.),
74 fPsiV0A(-999.),
75 fCent(-999.),
76 fOutput(0),
77 fCalib(0),
78 fRun(-1),
79 fMultV0(0),
80 fV0Cpol1(-1),
81 fV0Cpol2(-1),
82 fV0Cpol3(-1),
83 fV0Cpol4(-1),
84 fV0Apol1(-1),
85 fV0Apol2(-1),
86 fV0Apol3(-1),
87 fV0Apol4(-1),
88 fQvecIntList(0),
89 fQvecIntegral(0),
90 fSplineArrayV0A(0),
91 fSplineArrayV0C(0),
92 fSplineArrayTPC(0),
93 fQgenIntegral(0),
94 fSplineArrayV0Agen(0),
95 fSplineArrayV0Cgen(0),
96 fQvecMC(0),
97 fNch(0),
98 fQvecCalibType(0),
99 fV0Aeff(0)
100 {
101         // Constructor
102         fOutput=new TList();
103         fOutput->SetOwner();
104         fOutput->SetName("fOutput");
105
106         fCalib=new TList();
107         fCalib->SetOwner();
108         fCalib->SetName("fCalib");
109
110         fQvecIntList=new TList();
111         fQvecIntList->SetOwner();
112         fQvecIntList->SetName("fQvecIntList");
113
114         TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
115         TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
116         TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
117         TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
118         TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
119         TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
120         TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
121         TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
122         TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
123         TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
124         TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
125         TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
126         TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
127         TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
128         TH2F *fV0Mmc = new TH2F("fV0Mmc", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
129
130         fSplineArrayV0A = new TObjArray();
131         fSplineArrayV0A->SetOwner();
132         fSplineArrayV0C = new TObjArray();
133         fSplineArrayV0C->SetOwner();
134         fSplineArrayTPC = new TObjArray();
135         fSplineArrayTPC->SetOwner();
136         fSplineArrayV0Agen = new TObjArray();
137         fSplineArrayV0Agen->SetOwner();
138         fSplineArrayV0Cgen = new TObjArray();
139         fSplineArrayV0Cgen->SetOwner();
140
141         fOutput->Add(fHistoCuts);
142         fOutput->Add(fHistoVtxBefSel);
143         fOutput->Add(fHistoVtxAftSel);
144         fOutput->Add(fHistoEtaBefSel);
145         fOutput->Add(fHistoEtaAftSel);
146         fOutput->Add(fHistoNChAftSel);
147         fOutput->Add(fHistoQVector);
148         fOutput->Add(fHistoEP);
149         fOutput->Add(fPsiACor);
150         fOutput->Add(fPsiCCor);
151         fOutput->Add(fQVecACor);
152         fOutput->Add(fQVecCCor);
153         fOutput->Add(fV0M);
154         fOutput->Add(fV0MCor);
155         fOutput->Add(fV0Mmc);
156
157         for (Int_t i = 0; i<10; i++){
158                 fMeanQxa2[i] = -1;
159                 fMeanQya2[i] = -1;
160                 fMeanQxc2[i] = -1;
161                 fMeanQyc2[i] = -1;
162         }
163 }
164
165 //______________________________________________________
166 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
167 {
168         // Returns true if Event Cuts are selected and applied
169         fAOD = aod;
170         fTrackCuts = trackcuts; // FIXME: if track cuts is 0, do not set (use the pre-set member). Do we need to pass this here at all??
171         // FIXME: all those references by name are slow.
172         ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
173         Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
174         if(!IsPhysSel)return IsPhysSel;
175         ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
176         //loop on tracks, before event selection, filling QA histos
177         AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
178         if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
179         fIsSelected =kFALSE;
180         if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
181                 fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
182         }
183         if(fIsSelected){
184                 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
185                 if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
186         }
187         Int_t Nch=0;
188         for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
189                 AliAODTrack* track = fAOD->GetTrack(iTracks);
190                 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
191                 ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
192                 if(fIsSelected){
193                         ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
194                         Nch++;
195                 }
196         }
197         fNch = Nch;
198         if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
199         return fIsSelected;
200 }
201
202 //______________________________________________________
203 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
204 {
205         // reject events outside of range
206         AliAODVertex * vertex = fAOD->GetPrimaryVertex();
207         //when moving to 2011 w├Čone has to add a cut using SPD vertex.
208         //The point is that for events with |z|>20 the vertexer tracks is not working (only 2011!). One has to put a safety cut using SPD vertex large e.g. 15cm
209         if (!vertex)
210         {
211                 ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
212                 return kFALSE;
213         }
214         if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
215         {
216                 return kTRUE;
217         }
218         ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
219         return kFALSE;
220 }
221
222 //______________________________________________________
223 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
224 {
225         // Check centrality cut
226         fCent=-999.;
227         fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
228         if ( (fCent <= fCentralityCutMax)  &&  (fCent >= fCentralityCutMin) )  return kTRUE;
229         ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
230         return kFALSE;
231 }
232
233 //______________________________________________________
234 Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
235 {
236         // Check multiplicity cut
237         // FIXME: why this is not tracket in the track stats histos?
238         Int_t Ncharged=0;
239         for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
240                 AliAODTrack* track = fAOD->GetTrack(iTracks);
241                 if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
242                 Ncharged++;
243         }
244         if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
245
246         return kFALSE;
247 }
248
249 //______________________________________________________
250
251 Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
252
253         Double_t qVZERO = -999.;
254         Double_t psi=-999.;
255
256         if(fIsLHC10h){
257                 qVZERO=CalculateQVectorLHC10h();
258                 psi=fPsiV0A;
259         }else{
260                 qVZERO=CalculateQVector();
261                 psi=fPsiV0A;
262         }
263         //q vector from TPC
264         fqTPC=CalculateQVectorTPC();
265
266         //cut on q vector
267         if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
268         Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
269         ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
270         ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
271         ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
272
273
274         return kTRUE;
275 }
276
277 //______________________________________________________
278 Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
279
280         Int_t run = fAOD->GetRunNumber();
281         if(run != fRun){
282                 // Load the calibrations run dependent
283                 if(OpenInfoCalbration(run))fRun=run;
284                 else{
285                         fqV0C=-999.;
286                         fqV0A=-999.;
287                         return -999.;
288                 }
289         }
290
291         //V0 info
292         Double_t Qxa2 = 0, Qya2 = 0;
293         Double_t Qxc2 = 0, Qyc2 = 0;
294         Double_t sumMc = 0, sumMa = 0;
295
296         AliAODVZERO* aodV0 = fAOD->GetVZEROData();
297
298         for (Int_t iv0 = 0; iv0 < 64; iv0++) {
299
300                 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
301
302                 Float_t multv0 = aodV0->GetMultiplicity(iv0);
303                 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
304
305                 if (iv0 < 32){
306
307                         Double_t multCorC = -10;
308
309                         if (iv0 < 8)
310                                 multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
311                         else if (iv0 >= 8 && iv0 < 16)
312                                 multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
313                         else if (iv0 >= 16 && iv0 < 24)
314                                 multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
315                         else if (iv0 >= 24 && iv0 < 32)
316                                 multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
317
318                         if (multCorC < 0){
319                                 cout<<"Problem with multiplicity in V0C"<<endl;
320                                 fqV0C=-999.;
321                                 fqV0A=-999.;
322                                 return -999.;
323                         }
324
325                         Qxc2 += TMath::Cos(2*phiV0) * multCorC;
326                         Qyc2 += TMath::Sin(2*phiV0) * multCorC;
327
328                         sumMc += multCorC;
329                         ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
330
331                 } else {
332
333                         Double_t multCorA = -10;
334
335                         if (iv0 >= 32 && iv0 < 40)
336                                 multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
337                         else if (iv0 >= 40 && iv0 < 48)
338                                 multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
339                         else if (iv0 >= 48 && iv0 < 56)
340                                 multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
341                         else if (iv0 >= 56 && iv0 < 64)
342                                 multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
343
344                         if (multCorA < 0){
345                                 cout<<"Problem with multiplicity in V0A"<<endl;
346                                 fqV0C=-999.;
347                                 fqV0A=-999.;
348                                 return -999.;
349                         }
350
351                         Qxa2 += TMath::Cos(2*phiV0) * multCorA;
352                         Qya2 += TMath::Sin(2*phiV0) * multCorA;
353
354                         sumMa += multCorA;
355                         ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
356                 }
357         }
358
359         Short_t centrV0  = GetCentrCode(fAOD);
360
361         Double_t Qxamean2 = 0.;
362         Double_t Qyamean2 = 0.;
363         Double_t Qxcmean2 = 0.;
364         Double_t Qycmean2 = 0.;
365
366         if(centrV0!=-1){
367                 Qxamean2 = fMeanQxa2[centrV0];
368                 Qyamean2 = fMeanQya2[centrV0];
369                 Qxcmean2 = fMeanQxc2[centrV0];
370                 Qycmean2 = fMeanQyc2[centrV0];
371         }
372
373         Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
374         Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
375         Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
376         Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
377
378         fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
379         fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
380
381         ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
382         ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
383
384         fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
385         fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
386         fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
387         fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
388         fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
389         fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
390
391         ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
392         ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
393
394         return fqV0A; //FIXME we have to combine VZERO-A and C
395 }
396
397 //______________________________________________________
398
399 Double_t AliSpectraAODEventCuts::CalculateQVectorTPC(Double_t etaMin,Double_t etaMax){
400
401         Double_t Qx2 = 0, Qy2 = 0;
402         Int_t mult = 0;
403         for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
404                 AliAODTrack* aodTrack = fAOD->GetTrack(iT);
405                 if (!aodTrack->TestFilterBit(128)) continue;  //FIXME track type hard coded -> TPC only constrained to the vertex
406                 if (aodTrack->Eta() <etaMin || aodTrack->Eta() > etaMax)continue;
407                 mult++;
408                 Qx2 += TMath::Cos(2*aodTrack->Phi());
409                 Qy2 += TMath::Sin(2*aodTrack->Phi());
410         }
411         if(mult!=0)fqTPC= TMath::Sqrt((Qx2*Qx2 + Qy2*Qy2)/mult);
412         return fqTPC;
413 }
414
415 //______________________________________________________
416
417 Double_t AliSpectraAODEventCuts::CalculateQVector(){
418
419         //V0 info
420         Double_t Qxa2 = 0, Qya2 = 0;
421         Double_t Qxc2 = 0, Qyc2 = 0;
422
423         AliAODVZERO* aodV0 = fAOD->GetVZEROData();
424
425         for (Int_t iv0 = 0; iv0 < 64; iv0++) {
426
427                 Float_t multv0 = aodV0->GetMultiplicity(iv0);
428
429                 ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
430
431         }
432
433         fPsiV0A = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,8,2,Qxa2,Qya2); // V0A
434         fPsiV0C = fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,9,2,Qxc2,Qyc2); // V0C
435
436         ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
437         ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
438
439         fqV0A = TMath::Sqrt((Qxa2*Qxa2 + Qya2*Qya2));
440         fqV0C = TMath::Sqrt((Qxc2*Qxc2 + Qyc2*Qyc2));
441         fqV0Ax = Qxa2;
442         fqV0Cx = Qxc2;
443         fqV0Ay = Qya2;
444         fqV0Cy = Qyc2;
445
446         ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
447         ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
448
449         return fqV0A; //FIXME we have to combine VZERO-A and C
450
451 }
452
453 //______________________________________________________
454 Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
455 {
456
457         Short_t centrCode = -1;
458
459         AliCentrality* centrality = 0;
460         AliAODEvent* aod = (AliAODEvent*)ev;
461         centrality = aod->GetHeader()->GetCentralityP();
462
463         Float_t centV0 = centrality->GetCentralityPercentile("V0M");
464         Float_t centTrk = centrality->GetCentralityPercentile("TRK");
465
466
467         if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){
468
469                 if ((centV0 > 0) && (centV0 <= 5.0))
470                         centrCode = 0;
471                 else if ((centV0 > 5.0) && (centV0 <= 10.0))
472                         centrCode = 1;
473                 else if ((centV0 > 10.0) && (centV0 <= 20.0))
474                         centrCode = 2;
475                 else if ((centV0 > 20.0) && (centV0 <= 30.0))
476                         centrCode = 3;
477                 else if ((centV0 > 30.0) && (centV0 <= 40.0))
478                         centrCode = 4;
479                 else if ((centV0 > 40.0) && (centV0 <= 50.0))
480                         centrCode = 5;
481                 else if ((centV0 > 50.0) && (centV0 <= 60.0))
482                         centrCode = 6;
483                 else if ((centV0 > 60.0) && (centV0 <= 70.0))
484                         centrCode = 7;
485                 else if ((centV0 > 70.0) && (centV0 <= 80.0))
486                         centrCode = 8;
487         }
488
489         return centrCode;
490
491 }
492
493 //______________________________________________________
494 void AliSpectraAODEventCuts::PrintCuts()
495 {
496         // print info about event cuts
497         cout << "Event Stats" << endl;
498         cout << " > Trigger Selection: " << fSelectBit << endl;
499         cout << " > Centrality estimator: " << fCentralityMethod << endl;
500         cout << " > Number of accepted events: " << NumberOfEvents() << endl;
501         cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
502         cout << " > Number of PhysSel events: " <<  NumberOfPhysSelEvents()  << endl;
503         cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
504         cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
505         cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
506         cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
507         cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
508         cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
509         cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
510         cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
511         cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
512 }
513
514 //_____________________________________________________________________________
515 Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
516 {
517
518         AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
519         if(!cont){
520                 printf("OADB object hMultV0BefCorr is not available in the file\n");
521                 return 0;
522         }
523
524         if(!(cont->GetObject(run))){
525                 printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
526                 return 0;
527         }
528         fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
529
530         TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7);
531         fMultV0->Fit(fpolc1, "RN");
532         fV0Cpol1 = fpolc1->GetParameter(0);
533
534         TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15);
535         fMultV0->Fit(fpolc2, "RN");
536         fV0Cpol2 = fpolc2->GetParameter(0);
537
538         TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23);
539         fMultV0->Fit(fpolc3, "RN");
540         fV0Cpol3 = fpolc3->GetParameter(0);
541
542         TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31);
543         fMultV0->Fit(fpolc4, "RN");
544         fV0Cpol4 = fpolc4->GetParameter(0);
545
546         TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
547         fMultV0->Fit(fpola1, "RN");
548         fV0Apol1 = fpola1->GetParameter(0);
549
550         TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
551         fMultV0->Fit(fpola2, "RN");
552         fV0Apol2 = fpola2->GetParameter(0);
553
554         TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
555         fMultV0->Fit(fpola3, "RN");
556         fV0Apol3 = fpola3->GetParameter(0);
557
558         TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
559         fMultV0->Fit(fpola4, "RN");
560         fV0Apol4 = fpola4->GetParameter(0);
561
562         for(Int_t i=0; i < 10; i++){
563
564                 char nameQxa2[100];
565                 snprintf(nameQxa2,100, "hQxa2m_%i", i);
566
567                 char nameQya2[100];
568                 snprintf(nameQya2,100, "hQya2m_%i", i);
569
570                 char nameQxc2[100];
571                 snprintf(nameQxc2,100, "hQxc2m_%i", i);
572
573                 char nameQyc2[100];
574                 snprintf(nameQyc2,100, "hQyc2m_%i", i);
575
576                 AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
577                 if(!contQxa2){
578                         printf("OADB object %s is not available in the file\n", nameQxa2);
579                         return 0;
580                 }
581
582                 if(!(contQxa2->GetObject(run))){
583                         printf("OADB object %s is not available for run %i\n", nameQxa2, run);
584                         return 0;
585                 }
586
587                 fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
588
589
590                 AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
591                 if(!contQya2){
592                         printf("OADB object %s is not available in the file\n", nameQya2);
593                         return 0;
594                 }
595
596                 if(!(contQya2->GetObject(run))){
597                         printf("OADB object %s is not available for run %i\n", nameQya2, run);
598                         return 0;
599                 }
600
601                 fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
602
603
604                 AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
605                 if(!contQxc2){
606                         printf("OADB object %s is not available in the file\n", nameQxc2);
607                         return 0;
608                 }
609
610                 if(!(contQxc2->GetObject(run))){
611                         printf("OADB object %s is not available for run %i\n", nameQxc2, run);
612                         return 0;
613                 }
614
615                 fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
616
617
618                 AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
619                 if(!contQyc2){
620                         printf("OADB object %s is not available in the file\n", nameQyc2);
621                         return 0;
622                 }
623
624                 if(!(contQyc2->GetObject(run))){
625                         printf("OADB object %s is not available for run %i\n", nameQyc2, run);
626                         return 0;
627                 }
628
629                 fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
630
631         }
632         return 1;
633 }
634
635 //______________________________________________________
636
637
638 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
639 {
640         // Merge a list of AliSpectraAODEventCuts objects with this.
641         // Returns the number of merged objects (including this).
642
643         AliInfo("Merging");
644
645         if (!list)
646                 return 0;
647
648         if (list->IsEmpty())
649                 return 1;
650
651         TIterator* iter = list->MakeIterator();
652         TObject* obj;
653
654         // collections of all histograms
655         TList collections;
656
657         Int_t count = 0;
658
659         while ((obj = iter->Next())) {
660                 AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
661                 if (entry == 0)
662                         continue;
663
664                 TList * l = entry->GetOutputList();
665                 collections.Add(l);
666                 count++;
667         }
668
669         fOutput->Merge(&collections);
670
671         delete iter;
672
673         return count+1;
674 }
675
676 //______________________________________________________
677 Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
678
679         //check Qvec and Centrality consistency
680         if(fCent>90.) return -999.;
681         if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
682         if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
683         if(v0side==2/*TPC*/){ if(fqTPC== -999.) return -999.; }
684
685         fQvecIntegral = 0x0;
686
687         if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
688         if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
689         if(v0side==2/*TPC*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("TPC"); }
690
691
692         Int_t ic = -999;
693
694         if(fQvecCalibType==1){
695                 if(fNch<0.) return -999.;
696                 ic = GetNchBin(fQvecIntegral);
697         } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
698
699         TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
700
701         TSpline *spline = 0x0;
702
703         if(v0side==0/*V0A*/){
704                 if( CheckSplineArray(fSplineArrayV0A, ic) ) {
705                         spline = (TSpline*)fSplineArrayV0A->At(ic);
706                         //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
707                 } else {
708                         spline = new TSpline3(h1D,"sp3");
709                         fSplineArrayV0A->AddAtAndExpand(spline,ic);
710                         //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
711                 }
712         }
713         else if(v0side==1/*V0C*/){
714                 if( CheckSplineArray(fSplineArrayV0C, ic) ) {
715                         spline = (TSpline*)fSplineArrayV0C->At(ic);
716                 } else {
717                         spline = new TSpline3(h1D,"sp3");
718                         fSplineArrayV0C->AddAtAndExpand(spline,ic);
719                 }
720         }
721         else if(v0side==2/*TPC*/){
722                 if( CheckSplineArray(fSplineArrayTPC, ic) ) {
723                         spline = (TSpline*)fSplineArrayTPC->At(ic);
724                 } else {
725                         spline = new TSpline3(h1D,"sp3");
726                         fSplineArrayTPC->AddAtAndExpand(spline,ic);
727                 }
728         }
729
730         Double_t percentile =  -999.;
731         if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
732         if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
733
734         if(percentile>100. || percentile<0.) return -999.;
735
736         return percentile;
737 }
738
739 //______________________________________________________
740 Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side, Int_t type=1){
741
742         if(!fIsMC) return -999.;
743
744         // V0A efficiecy
745         if(type==1){
746                 fV0Aeff = 0x0;
747                 fV0Aeff = (TH1F*)fQvecIntList->FindObject("vzeroa_efficiency_prim_plus_sec_over_gen");
748                 if(!fV0Aeff) return -999.;
749         }
750
751         // get MC array
752         TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
753
754         if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
755
756         Double_t Qx2mc = 0., Qy2mc = 0.;
757         Double_t mult2mc = 0;
758
759         Int_t nMC = arrayMC->GetEntries();
760
761         if(type==0){ // type==0: q-vec from tracks in vzero acceptance
762
763                 // loop on generated
764                 for (Int_t iMC = 0; iMC < nMC; iMC++){
765                         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
766                         if(!partMC->Charge()) continue;//Skip neutrals
767
768                         // check vzero side
769                         if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
770
771                         // Calculate Qvec components
772                         Qx2mc += TMath::Cos(2.*partMC->Phi());
773                         Qy2mc += TMath::Sin(2.*partMC->Phi());
774                         mult2mc++;
775
776                 }// end loop on generated.
777
778         }//end if on type==0
779
780         else if(type==1){ // type==1 (default): q-vec from vzero
781
782                 // only used in qgen_vzero
783                 Double_t multv0mc[64];
784                 for(Int_t iCh=0;iCh<64;iCh++) multv0mc[iCh]=0; // initialize multv0mc to zero
785
786                 // loop on generated
787                 for (Int_t iMC = 0; iMC < nMC; iMC++){
788
789                         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
790                         if(!partMC->Charge()) continue;//Skip neutrals
791
792                         // check vzero side
793                         if( CheckVZEROacceptance(partMC->Eta()) != v0side ) continue;
794
795                         //get v0 channel of gen track
796                         Int_t iv0 = CheckVZEROchannel(v0side, partMC->Eta(), partMC->Phi());
797
798                         //use the efficiecy as a weigth
799                         multv0mc[iv0] += fV0Aeff->GetFunction("f")->Eval(partMC->Pt());
800
801                         //calculate multiplicity for each vzero channell
802                         //multv0mc[iv0]++;
803
804                 }// end loop on generated.
805
806                 Int_t firstCh=-1,lastCh=-1;
807                 if (v0side == 0) {
808                         firstCh = 32;
809                         lastCh = 64;
810                 }
811                 else if (v0side == 1) {
812                         firstCh = 0;
813                         lastCh = 32;
814                 }
815                 for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
816                         Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
817                         Double_t mult = multv0mc[iCh];
818                         Qx2mc += mult*TMath::Cos(2*phi);
819                         Qy2mc += mult*TMath::Sin(2*phi);
820                         mult2mc += mult;
821
822                         ((TH2F*)fOutput->FindObject("fV0Mmc"))->Fill(iCh,multv0mc[iCh]);
823                 }
824
825         }//end if on type==1
826
827         // return q vector
828         fQvecMC = TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
829
830         return fQvecMC;
831
832 }
833
834 //______________________________________________________
835 Int_t AliSpectraAODEventCuts::CheckVZEROchannel(Int_t vzeroside, Double_t eta, Double_t phi){
836
837         //VZEROA eta acceptance
838         Int_t ring = -1.;
839
840         if ( vzeroside == 0){
841                 if (eta > 4.5 && eta <= 5.1 ) ring = 0;
842                 if (eta > 3.9 && eta <= 4.5 ) ring = 1;
843                 if (eta > 3.4 && eta <= 3.9 ) ring = 2;
844                 if (eta > 2.8 && eta <= 3.4 ) ring = 3;
845         } else if (vzeroside == 1){
846                 if (eta > -3.7 && eta <= -3.2 ) ring = 0;
847                 if (eta > -3.2 && eta <= -2.7 ) ring = 1;
848                 if (eta > -2.7 && eta <= -2.2 ) ring = 2;
849                 if (eta > -2.2 && eta <= -1.7 ) ring = 3;
850         }
851
852
853         //VZEROAC phi acceptance
854         Int_t phisector= -1;
855
856
857         if ( phi > 0. && phi <= TMath::Pi()/4. ) phisector = 0;
858
859         else if (phi > TMath::Pi()/4. && phi <= TMath::Pi()/2.) phisector = 1;
860
861         else if (phi > TMath::Pi()/2 && phi <= (3./4.)*TMath::Pi() ) phisector =2;
862
863         else if (phi > (3./4.)*TMath::Pi() && phi <= TMath::Pi() ) phisector = 3;
864
865         else if (phi > TMath::Pi() && phi <= (5./4.)*TMath::Pi() ) phisector = 4;
866
867         else if (phi > (5./4.)*TMath::Pi() && phi <= (3./2.)*TMath::Pi() ) phisector = 5;
868
869         else if (phi > (3./2.)*TMath::Pi() && phi <= (7./4.)*TMath::Pi() ) phisector = 6;
870
871         else if (phi > (7./4.)*TMath::Pi() && phi <= TMath::TwoPi() ) phisector = 7;
872
873         if (vzeroside ==0) return Int_t(32+(phisector+(ring*8.))); // iv0 >= 32 -> V0A
874         if (vzeroside ==1) return Int_t(phisector+(ring*8.));  // iv0 < 32 -> V0C
875
876         else return -999.;
877 }
878
879 //______________________________________________________
880 Int_t AliSpectraAODEventCuts::CheckVZEROacceptance(Double_t eta){
881
882         // eval VZERO side - FIXME Add TPC!
883
884         if(eta > 2.8  && eta < 5.1) return 0; //VZEROA
885
886         else if (eta > -3.7  && eta < -1.7) return 1; //VZEROC
887
888         return -999.;
889
890 }
891
892 //______________________________________________________
893 Double_t AliSpectraAODEventCuts::GetQvecPercentileMC(Int_t v0side, Int_t type=1){
894
895         //check Qvec and Centrality consistency
896         if(fCent>90.) return -999.;
897
898         Double_t qvec = CalculateQVectorMC(v0side, type);
899         if(qvec==-999.) return -999.;
900
901         fQgenIntegral = 0x0;
902
903         if(type==0){
904                 if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen"); }
905                 if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen"); }
906         } else if (type==1){
907                 if(v0side==0/*V0A*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROAgen_vzero"); }
908                 if(v0side==1/*V0C*/){ fQgenIntegral = (TH2D*)fQvecIntList->FindObject("VZEROCgen_vzero"); }
909         }
910
911
912         //FIXME you need a check on the TH2D, add AliFatal or a return.
913
914         Int_t ic = -999;
915
916         if(fQvecCalibType==1){
917                 if(fNch<0.) return -999.;
918                 ic = GetNchBin(fQgenIntegral);
919         } else ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
920
921         TH1D *h1D = (TH1D*)fQgenIntegral->ProjectionY("h1Dgen",ic+1,ic+1);
922
923         TSpline *spline = 0x0;
924
925         if(v0side==0/*V0A*/){
926                 if( CheckSplineArray(fSplineArrayV0Agen, ic) ) {
927                         spline = (TSpline*)fSplineArrayV0Agen->At(ic);
928                         //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
929                 } else {
930                         spline = new TSpline3(h1D,"sp3");
931                         fSplineArrayV0Agen->AddAtAndExpand(spline,ic);
932                 }
933         }
934         else if(v0side==1/*V0C*/){
935                 if( CheckSplineArray(fSplineArrayV0Cgen, ic) ) {
936                         spline = (TSpline*)fSplineArrayV0Cgen->At(ic);
937                 } else {
938                         spline = new TSpline3(h1D,"sp3");
939                         fSplineArrayV0Cgen->AddAtAndExpand(spline,ic);
940                 }
941         }
942
943         Double_t percentile = 100*spline->Eval(qvec);
944
945         if(percentile>100. || percentile<0.) return -999.;
946
947         return percentile;
948
949 }
950
951 //______________________________________________________
952 Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr, Int_t n){
953         //check TSpline array consistency
954
955         //   Int_t n = (Int_t)fCent; //FIXME should be ok for icentr and nch
956
957         if(splarr->GetSize()<n) return kFALSE;
958
959         if(!splarr->At(n)) return kFALSE;
960
961         return kTRUE;
962
963 }
964
965 //______________________________________________________
966 Int_t AliSpectraAODEventCuts::GetNchBin(TH2D * h){
967
968         Double_t xmax = h->GetXaxis()->GetXmax();
969
970         if(fNch>xmax) return (Int_t)h->GetNbinsX();
971
972         return (fNch*h->GetNbinsX())/h->GetXaxis()->GetXmax();
973
974 }
975