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