Q vector selection in MC
[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   fqV0C(-999.),
67   fqV0A(-999.),
68   fqV0Cx(-999.),
69   fqV0Ax(-999.),
70   fqV0Cy(-999.),
71   fqV0Ay(-999.),
72   fPsiV0C(-999.),
73   fPsiV0A(-999.),
74   fCent(-999.),
75   fOutput(0),
76   fCalib(0),
77   fRun(-1),
78   fMultV0(0),
79   fV0Cpol1(-1),
80   fV0Cpol2(-1),
81   fV0Cpol3(-1),
82   fV0Cpol4(-1),
83   fV0Apol1(-1),
84   fV0Apol2(-1),
85   fV0Apol3(-1),
86   fV0Apol4(-1),
87   fQvecIntList(0),
88   fQvecIntegral(0), 
89   fSplineArrayV0A(0),
90   fSplineArrayV0C(0)
91 {
92   // Constructor
93   fOutput=new TList();
94   fOutput->SetOwner();
95   fOutput->SetName("fOutput");
96   
97   fCalib=new TList();
98   fCalib->SetOwner();
99   fCalib->SetName("fCalib");
100   
101   fQvecIntList=new TList();
102   fQvecIntList->SetOwner();
103   fQvecIntList->SetName("fQvecIntList");
104   
105   TH1I *fHistoCuts = new TH1I("fHistoCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
106   TH1F *fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection;z (cm)",500,-15,15);
107   TH1F *fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection;z (cm)",500,-15,15);
108   TH1F *fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection;eta",500,-2,2);
109   TH1F *fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection;eta",500,-2,2);
110   TH1F *fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection;Nch",2000,-0.5,1999.5);
111   TH2F *fHistoQVector = new TH2F("fHistoQVector", "QVector with VZERO distribution;centrality;Q vector from EP task",20,0,100,100,0,10);
112   TH2F *fHistoEP = new TH2F("fHistoEP", "EP with VZERO distribution;centrality;Psi_{EP} from EP task",20,0,100,100,-2,2);
113   TH2F *fPsiACor = new TH2F("fPsiACor", "EP with VZERO A distribution;centrality;Psi_{EP} VZERO-A",20,0,100,100,-2,2);
114   TH2F *fPsiCCor = new TH2F("fPsiCCor", "EP with VZERO C distribution;centrality;Psi_{EP} VZERO-C",20,0,100,100,-2,2);
115   TH2F *fQVecACor = new TH2F("fQVecACor", "QVec VZERO A;centrality;Qvector VZERO-A",20,0,100,100,0,10);
116   TH2F *fQVecCCor = new TH2F("fQVecCCor", "QVec VZERO C;centrality;Qvector VZERO-C",20,0,100,100,0,10);
117   TH2F *fV0M = new TH2F("fV0M", "V0 Multiplicity, before correction;V0 sector",64,-.5,63.5,500,0,1000);
118   TH2F *fV0MCor = new TH2F("fV0MCor", "V0 Multiplicity, after correction;V0 sector",64,-.5,63.5,500,0,1000);
119   
120   fSplineArrayV0A = new TObjArray();
121   fSplineArrayV0A->SetOwner();
122   fSplineArrayV0C = new TObjArray();
123   fSplineArrayV0C->SetOwner();
124   
125   fOutput->Add(fHistoCuts);
126   fOutput->Add(fHistoVtxBefSel);
127   fOutput->Add(fHistoVtxAftSel);
128   fOutput->Add(fHistoEtaBefSel);
129   fOutput->Add(fHistoEtaAftSel);
130   fOutput->Add(fHistoNChAftSel);
131   fOutput->Add(fHistoQVector);
132   fOutput->Add(fHistoEP);
133   fOutput->Add(fPsiACor);
134   fOutput->Add(fPsiCCor);
135   fOutput->Add(fQVecACor);
136   fOutput->Add(fQVecCCor);
137   fOutput->Add(fV0M);
138   fOutput->Add(fV0MCor);
139   
140   for (Int_t i = 0; i<10; i++){
141     fMeanQxa2[i] = -1;
142     fMeanQya2[i] = -1;
143     fMeanQxc2[i] = -1;
144     fMeanQyc2[i] = -1;   
145   }
146 }
147
148 //______________________________________________________
149 Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod,AliSpectraAODTrackCuts     *trackcuts)
150 {
151   // Returns true if Event Cuts are selected and applied
152   fAOD = aod;
153   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??
154   // FIXME: all those references by name are slow.
155   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kProcessedEvents);
156   Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fSelectBit);
157   if(!IsPhysSel)return IsPhysSel;
158   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kPhysSelEvents);
159   //loop on tracks, before event selection, filling QA histos
160   AliAODVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
161   if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxBefSel"))->Fill(vertex->GetZ());
162   fIsSelected =kFALSE;
163   if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut()){ //selection on vertex and Centrality
164     fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
165   }
166   if(fIsSelected){
167     ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kAcceptedEvents);
168     if(vertex)((TH1F*)fOutput->FindObject("fHistoVtxAftSel"))->Fill(vertex->GetZ());
169   }
170   Int_t Nch=0;
171   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
172     AliAODTrack* track = fAOD->GetTrack(iTracks);
173     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
174     ((TH1F*)fOutput->FindObject("fHistoEtaBefSel"))->Fill(track->Eta());
175     if(fIsSelected){
176       ((TH1F*)fOutput->FindObject("fHistoEtaAftSel"))->Fill(track->Eta());
177       Nch++;
178     }
179   }
180   if(fIsSelected)((TH1F*)fOutput->FindObject("fHistoNChAftSel"))->Fill(Nch);
181   return fIsSelected;
182 }
183
184 //______________________________________________________
185 Bool_t AliSpectraAODEventCuts::CheckVtxRange()
186 {
187   // reject events outside of range
188   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
189   //when moving to 2011 w├Čone has to add a cut using SPD vertex.
190   //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
191   if (!vertex)
192     {
193       ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxNoEvent);
194       return kFALSE;
195     }
196   if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
197     {
198       return kTRUE;
199     }
200   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxRange);
201   return kFALSE;
202 }
203
204 //______________________________________________________
205 Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
206 {
207   // Check centrality cut
208   fCent=-999.;
209   fCent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
210   if ( (fCent <= fCentralityCutMax)  &&  (fCent >= fCentralityCutMin) )  return kTRUE;   
211   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kVtxCentral);
212   return kFALSE;
213 }
214
215 //______________________________________________________
216 Bool_t AliSpectraAODEventCuts::CheckMultiplicityCut()
217 {
218   // Check multiplicity cut
219   // FIXME: why this is not tracket in the track stats histos?
220   Int_t Ncharged=0;
221   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
222     AliAODTrack* track = fAOD->GetTrack(iTracks);
223     if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
224     Ncharged++;
225   }
226   if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
227   
228   return kFALSE;
229 }
230
231 //______________________________________________________
232 Bool_t AliSpectraAODEventCuts::CheckQVectorCut()
233
234   Double_t qxEPVZERO = -999., qyEPVZERO = -999.;
235   Double_t qVZERO = -999.;
236   Double_t psi=-999.;
237   
238   if(fIsLHC10h){
239     qVZERO=CalculateQVectorLHC10h();
240     psi=fPsiV0A;
241   }else{
242     psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);//FIXME we can a flag for 2010 and 2011
243     qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
244   }
245   
246   //cut on q vector
247   if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
248   Double_t cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentralityMethod.Data());
249   ((TH2F*)fOutput->FindObject("fHistoQVector"))->Fill(cent,qVZERO);
250   ((TH2F*)fOutput->FindObject("fHistoEP"))->Fill(cent,psi);
251   ((TH1I*)fOutput->FindObject("fHistoCuts"))->Fill(kQVector);
252   
253
254   return kTRUE;
255 }
256
257 //______________________________________________________
258 Double_t AliSpectraAODEventCuts::CalculateQVectorLHC10h(){
259   
260   Int_t run = fAOD->GetRunNumber();
261   if(run != fRun){
262     // Load the calibrations run dependent
263     if(OpenInfoCalbration(run))fRun=run;
264     else{
265       fqV0C=-999.;
266       fqV0A=-999.;
267       return -999.;
268     }
269   }
270   
271   //V0 info    
272   Double_t Qxa2 = 0, Qya2 = 0;
273   Double_t Qxc2 = 0, Qyc2 = 0;
274   Double_t sumMc = 0, sumMa = 0;
275   
276   AliAODVZERO* aodV0 = fAOD->GetVZEROData();
277   
278   for (Int_t iv0 = 0; iv0 < 64; iv0++) {
279     
280     Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
281     
282     Float_t multv0 = aodV0->GetMultiplicity(iv0);
283     ((TH2F*)fOutput->FindObject("fV0M"))->Fill(iv0,multv0);
284     
285     if (iv0 < 32){
286       
287       Double_t multCorC = -10;
288       
289       if (iv0 < 8)
290         multCorC = multv0*fV0Cpol1/fMultV0->GetBinContent(iv0+1);
291       else if (iv0 >= 8 && iv0 < 16)
292         multCorC = multv0*fV0Cpol2/fMultV0->GetBinContent(iv0+1);
293       else if (iv0 >= 16 && iv0 < 24)
294         multCorC = multv0*fV0Cpol3/fMultV0->GetBinContent(iv0+1);
295       else if (iv0 >= 24 && iv0 < 32)
296         multCorC = multv0*fV0Cpol4/fMultV0->GetBinContent(iv0+1);
297       
298       if (multCorC < 0){
299         cout<<"Problem with multiplicity in V0C"<<endl;
300         fqV0C=-999.;
301         fqV0A=-999.;
302         return -999.;
303       }
304       
305       Qxc2 += TMath::Cos(2*phiV0) * multCorC;
306       Qyc2 += TMath::Sin(2*phiV0) * multCorC;
307       
308       sumMc += multCorC;
309       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorC);
310       
311     } else {
312       
313       Double_t multCorA = -10;
314       
315       if (iv0 >= 32 && iv0 < 40)
316         multCorA = multv0*fV0Apol1/fMultV0->GetBinContent(iv0+1);
317       else if (iv0 >= 40 && iv0 < 48)
318         multCorA = multv0*fV0Apol2/fMultV0->GetBinContent(iv0+1);
319       else if (iv0 >= 48 && iv0 < 56)
320         multCorA = multv0*fV0Apol3/fMultV0->GetBinContent(iv0+1);
321       else if (iv0 >= 56 && iv0 < 64)
322         multCorA = multv0*fV0Apol4/fMultV0->GetBinContent(iv0+1);
323       
324       if (multCorA < 0){
325         cout<<"Problem with multiplicity in V0A"<<endl;
326         fqV0C=-999.;
327         fqV0A=-999.;
328         return -999.;
329       }
330       
331       Qxa2 += TMath::Cos(2*phiV0) * multCorA;
332       Qya2 += TMath::Sin(2*phiV0) * multCorA;
333       
334       sumMa += multCorA;
335       ((TH2F*)fOutput->FindObject("fV0MCor"))->Fill(iv0,multCorA);
336     }
337   }
338   
339   Short_t centrV0  = GetCentrCode(fAOD);
340   
341   Double_t Qxamean2 = 0.;
342   Double_t Qyamean2 = 0.;
343   Double_t Qxcmean2 = 0.;
344   Double_t Qycmean2 = 0.;
345   
346   if(centrV0!=-1){
347     Qxamean2 = fMeanQxa2[centrV0];
348     Qyamean2 = fMeanQya2[centrV0];
349     Qxcmean2 = fMeanQxc2[centrV0];
350     Qycmean2 = fMeanQyc2[centrV0];
351   }
352   
353   Double_t QxaCor2 = Qxa2 - Qxamean2*sumMa;
354   Double_t QyaCor2 = Qya2 - Qyamean2*sumMa;
355   Double_t QxcCor2 = Qxc2 - Qxcmean2*sumMc;
356   Double_t QycCor2 = Qyc2 - Qycmean2*sumMc;
357   
358   fPsiV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
359   fPsiV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
360   
361   ((TH2F*)fOutput->FindObject("fPsiACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0A);
362   ((TH2F*)fOutput->FindObject("fPsiCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fPsiV0C);
363   
364   fqV0A = TMath::Sqrt((QxaCor2*QxaCor2 + QyaCor2*QyaCor2)/sumMa);
365   fqV0C = TMath::Sqrt((QxcCor2*QxcCor2 + QycCor2*QycCor2)/sumMc);
366   fqV0Ax = QxaCor2*TMath::Sqrt(1./sumMa);
367   fqV0Cx = QxcCor2*TMath::Sqrt(1./sumMc);
368   fqV0Ay = QyaCor2*TMath::Sqrt(1./sumMa);
369   fqV0Cy = QycCor2*TMath::Sqrt(1./sumMc);
370   
371   ((TH2F*)fOutput->FindObject("fQVecACor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0A);
372   ((TH2F*)fOutput->FindObject("fQVecCCor"))->Fill((Float_t)fAOD->GetCentrality()->GetCentralityPercentile("V0M"), fqV0C);
373   
374   return fqV0A; //FIXME we have to combine VZERO-A and C
375 }
376
377 //______________________________________________________
378 Short_t AliSpectraAODEventCuts::GetCentrCode(AliVEvent* ev)
379 {
380   
381   Short_t centrCode = -1;
382   
383   AliCentrality* centrality = 0;
384   AliAODEvent* aod = (AliAODEvent*)ev;  
385   centrality = aod->GetHeader()->GetCentralityP();
386   
387   Float_t centV0 = centrality->GetCentralityPercentile("V0M");
388   Float_t centTrk = centrality->GetCentralityPercentile("TRK"); 
389   
390   
391   if (TMath::Abs(centV0 - centTrk) < 5.0 && centV0 <= 80 && centV0 > 0){  
392     
393      if ((centV0 > 0) && (centV0 <= 5.0))
394        centrCode = 0; 
395      else if ((centV0 > 5.0) && (centV0 <= 10.0))
396        centrCode = 1;
397      else if ((centV0 > 10.0) && (centV0 <= 20.0))
398        centrCode = 2;
399      else if ((centV0 > 20.0) && (centV0 <= 30.0))
400        centrCode = 3;   
401      else if ((centV0 > 30.0) && (centV0 <= 40.0))
402        centrCode = 4; 
403      else if ((centV0 > 40.0) && (centV0 <= 50.0))
404        centrCode = 5;  
405      else if ((centV0 > 50.0) && (centV0 <= 60.0))
406        centrCode = 6;
407      else if ((centV0 > 60.0) && (centV0 <= 70.0))
408        centrCode = 7;                                  
409      else if ((centV0 > 70.0) && (centV0 <= 80.0))
410        centrCode = 8;  
411   }
412   
413   return centrCode;
414   
415 }
416
417 //______________________________________________________
418 void AliSpectraAODEventCuts::PrintCuts()
419 {
420   // print info about event cuts
421   cout << "Event Stats" << endl;
422   cout << " > Trigger Selection: " << fSelectBit << endl;
423   cout << " > Centrality estimator: " << fCentralityMethod << endl;
424   cout << " > Number of accepted events: " << NumberOfEvents() << endl;
425   cout << " > Number of processed events: " << NumberOfProcessedEvents() << endl;
426   cout << " > Number of PhysSel events: " <<  NumberOfPhysSelEvents()  << endl;
427   cout << " > Vertex out of range: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxRange + 1) << endl;
428   cout << " > Events cut by centrality: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxCentral + 1) << endl;
429   cout << " > Events without vertex: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kVtxNoEvent + 1) << endl;
430   cout << " > QVector cut: " << ((TH1I*)fOutput->FindObject("fHistoCuts"))->GetBinContent(kQVector + 1) << endl;
431   cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
432   cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
433   cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
434   cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
435   cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
436 }
437
438 //_____________________________________________________________________________
439 Bool_t AliSpectraAODEventCuts::OpenInfoCalbration(Int_t run)
440 {
441   
442   AliOADBContainer* cont = (AliOADBContainer*) fCalib->FindObject("hMultV0BefCorr");
443   if(!cont){
444     printf("OADB object hMultV0BefCorr is not available in the file\n");
445     return 0;   
446   }
447   
448   if(!(cont->GetObject(run))){
449     printf("OADB object hMultV0BefCorr is not available for run %i\n",run);
450     return 0;   
451   }
452   fMultV0 = ((TH2F*) cont->GetObject(run))->ProfileX();
453   
454   TF1* fpolc1 = new TF1("fpolc1","pol0", 0, 7); 
455   fMultV0->Fit(fpolc1, "RN");
456   fV0Cpol1 = fpolc1->GetParameter(0);
457   
458   TF1* fpolc2 = new TF1("fpolc2","pol0", 8, 15); 
459   fMultV0->Fit(fpolc2, "RN");
460   fV0Cpol2 = fpolc2->GetParameter(0);
461   
462   TF1* fpolc3 = new TF1("fpolc3","pol0", 16, 23); 
463   fMultV0->Fit(fpolc3, "RN");
464   fV0Cpol3 = fpolc3->GetParameter(0);
465   
466   TF1* fpolc4 = new TF1("fpolc4","pol0", 24, 31); 
467   fMultV0->Fit(fpolc4, "RN");
468   fV0Cpol4 = fpolc4->GetParameter(0);
469   
470   TF1* fpola1 = new TF1("fpola1","pol0", 32, 39);
471   fMultV0->Fit(fpola1, "RN");
472   fV0Apol1 = fpola1->GetParameter(0);
473
474   TF1* fpola2 = new TF1("fpola2","pol0", 40, 47);
475   fMultV0->Fit(fpola2, "RN");
476   fV0Apol2 = fpola2->GetParameter(0);
477   
478   TF1* fpola3 = new TF1("fpola3","pol0", 48, 55);
479   fMultV0->Fit(fpola3, "RN");
480   fV0Apol3 = fpola3->GetParameter(0);
481   
482   TF1* fpola4 = new TF1("fpola4","pol0", 56, 63);
483   fMultV0->Fit(fpola4, "RN");
484   fV0Apol4 = fpola4->GetParameter(0);
485
486   for(Int_t i=0; i < 10; i++){
487
488     char nameQxa2[100];
489     snprintf(nameQxa2,100, "hQxa2m_%i", i);
490
491     char nameQya2[100];
492     snprintf(nameQya2,100, "hQya2m_%i", i);
493
494     char nameQxc2[100];
495     snprintf(nameQxc2,100, "hQxc2m_%i", i);
496
497     char nameQyc2[100];
498     snprintf(nameQyc2,100, "hQyc2m_%i", i);
499         
500     AliOADBContainer* contQxa2 = (AliOADBContainer*) fCalib->FindObject(nameQxa2);
501     if(!contQxa2){
502       printf("OADB object %s is not available in the file\n", nameQxa2);
503       return 0; 
504     }
505         
506     if(!(contQxa2->GetObject(run))){
507       printf("OADB object %s is not available for run %i\n", nameQxa2, run);
508       return 0; 
509     }
510
511     fMeanQxa2[i] = ((TH1F*) contQxa2->GetObject(run))->GetMean();
512
513  
514     AliOADBContainer* contQya2 = (AliOADBContainer*) fCalib->FindObject(nameQya2);
515     if(!contQya2){
516       printf("OADB object %s is not available in the file\n", nameQya2);
517       return 0; 
518     }
519         
520     if(!(contQya2->GetObject(run))){
521       printf("OADB object %s is not available for run %i\n", nameQya2, run);
522       return 0; 
523     }
524
525     fMeanQya2[i] = ((TH1F*) contQya2->GetObject(run))->GetMean();
526
527
528     AliOADBContainer* contQxc2 = (AliOADBContainer*) fCalib->FindObject(nameQxc2);
529     if(!contQxc2){
530       printf("OADB object %s is not available in the file\n", nameQxc2);
531       return 0; 
532     }
533         
534     if(!(contQxc2->GetObject(run))){
535       printf("OADB object %s is not available for run %i\n", nameQxc2, run);
536       return 0; 
537     }
538
539     fMeanQxc2[i] = ((TH1F*) contQxc2->GetObject(run))->GetMean();
540  
541
542     AliOADBContainer* contQyc2 = (AliOADBContainer*) fCalib->FindObject(nameQyc2);
543     if(!contQyc2){
544       printf("OADB object %s is not available in the file\n", nameQyc2);
545       return 0; 
546     }
547         
548     if(!(contQyc2->GetObject(run))){
549       printf("OADB object %s is not available for run %i\n", nameQyc2, run);
550       return 0; 
551     }
552
553     fMeanQyc2[i] = ((TH1F*) contQyc2->GetObject(run))->GetMean();
554
555   }
556   return 1;
557 }
558
559 //______________________________________________________
560
561
562 Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
563 {
564   // Merge a list of AliSpectraAODEventCuts objects with this.
565   // Returns the number of merged objects (including this).
566   
567   AliInfo("Merging");
568   
569   if (!list)
570     return 0;
571
572   if (list->IsEmpty())
573     return 1;
574   
575   TIterator* iter = list->MakeIterator();
576   TObject* obj;
577   
578   // collections of all histograms
579   TList collections;
580   
581   Int_t count = 0;
582
583   while ((obj = iter->Next())) {
584     AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
585     if (entry == 0) 
586       continue;
587
588     TList * l = entry->GetOutputList();      
589     collections.Add(l);
590     count++;
591   }
592   
593   fOutput->Merge(&collections);
594   
595   delete iter;
596
597   return count+1;
598 }
599
600 //______________________________________________________
601 Double_t AliSpectraAODEventCuts::GetQvecPercentile(Int_t v0side){
602  
603   //check Qvec and Centrality consistency
604   if(fCent>90.) return -999.;
605   if(v0side==0/*V0A*/){ if(fqV0A== -999.) return -999.; }
606   if(v0side==1/*V0C*/){ if(fqV0C== -999.) return -999.; }
607   
608   fQvecIntegral = 0x0;
609
610   if(v0side==0/*V0A*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROA"); }
611   if(v0side==1/*V0C*/){ fQvecIntegral = (TH2D*)fQvecIntList->FindObject("VZEROC"); }
612
613   Double_t ic = (Int_t)fCent; //fQvecIntegral: 1% centrality bin
614   
615   TH1D *h1D = (TH1D*)fQvecIntegral->ProjectionY("h1D",ic+1,ic+1);
616   
617   TSpline *spline = 0x0;
618   
619   if(v0side==0/*V0A*/){
620     if( CheckSplineArray(fSplineArrayV0A) ) {
621       spline = (TSpline*)fSplineArrayV0A->At(ic);
622       //cout<<"Qvec V0A - ic: "<<ic<<" - Found TSpline..."<<endl;
623     } else {
624       spline = new TSpline3(h1D,"sp3");
625       fSplineArrayV0A->AddAtAndExpand(spline,ic);
626       //cout<<"Qvec V0A - ic: "<<ic<<" - TSpline not found. Creating..."<<endl;
627     }
628   }
629   else if(v0side==1/*V0C*/){
630     if( CheckSplineArray(fSplineArrayV0C) ) {
631       spline = (TSpline*)fSplineArrayV0C->At(ic);
632     } else {
633       spline = new TSpline3(h1D,"sp3");
634       fSplineArrayV0C->AddAtAndExpand(spline,ic);
635     }
636   }
637   
638   Double_t percentile =  -999.;
639   if(v0side==0/*V0A*/){ percentile = 100*spline->Eval(fqV0A); }
640   if(v0side==1/*V0C*/){ percentile = 100*spline->Eval(fqV0C); }
641   
642   if(percentile>100. || percentile<0.) return -999.;
643
644   return percentile;
645 }
646
647 //______________________________________________________
648 Double_t AliSpectraAODEventCuts::CalculateQVectorMC(Int_t v0side){
649   
650   if(!fIsMC) return -999.;
651   
652   // 1. get MC array
653   TClonesArray *arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
654   
655   if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
656   
657   Double_t Qx2mc = 0., Qy2mc = 0.;
658   Int_t mult2mc = 0;
659   
660   Int_t nMC = arrayMC->GetEntries();
661       
662   // 2. loop on generated
663   for (Int_t iMC = 0; iMC < nMC; iMC++){
664     AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
665   
666     // 3. set VZERO side - FIXME Add TPC!
667     Double_t EtaVZERO[2] = {-999.,-999.};
668     if(v0side==0/*V0A*/){ EtaVZERO[0] = 2.8; EtaVZERO[1] = 5.1; } 
669     if(v0side==1/*V0C*/){ EtaVZERO[0] = -3.7; EtaVZERO[1] = -1.7; } 
670     
671     if(partMC->Eta()<EtaVZERO[0] || partMC->Eta()>EtaVZERO[1]) continue;
672     
673     // 4. Calculate Qvec components
674     
675     Qx2mc += TMath::Cos(2.*partMC->Phi());
676     Qy2mc += TMath::Sin(2.*partMC->Phi());
677     mult2mc++;
678     
679   }
680   
681   // 5. return q vector
682   return TMath::Sqrt((Qx2mc*Qx2mc + Qy2mc*Qy2mc)/mult2mc);
683   
684 }
685
686 //______________________________________________________
687 Bool_t AliSpectraAODEventCuts::CheckSplineArray(TObjArray * splarr){
688   //check TSpline array consistency
689   
690   Int_t n = (Int_t)fCent;
691   
692   if(splarr->GetSize()<n) return kFALSE;
693   
694   if(!splarr->At(n)) return kFALSE;
695     
696   return kTRUE;
697   
698 }