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                 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
 390                 if(!header) AliFatal("Not a standard AOD");
 391
 392                 if(TMath::Abs(0.8-fetarangeofmultiplicitycut)<0.1)
 393                         Ncharged=header->GetRefMultiplicityComb08();
 394                 else if (TMath::Abs(0.5-fetarangeofmultiplicitycut)<0.1)
 395                         Ncharged=header->GetRefMultiplicityComb05();
 396                 else
 397                         Ncharged=-1;
 398         }
 399         else
 400                 return kFALSE;
 401
 402    fHistoMultiplicty->Fill(0.5,Ncharged);
 403    if(Ncharged>=fMultiplicityCutMin && Ncharged<fMultiplicityCutMax&& Ncharged>0)
 404    {
 405         fHistoMultiplicty->Fill(1.5,Ncharged);
 406         return kTRUE;
 407    }
 408   fHistoCuts->Fill(kVtxCentral);
 409   return kFALSE;
 410 }
 411
 412 //______________________________________________________
 413 Bool_t AliSpectraBothEventCuts::CheckQVectorCut()
 414 {
 415          if(fQVectorCutMin<0.0 && fQVectorCutMax<0.0)return kTRUE;
 416   // Check qvector cut
 417   /// FIXME: Q vector
 418   // //Selection on QVector, before ANY other selection on the event
 419   // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
 420   // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
 421   // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
 422
 423   // Int_t multPos = 0;
 424   // Int_t multNeg = 0;
 425   // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
 426   //   AliAODTrack* aodTrack = fAOD->GetTrack(iT);
 427   //   if (!aodTrack->TestFilterBit(fTrackBits)) continue;
 428   //   if (aodTrack->Eta() >= 0){
 429   //     multPos++;
 430   //     Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
 431   //     Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
 432   //   } else {
 433   //     multNeg++;
 434   //     Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
 435   //     Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
 436   //   }
 437   // }
 438   // Double_t qPos=-999;
 439   // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
 440   // Double_t qNeg=-999;
 441   // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
 442   //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
 443
 444   Double_t qxEPVZERO = 0, qyEPVZERO = 0;
 445   Double_t qVZERO = -999;
 446   Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
 447
 448   qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
 449   if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
 450   fHistoQVector->Fill(qVZERO);
 451   fHistoEP->Fill(psi);
 452
 453   fHistoCuts->Fill(kQVector);
 454   // fHistoQVectorPos->Fill(qPos);
 455   // fHistoQVectorNeg->Fill(qNeg);
 456   return kTRUE;
 457 }
 458 //____________________________________________________________
 459  Bool_t AliSpectraBothEventCuts::CheckVtxChi2perNDF()
 460  {
 461         if(fMaxChi2perNDFforVertex<0)
 462                 return kTRUE;
 463          const AliVVertex * vertex = fAOD->GetPrimaryVertex();
 464         if(TMath::Abs(vertex->GetChi2perNDF())>fMaxChi2perNDFforVertex)
 465                 return kFALSE;
 466         return kTRUE;
 467  }
 468
 469
 470
 471 //______________________________________________________
 472 void AliSpectraBothEventCuts::PrintCuts()
 473 {
 474   // print info about event cuts
 475   cout << "Event Stats" << endl;
 476   if(fHistoCuts)
 477   {
 478          cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
 479          cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
 480          cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
 481          cout << " > With good veretx: " << fHistoCuts->GetBinContent(kGoodVtx + 1) << endl;
 482          cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
 483          cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
 484           cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
 485   }
 486   cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
 487   // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
 488   // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
 489   cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
 490   cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
 491   cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
 492   cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
 493 }
 494 //______________________________________________________
 495
 496 Double_t AliSpectraBothEventCuts::ApplyCentralityPatchAOD049()
 497 {
 498    //
 499    //Apply centrality patch for AOD049 outliers
 500    //
 501   // if (fCentralityType!="V0M") {
 502   //   AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
 503   //   return -999.0;
 504   // }
 505   AliCentrality *centrality = fAOD->GetCentrality();
 506    if (!centrality) {
 507      AliWarning("Cannot get centrality from AOD event.");
 508      return -999.0;
 509    }
 510
 511    Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
 512    /*
 513      Bool_t isSelRun = kFALSE;
 514      Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
 515      if(cent<0){
 516      Int_t quality = centrality->GetQuality();
 517      if(quality<=1){
 518      cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 519      } else {
 520      Int_t runnum=aodEvent->GetRunNumber();
 521      for(Int_t ir=0;ir<5;ir++){
 522      if(runnum==selRun[ir]){
 523      isSelRun=kTRUE;
 524      break;
 525      }
 526      }
 527      if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
 528      }
 529      }
 530    */
 531    if(cent>=0.0) {
 532      Float_t v0 = 0.0;
 533      AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
 534      AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
 535      v0+=aodV0->GetMTotV0A();
 536      v0+=aodV0->GetMTotV0C();
 537      if ( (cent==0) && (v0<19500) ) {
 538        AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
 539        return -999.0;
 540      }
 541      Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
 542      Float_t val = 1.30552 +  0.147931 * v0;
 543
 544      Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
 545                               120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
 546                               92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
 547                               68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
 548                               51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
 549                               37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
 550                               26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
 551                               19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
 552                               12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
 553                               13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
 554      };
 555
 556      if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] )  {
 557        AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
 558        return -999.0;
 559      }
 560    } else {
 561      //force it to be -999. whatever the negative value was
 562      cent = -999.;
 563    }
 564    return cent;
 565 }
 566
 567 //______________________________________________________
 568
 569
 570 Long64_t AliSpectraBothEventCuts::Merge(TCollection* list)
 571 {
 572   // Merge a list of AliSpectraBothEventCuts objects with this.
 573   // Returns the number of merged objects (including this).
 574
 575   //  AliInfo("Merging");
 576
 577   if (!list)
 578     return 0;
 579
 580   if (list->IsEmpty())
 581     return 1;
 582
 583   TIterator* iter = list->MakeIterator();
 584   TObject* obj;
 585
 586   // collections of all histograms
 587   TList collections;//FIXME we should use only 1 collection
 588   TList collections_histoVtxBefSel;
 589   TList collections_histoVtxAftSel;
 590   TList collections_histoEtaBefSel;
 591   TList collections_histoEtaAftSel;
 592   TList collections_histoNChAftSel;
 593   // TList collections_histoQVectorPos;
 594   // TList collections_histoQVectorNeg;
 595   TList collections_histoQVector;
 596   TList collections_histoEP;
 597   TList collections_histoVtxAftSelwithoutZvertexCut;
 598   TList collections_histoVtxalltriggerEventswithMCz;
 599   TList collections_histoVtxAftSelwithoutZvertexCutusingMCz;
 600   TList collections_histoRunNumbers;
 601   TList collections_histoCentrality;
 602   TList collections_histoMultiplicty;
 603
 604   Int_t count = 0;
 605
 606   while ((obj = iter->Next())) {
 607     AliSpectraBothEventCuts* entry = dynamic_cast<AliSpectraBothEventCuts*> (obj);
 608     if (entry == 0)
 609       continue;
 610
 611     TH1I * histo = entry->GetHistoCuts();
 612     collections.Add(histo);
 613     TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
 614     collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
 615     TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
 616     collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
 617     TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
 618     collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
 619     TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
 620     collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
 621     TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
 622     collections_histoNChAftSel.Add(histo_histoNChAftSel);
 623     // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
 624     // collections_histoQVectorPos.Add(histo_histoQVectorPos);
 625     // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
 626     // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
 627     TH1F * histo_histoQVector = entry->GetHistoQVector();
 628     collections_histoQVector.Add(histo_histoQVector);
 629     TH1F * histo_histoEP = entry->GetHistoEP();
 630     collections_histoEP.Add(histo_histoEP);
 631     TH1F* histo_histoVtxAftSelwithoutZvertexCut=entry->GetHistoVtxAftSelwithoutZvertexCut();
 632     collections_histoVtxAftSelwithoutZvertexCut.Add(histo_histoVtxAftSelwithoutZvertexCut);
 633     TH1F* histo_histoVtxalltriggerEventswithMCz=entry->GetHistoVtxGenerated();
 634     collections_histoVtxalltriggerEventswithMCz.Add(histo_histoVtxalltriggerEventswithMCz);
 635
 636    TH1F* histo_histoVtxAftSelwithoutZvertexCutusingMCz=entry->GetHistoVtxAftSelwithoutZvertexCutusingMCz();
 637     collections_histoVtxAftSelwithoutZvertexCutusingMCz.Add(histo_histoVtxAftSelwithoutZvertexCutusingMCz);
 638
 639     TH1F* histo_histoRunNumbers=entry->GetHistoRunNumbers();
 640     collections_histoRunNumbers.Add(histo_histoRunNumbers);
 641
 642    TH2F* histo_histoCentrality=entry->GetHistoCentrality();
 643   collections_histoCentrality.Add(histo_histoCentrality);
 644
 645 TH2F* histo_histoMultiplicty=entry->GetHistoMultiplicty();
 646   collections_histoMultiplicty.Add(histo_histoMultiplicty);
 647
 648
 649     count++;
 650   }
 651
 652   fHistoCuts->Merge(&collections);
 653   fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
 654   fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
 655   fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
 656   fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
 657   fHistoNChAftSel->Merge(&collections_histoNChAftSel);
 658   // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
 659   // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
 660   fHistoQVector->Merge(&collections_histoQVector);
 661   fHistoEP->Merge(&collections_histoEP);
 662
 663   fHistoVtxAftSelwithoutZvertexCut->Merge(&collections_histoVtxAftSelwithoutZvertexCut);
 664   fHistoVtxalltriggerEventswithMCz->Merge(&collections_histoVtxalltriggerEventswithMCz);
 665   fHistoVtxAftSelwithoutZvertexCutusingMCz->Merge(&collections_histoVtxAftSelwithoutZvertexCutusingMCz);
 666   fHistoRunNumbers->Merge(&collections_histoRunNumbers);
 667   fHistoCentrality->Merge(&collections_histoCentrality);
 668   fHistoMultiplicty->Merge(&collections_histoMultiplicty);
 669
 670
 671   delete iter;
 672
 673   return count+1;
 674 }
 675 //__________________________________________________________________________________________________________
 676 void AliSpectraBothEventCuts::SetRunNumberRange(Int_t min, Int_t max)
 677 {
 678         if(max>min&&min>=0)
 679         {
 680                  fMinRun=min;
 681                  fMaxRun=max;
 682         }
 683 }