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