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