adding AliPPVsMultUtils to event cuts
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / AliSpectraBothEventCuts.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 //         AliSpectraBothEventCuts 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 "AliAnalysisTask.h"
29 #include "AliAnalysisManager.h"
30 #include "AliAODTrack.h"
31 #include "AliAODMCParticle.h"
32 #include "AliAODEvent.h"
33 #include "AliESDEvent.h"
34 #include "AliAODInputHandler.h"
35 #include "AliAnalysisTaskESDfilter.h"
36 #include "AliAnalysisDataContainer.h"
37 #include "AliSpectraBothEventCuts.h"
38 #include "AliSpectraBothTrackCuts.h"
39 #include "AliAnalysisUtils.h"
40 #include "AliPPVsMultUtils.h"   
41 //#include "AliSpectraBothHistoManager.h"
42 #include <iostream>
43
44 using namespace std;
45
46 ClassImp(AliSpectraBothEventCuts)
47
48 AliSpectraBothEventCuts::AliSpectraBothEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0),fAODEvent(AliSpectraBothTrackCuts::kAODobject), fTrackBits(0),fIsMC(0),fCentEstimator(""), fUseCentPatchAOD049(0), fUseSDDPatchforLHC11a(kDoNotCheckforSDD),fTriggerSettings(AliVEvent::kMB),fTrackCuts(0),
49 fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(0),fMaxChi2perNDFforVertex(0),
50 fMinRun(0),fMaxRun(0),fetarangeofmultiplicitycut(0.0),fUseAliPPVsMultUtils(false),
51 fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
52 ,fHistoEP(0),fHistoVtxAftSelwithoutZvertexCut(0),fHistoVtxalltriggerEventswithMCz(0),fHistoVtxAftSelwithoutZvertexCutusingMCz(0),fHistoRunNumbers(0),
53 fHistoCentrality(0),fHistoMultiplicty(0),fAnalysisUtils(0),fAliPPVsMultUtils(0)
54
55 {
56   // Constructori
57  // Bool_t oldStatus = TH1::AddDirectoryStatus();
58   //TH1::AddDirectory(kFALSE);  
59  // fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
60  // fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",300,-15,15);
61 //  fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",300,-15,15);
62 //  fHistoVtxAftSelwithoutZvertexCut=new TH1F("fHistoVtxAftSelwithoutZvertexcut", "Vtx distr after event selection without Z vertex cut",300,-15,15);
63 //  fHistoVtxalltriggerEventswithMCz=new TH1F("fHistoVtxalltriggerEventswithMCz", "generated z vertex position",300,-15,15);
64 //  fHistoVtxAftSelwithoutZvertexCutusingMCz=new TH1F("fHistoVtxAftSelwithoutZvertexCutusingMCz", "Vtx distr after event selection without Z vertex cut using MC z",300,-15,15);
65 //  fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
66 //  fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
67  // fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
68   //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
69   //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
70  // fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
71  // fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
72   fCentralityCutMin = 0.0;      // default value of centrality cut minimum, 0 ~ no cut
73   fCentralityCutMax = 10000.0;  // default value of centrality cut maximum,  ~ no cut
74   // fQVectorPosCutMin=0.0;
75   // fQVectorPosCutMax=10000.0;
76   // fQVectorNegCutMin=0.0;
77   // fQVectorNegCutMax=10000.0;
78   fQVectorCutMin=0.0;
79   fQVectorCutMax=10000.0;
80   fVertexCutMin=-10.0;
81   fVertexCutMax=10.0;
82   fMultiplicityCutMin=-1.0;
83   fMultiplicityCutMax=-1.0;
84   fTrackBits=1;
85   fCentEstimator="V0M";
86   fMaxChi2perNDFforVertex=-1;
87  // TH1::AddDirectory(oldStatus);       
88 }
89 //______________________________________________________
90
91 AliSpectraBothEventCuts::~AliSpectraBothEventCuts()
92 {
93         if(fHistoCuts)
94                 delete fHistoCuts;
95         if(fHistoVtxBefSel)
96                 delete fHistoVtxBefSel;
97         if(fHistoVtxAftSel)
98                 delete fHistoVtxAftSel;
99         if(fHistoVtxAftSelwithoutZvertexCut)
100                 delete fHistoVtxAftSelwithoutZvertexCut;
101         if(fHistoVtxalltriggerEventswithMCz)
102                 delete fHistoVtxalltriggerEventswithMCz;
103         if(fHistoVtxAftSelwithoutZvertexCutusingMCz)
104                 delete fHistoVtxAftSelwithoutZvertexCutusingMCz;
105         if(fHistoEtaBefSel)
106                 delete fHistoEtaBefSel;
107         if(fHistoEtaAftSel)
108                 delete fHistoEtaAftSel ;
109         if(fHistoNChAftSel)
110                 delete fHistoNChAftSel;
111         if(fHistoQVector)
112                 delete fHistoQVector;
113         if(fHistoEP)
114                 delete fHistoEP;
115         if(fHistoRunNumbers)
116                 delete fHistoRunNumbers;
117         if(fHistoCentrality)
118                 delete fHistoCentrality;
119         if(fHistoMultiplicty)
120                 delete fHistoMultiplicty;
121         if(fAnalysisUtils)
122                 delete fAnalysisUtils;
123         if(fAliPPVsMultUtils)
124                 delete fAliPPVsMultUtils;
125
126 }
127 //______________________________________________________
128 void AliSpectraBothEventCuts::InitHisto()
129 {
130         Bool_t oldStatus = TH1::AddDirectoryStatus();
131         TH1::AddDirectory(kFALSE);
132         if(!fHistoCuts) 
133                 fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
134         if(!fHistoVtxBefSel )
135                 fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",300,-15,15);
136         if(!fHistoVtxAftSel)
137                 fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",300,-15,15);
138         if(!fHistoVtxAftSelwithoutZvertexCut)
139                 fHistoVtxAftSelwithoutZvertexCut=new TH1F("fHistoVtxAftSelwithoutZvertexcut", "Vtx distr after event selection without Z vertex cut",300,-15,15);
140         if(!fHistoVtxalltriggerEventswithMCz)
141                 fHistoVtxalltriggerEventswithMCz=new TH1F("fHistoVtxalltriggerEventswithMCz", "generated z vertex position",300,-15,15);
142         if(!fHistoVtxAftSelwithoutZvertexCutusingMCz)
143                 fHistoVtxAftSelwithoutZvertexCutusingMCz=new TH1F("fHistoVtxAftSelwithoutZvertexCutusingMCz", "Vtx distr after event selection without Z vertex cut using MC z",300,-15,15);
144         if(!fHistoEtaBefSel)
145                 fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
146         if(!fHistoEtaAftSel)
147                 fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
148         if(!fHistoNChAftSel)
149                 fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
150   //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
151   //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
152         if(!fHistoQVector)
153                 fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
154         if(!fHistoEP)
155                 fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
156         if(!fHistoRunNumbers)
157         {
158                 if(fMaxRun>fMinRun&&fMinRun>=0)
159                         fHistoRunNumbers=new TH1F("fHistoRunNumbers","Run numbers",fMaxRun-fMinRun+1,fMinRun-0.5,fMaxRun+0.5);
160                 else
161                         fHistoRunNumbers=new TH1F("fHistoRunNumbers","Run numbers",1001,120000-.5,121000+0.5);
162
163         }
164         if(!fHistoCentrality)
165                 fHistoCentrality = new TH2F("fHistoCentrality", "centrality",2,0,2,100,0.0,100);
166
167         if(!fHistoMultiplicty)
168                 fHistoMultiplicty= new TH2F("fHistoMultiplicty", "multiplicty estimator",2,0,2,155,-4.5,150.5);
169
170         TH1::AddDirectory(oldStatus);           
171 }
172 //______________________________________________________
173 Bool_t AliSpectraBothEventCuts::IsSelected(AliVEvent * aod,AliSpectraBothTrackCuts* trackcuts,Bool_t isMC,Double_t mcZ)
174 {
175   // Returns true if Event Cuts are selected and applied
176   fAOD = aod;
177   fTrackCuts = trackcuts;
178   fHistoCuts->Fill(kProcessedEvents);
179   fHistoRunNumbers->Fill(aod->GetRunNumber());
180   Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerSettings);//FIXME we can add the trigger mask here
181   if(!IsPhysSel)
182         return IsPhysSel;
183   if(fAnalysisUtils) // we check for pile-up
184                  IsPhysSel = (!fAnalysisUtils->IsPileUpEvent(fAOD));
185   if(!IsPhysSel)
186         return IsPhysSel;
187         
188   if(isMC)      
189         fHistoVtxalltriggerEventswithMCz->Fill(mcZ);
190    //loop on tracks, before event selection, filling QA histos
191   AliESDEvent* esdevent=0x0;
192   AliAODEvent* aodevent=0x0;
193   Bool_t isSDD=kFALSE;
194   TString nameoftrack(fAOD->ClassName());  
195     if(!nameoftrack.CompareTo("AliESDEvent"))
196     {
197                 fAODEvent=AliSpectraBothTrackCuts::kESDobject;
198                 
199                 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
200                 {
201                         esdevent=dynamic_cast<AliESDEvent*>(fAOD);
202                         if(!esdevent)
203                                 return kFALSE;
204                         if(esdevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
205                                 isSDD=kTRUE;
206                 }       
207         }
208         else if(!nameoftrack.CompareTo("AliAODEvent"))
209         {
210                 fAODEvent=AliSpectraBothTrackCuts::kAODobject;
211                 if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
212                 {
213                         aodevent=dynamic_cast<AliAODEvent*>(fAOD);
214                         if(!aodevent)
215                                 return kFALSE;
216                         if(aodevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
217                                 isSDD=kTRUE;    
218                 }       
219         }
220         else
221                 return false; 
222       if(fUseSDDPatchforLHC11a==kwithSDD&&isSDD==kFALSE)
223                 return false;
224      if(fUseSDDPatchforLHC11a==kwithoutSDD&&isSDD==kTRUE)
225                 return false;
226
227     fHistoCuts->Fill(kPhysSelEvents);
228
229
230
231
232    const AliVVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated     
233
234   if(vertex)
235         fHistoVtxBefSel->Fill(vertex->GetZ());
236   fIsSelected =kFALSE;
237   if(CheckVtx())
238    { //selection on vertex
239
240     fIsSelected=kTRUE;
241   }
242   if(fIsSelected&&vertex)
243   {
244         fHistoVtxAftSelwithoutZvertexCut->Fill(vertex->GetZ());
245         if(isMC)
246                 fHistoVtxAftSelwithoutZvertexCutusingMCz->Fill(mcZ);    
247         if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
248         {
249                 fHistoCuts->Fill(kAcceptedEvents);
250                 fIsSelected=kTRUE;
251                 fHistoVtxAftSel->Fill(vertex->GetZ());
252         }
253         else    
254         {
255                 fIsSelected=kFALSE;
256         }       
257   }
258   if(fIsSelected)
259   {
260         if( CheckCentralityCut() && CheckMultiplicityCut() && CheckQVectorCut())
261                 fIsSelected=kTRUE;
262   }     
263
264
265
266   Int_t Nch=0;
267   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) 
268   {
269          AliVTrack* track =dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
270    /* if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
271                 track=dynamic_cast<AliVTrack*>(esdevent->GetTrack(iTracks));
272      else if (fAODEvent==AliSpectraBothTrackCuts::kAODobject)
273                 track=dynamic_cast<AliVTrack*>(aodevent->GetTrack(iTracks));
274      else return false;*/
275      
276          if (!fTrackCuts->IsSelected(track,kFALSE)) 
277                 continue;
278         fHistoEtaBefSel->Fill(track->Eta());
279          if(fIsSelected)
280         {
281                 fHistoEtaAftSel->Fill(track->Eta());
282                 Nch++;
283         }
284   }
285   if(fIsSelected)
286         fHistoNChAftSel->Fill(Nch);
287   return fIsSelected;
288 }
289
290 //______________________________________________________
291 Bool_t AliSpectraBothEventCuts::CheckVtx()
292 {
293   // reject events outside of range
294   const AliVVertex * vertex = fAOD->GetPrimaryVertex();
295   //when moving to 2011 w├Čone has to add a cut using SPD vertex.
296   //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
297   if (!vertex)
298     {
299       fHistoCuts->Fill(kVtxNoEvent);
300       return kFALSE;
301     }
302     if(vertex->GetNContributors()<1)
303    {
304
305       fHistoCuts->Fill(kZeroCont);
306       return kFALSE;
307
308    }
309         
310    TString tmp(vertex->GetTitle());
311    if(tmp.Contains("NoConstraint"))
312    {
313         fHistoCuts->Fill(kTPCasPV);
314         return kFALSE;
315    }
316
317
318  // if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
319    // {
320     //  return kTRUE;
321    // }
322    if(!CheckVtxChi2perNDF())      
323         return kFALSE;
324                 
325  fHistoCuts->Fill(kGoodVtx);
326   //return kFALSE;
327    return kTRUE;
328 }
329
330 //______________________________________________________
331 Bool_t AliSpectraBothEventCuts::CheckCentralityCut()
332 {
333   // Check centrality cut
334  if ( fCentralityCutMax<0.0  &&  fCentralityCutMin<0.0 )  
335         return kTRUE;
336   Double_t cent=0;      
337   if(fUseAliPPVsMultUtils)
338   {
339         if(!fAliPPVsMultUtils)
340                 fAliPPVsMultUtils=new AliPPVsMultUtils();
341         cent=fAliPPVsMultUtils->GetMultiplicityPercentile(fAOD,fCentEstimator.Data());
342   }     
343   else
344   {
345         if(!fUseCentPatchAOD049)
346                 cent=fAOD->GetCentrality()->GetCentralityPercentile(fCentEstimator.Data());
347         else 
348                 cent=ApplyCentralityPatchAOD049();
349   }     
350   fHistoCentrality->Fill(0.5,cent);     
351   if ( (cent <= fCentralityCutMax)  &&  (cent >= fCentralityCutMin) )  
352   {
353          fHistoCentrality->Fill(1.5,cent);      
354         return kTRUE;
355   }        
356   fHistoCuts->Fill(kVtxCentral);
357
358   return kFALSE;
359 }
360
361 //______________________________________________________
362 Bool_t AliSpectraBothEventCuts::CheckMultiplicityCut()
363 {
364   // Check multiplicity cut
365         if(fMultiplicityCutMin<0 && fMultiplicityCutMax<0)
366                 return kTRUE;
367         Int_t Ncharged=-1;
368         if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
369         {       
370                   AliESDEvent* esdevent=dynamic_cast<AliESDEvent*>(fAOD);
371                  AliESDtrackCuts::MultEstTrackType estType = esdevent->GetPrimaryVertexTracks()->GetStatus() ? AliESDtrackCuts::kTrackletsITSTPC : AliESDtrackCuts::kTracklets;
372                 Ncharged=AliESDtrackCuts::GetReferenceMultiplicity(esdevent,estType,fetarangeofmultiplicitycut);
373         }
374         else if(fAODEvent==AliSpectraBothTrackCuts::kAODobject)
375         {
376                 AliAODEvent* aodevent=0x0;
377                 aodevent=dynamic_cast<AliAODEvent*>(fAOD);
378                 if(TMath::Abs(0.8-fetarangeofmultiplicitycut)<0.1)
379                         Ncharged=aodevent->GetHeader()->GetRefMultiplicityComb08();
380                 else if (TMath::Abs(0.5-fetarangeofmultiplicitycut)<0.1)
381                         Ncharged=aodevent->GetHeader()->GetRefMultiplicityComb05();
382                 else 
383                         Ncharged=-1;
384         }
385         else
386                 return kFALSE;   
387
388    fHistoMultiplicty->Fill(0.5,Ncharged);
389    if(Ncharged>=fMultiplicityCutMin && Ncharged<fMultiplicityCutMax&& Ncharged>0)
390    { 
391         fHistoMultiplicty->Fill(1.5,Ncharged);
392         return kTRUE;
393    }
394   fHistoCuts->Fill(kVtxCentral);
395   return kFALSE;
396 }
397
398 //______________________________________________________
399 Bool_t AliSpectraBothEventCuts::CheckQVectorCut()
400 {
401          if(fQVectorCutMin<0.0 && fQVectorCutMax<0.0)return kTRUE;
402   // Check qvector cut
403   /// FIXME: Q vector
404   // //Selection on QVector, before ANY other selection on the event
405   // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
406   // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
407   // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
408  
409   // Int_t multPos = 0;
410   // Int_t multNeg = 0;
411   // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
412   //   AliAODTrack* aodTrack = fAOD->GetTrack(iT);
413   //   if (!aodTrack->TestFilterBit(fTrackBits)) continue;
414   //   if (aodTrack->Eta() >= 0){
415   //     multPos++;
416   //     Qx2EtaPos += TMath::Cos(2*aodTrack->Phi()); 
417   //     Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
418   //   } else {
419   //     multNeg++;
420   //     Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi()); 
421   //     Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
422   //   }
423   // } 
424   // Double_t qPos=-999;
425   // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
426   // Double_t qNeg=-999;
427   // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
428   //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
429   
430   Double_t qxEPVZERO = 0, qyEPVZERO = 0;
431   Double_t qVZERO = -999;
432   Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
433   
434   qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
435   if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
436   fHistoQVector->Fill(qVZERO);
437   fHistoEP->Fill(psi);
438   
439   fHistoCuts->Fill(kQVector);
440   // fHistoQVectorPos->Fill(qPos);
441   // fHistoQVectorNeg->Fill(qNeg);
442   return kTRUE;
443 }
444 //____________________________________________________________
445  Bool_t AliSpectraBothEventCuts::CheckVtxChi2perNDF()
446  {
447         if(fMaxChi2perNDFforVertex<0)
448                 return kTRUE;
449          const AliVVertex * vertex = fAOD->GetPrimaryVertex();
450         if(TMath::Abs(vertex->GetChi2perNDF())>fMaxChi2perNDFforVertex) 
451                 return kFALSE;
452         return kTRUE;
453  }
454
455
456
457 //______________________________________________________
458 void AliSpectraBothEventCuts::PrintCuts()
459 {
460   // print info about event cuts
461   cout << "Event Stats" << endl;
462   if(fHistoCuts)
463   {             
464          cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
465          cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
466          cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
467          cout << " > With good veretx: " << fHistoCuts->GetBinContent(kGoodVtx + 1) << endl;
468          cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
469          cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
470           cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
471   }     
472   cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
473   // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
474   // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
475   cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
476   cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
477   cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
478   cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
479 }
480 //______________________________________________________
481
482 Double_t AliSpectraBothEventCuts::ApplyCentralityPatchAOD049()
483 {
484    //
485    //Apply centrality patch for AOD049 outliers
486    //
487   // if (fCentralityType!="V0M") {
488   //   AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
489   //   return -999.0;
490   // }
491   AliCentrality *centrality = fAOD->GetCentrality();
492    if (!centrality) {
493      AliWarning("Cannot get centrality from AOD event.");
494      return -999.0;
495    }
496    
497    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
498    /*
499      Bool_t isSelRun = kFALSE;
500      Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
501      if(cent<0){
502      Int_t quality = centrality->GetQuality();
503      if(quality<=1){
504      cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
505      } else {
506      Int_t runnum=aodEvent->GetRunNumber();
507      for(Int_t ir=0;ir<5;ir++){
508      if(runnum==selRun[ir]){
509      isSelRun=kTRUE;
510      break;
511      }
512      }
513      if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
514      }
515      }
516    */
517    if(cent>=0.0) {
518      Float_t v0 = 0.0;
519      AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
520      AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
521      v0+=aodV0->GetMTotV0A();
522      v0+=aodV0->GetMTotV0C();
523      if ( (cent==0) && (v0<19500) ) {
524        AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
525        return -999.0;
526      }
527      Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
528      Float_t val = 1.30552 +  0.147931 * v0;
529      
530      Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
531                               120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
532                               92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
533                               68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
534                               51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
535                               37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
536                               26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
537                               19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
538                               12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
539                               13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
540      };
541      
542      if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
543        AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
544        return -999.0;
545      }
546    } else {
547      //force it to be -999. whatever the negative value was
548      cent = -999.;
549    }
550    return cent;
551 }
552
553 //______________________________________________________
554
555
556 Long64_t AliSpectraBothEventCuts::Merge(TCollection* list)
557 {
558   // Merge a list of AliSpectraBothEventCuts objects with this.
559   // Returns the number of merged objects (including this).
560
561   //  AliInfo("Merging");
562
563   if (!list)
564     return 0;
565
566   if (list->IsEmpty())
567     return 1;
568
569   TIterator* iter = list->MakeIterator();
570   TObject* obj;
571
572   // collections of all histograms
573   TList collections;//FIXME we should use only 1 collection
574   TList collections_histoVtxBefSel;
575   TList collections_histoVtxAftSel;
576   TList collections_histoEtaBefSel;
577   TList collections_histoEtaAftSel;
578   TList collections_histoNChAftSel;
579   // TList collections_histoQVectorPos;
580   // TList collections_histoQVectorNeg;
581   TList collections_histoQVector;
582   TList collections_histoEP;
583   TList collections_histoVtxAftSelwithoutZvertexCut;
584   TList collections_histoVtxalltriggerEventswithMCz;
585   TList collections_histoVtxAftSelwithoutZvertexCutusingMCz;                    
586   TList collections_histoRunNumbers;
587   TList collections_histoCentrality;                    
588   TList collections_histoMultiplicty;
589
590   Int_t count = 0;
591
592   while ((obj = iter->Next())) {
593     AliSpectraBothEventCuts* entry = dynamic_cast<AliSpectraBothEventCuts*> (obj);
594     if (entry == 0) 
595       continue;
596
597     TH1I * histo = entry->GetHistoCuts();      
598     collections.Add(histo);
599     TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();      
600     collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
601     TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();      
602     collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
603     TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();      
604     collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
605     TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();      
606     collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
607     TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();      
608     collections_histoNChAftSel.Add(histo_histoNChAftSel);
609     // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();      
610     // collections_histoQVectorPos.Add(histo_histoQVectorPos);
611     // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();      
612     // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
613     TH1F * histo_histoQVector = entry->GetHistoQVector();      
614     collections_histoQVector.Add(histo_histoQVector);
615     TH1F * histo_histoEP = entry->GetHistoEP();      
616     collections_histoEP.Add(histo_histoEP);
617     TH1F* histo_histoVtxAftSelwithoutZvertexCut=entry->GetHistoVtxAftSelwithoutZvertexCut();
618     collections_histoVtxAftSelwithoutZvertexCut.Add(histo_histoVtxAftSelwithoutZvertexCut);
619     TH1F* histo_histoVtxalltriggerEventswithMCz=entry->GetHistoVtxGenerated();
620     collections_histoVtxalltriggerEventswithMCz.Add(histo_histoVtxalltriggerEventswithMCz);
621     
622    TH1F* histo_histoVtxAftSelwithoutZvertexCutusingMCz=entry->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
623     collections_histoVtxAftSelwithoutZvertexCutusingMCz.Add(histo_histoVtxAftSelwithoutZvertexCutusingMCz);     
624     
625     TH1F* histo_histoRunNumbers=entry->GetHistoRunNumbers();
626     collections_histoRunNumbers.Add(histo_histoRunNumbers);
627  
628    TH2F* histo_histoCentrality=entry->GetHistoCentrality();
629   collections_histoCentrality.Add(histo_histoCentrality);               
630
631 TH2F* histo_histoMultiplicty=entry->GetHistoMultiplicty();
632   collections_histoMultiplicty.Add(histo_histoMultiplicty);
633
634
635     count++;
636   }
637   
638   fHistoCuts->Merge(&collections);
639   fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
640   fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
641   fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
642   fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
643   fHistoNChAftSel->Merge(&collections_histoNChAftSel);
644   // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
645   // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
646   fHistoQVector->Merge(&collections_histoQVector);
647   fHistoEP->Merge(&collections_histoEP);
648
649   fHistoVtxAftSelwithoutZvertexCut->Merge(&collections_histoVtxAftSelwithoutZvertexCut);
650   fHistoVtxalltriggerEventswithMCz->Merge(&collections_histoVtxalltriggerEventswithMCz);
651   fHistoVtxAftSelwithoutZvertexCutusingMCz->Merge(&collections_histoVtxAftSelwithoutZvertexCutusingMCz);
652   fHistoRunNumbers->Merge(&collections_histoRunNumbers);
653   fHistoCentrality->Merge(&collections_histoCentrality);
654   fHistoMultiplicty->Merge(&collections_histoMultiplicty);
655
656
657   delete iter;
658
659   return count+1;
660 }
661 //__________________________________________________________________________________________________________
662 void AliSpectraBothEventCuts::SetRunNumberRange(Int_t min, Int_t max)
663 {
664         if(max>min&&min>=0)
665         {
666                  fMinRun=min;                    
667                  fMaxRun=max;
668         }
669 }