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