PWGGA/GammaConv/AliConvEventCuts.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Authors: Svein Lindal, Daniel Lohner                                   *
   5  * Version 1.0                                                            *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////
  17 //---------------------------------------------
  18 // Class handling all kinds of selection cuts for
  19 // Gamma Conversion analysis
  20 //---------------------------------------------
  21 ////////////////////////////////////////////////
  22
  23 #include "AliConvEventCuts.h"
  24
  25 #include "AliAODTrack.h"
  26 #include "AliESDtrack.h"
  27 #include "AliAnalysisManager.h"
  28 #include "AliInputEventHandler.h"
  29 #include "AliMCEventHandler.h"
  30 #include "AliAODHandler.h"
  31 #include "TH1.h"
  32 #include "TH2.h"
  33 #include "TF1.h"
  34 #include "AliStack.h"
  35 #include "TObjString.h"
  36 #include "AliAODEvent.h"
  37 #include "AliESDEvent.h"
  38 #include "AliCentrality.h"
  39 #include "TList.h"
  40 #include "TFile.h"
  41 #include "AliLog.h"
  42 #include "AliGenCocktailEventHeader.h"
  43 #include "AliGenDPMjetEventHeader.h"
  44 #include "AliGenPythiaEventHeader.h"
  45 #include "AliGenHijingEventHeader.h"
  46 #include "AliTriggerAnalysis.h"
  47 #include "AliV0ReaderV1.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODMCHeader.h"
  50
  51 class iostream;
  52
  53 using namespace std;
  54
  55 ClassImp(AliConvEventCuts)
  56
  57
  58 const char* AliConvEventCuts::fgkCutNames[AliConvEventCuts::kNCuts] = {
  59    "HeavyIon",                                          //0
  60    "CentralityMin",                                     //1
  61    "CentralityMax",                                     //2
  62    "SelectSpecialTrigger",                      //3
  63    "SelectSpecialSubTriggerClass",      //4
  64    "RemovePileUp",                                      //5
  65    "RejectExtraSignals",                        //6
  66 };
  67
  68
  69 //________________________________________________________________________
  70 AliConvEventCuts::AliConvEventCuts(const char *name,const char *title) :
  71    AliAnalysisCuts(name,title),
  72    fHistograms(NULL),
  73    fHeaderList(NULL),
  74    fEventQuality(-1),
  75    fIsHeavyIon(0),
  76    fDetectorCentrality(0),
  77    fModCentralityClass(0),
  78    fMaxVertexZ(10),
  79    fCentralityMin(0),
  80    fCentralityMax(0),
  81    fMultiplicityMethod(0),
  82    fSpecialTrigger(0),
  83    fSpecialSubTrigger(0),
  84    fRemovePileUp(kFALSE),
  85    fRejectExtraSignals(0),
  86    fOfflineTriggerMask(0),
  87    fHasV0AND(kTRUE),
  88    fIsSDDFired(kTRUE),
  89    fRandom(0),
  90    fnHeaders(0),
  91    fNotRejectedStart(NULL),
  92    fNotRejectedEnd(NULL),
  93    fGeneratorNames(NULL),
  94    fCutString(NULL),
  95    fUtils(NULL),
  96    fEtaShift(0.0),
  97    fDoEtaShift(kFALSE),
  98    fDoReweightHistoMCPi0(kFALSE),
  99    fDoReweightHistoMCEta(kFALSE),
 100    fDoReweightHistoMCK0s(kFALSE),
 101    fPathTrFReweighting(""),
 102    fNameHistoReweightingPi0(""),
 103    fNameHistoReweightingEta(""),
 104    fNameHistoReweightingK0s(""),
 105    fNameFitDataPi0(""),
 106    fNameFitDataEta(""),
 107    fNameFitDataK0s(""),
 108    fHistoEventCuts(NULL),
 109    hCentrality(NULL),
 110    hCentralityVsNumberOfPrimaryTracks(NULL),
 111    hVertexZ(NULL),
 112    hTriggerClass(NULL),
 113    hTriggerClassSelected(NULL),
 114    hReweightMCHistPi0(NULL),
 115    hReweightMCHistEta(NULL),
 116    hReweightMCHistK0s(NULL),
 117    fFitDataPi0(NULL),
 118    fFitDataEta(NULL),
 119    fFitDataK0s(NULL),
 120    fAddedSignalPDGCode(0),
 121    fPreSelCut(kFALSE),
 122    fTriggerSelectedManually(kFALSE),
 123    fSpecialTriggerName(""),
 124    fSpecialSubTriggerName(""),
 125    fNSpecialSubTriggerOptions(0),
 126    fV0ReaderName("")
 127
 128 {
 129    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
 130    fCutString=new TObjString((GetCutNumber()).Data());
 131
 132    fUtils = new AliAnalysisUtils();
 133    //if you do not want to apply the cut on the distance between the SPD and TRK vertex:
 134    //fUtils->SetCutOnZVertexSPD(kFALSE);
 135
 136
 137 }
 138
 139 //________________________________________________________________________
 140 AliConvEventCuts::AliConvEventCuts(const AliConvEventCuts &ref) :
 141    AliAnalysisCuts(ref),
 142    fHistograms(NULL),
 143    fHeaderList(ref.fHeaderList),
 144    fEventQuality(ref.fEventQuality),
 145    fIsHeavyIon(ref.fIsHeavyIon),
 146    fDetectorCentrality(ref.fDetectorCentrality),
 147    fModCentralityClass(ref.fModCentralityClass),
 148    fMaxVertexZ(ref.fMaxVertexZ),
 149    fCentralityMin(ref.fCentralityMin),
 150    fCentralityMax(ref.fCentralityMax),
 151    fMultiplicityMethod(ref.fMultiplicityMethod),
 152    fSpecialTrigger(ref.fSpecialTrigger),
 153    fSpecialSubTrigger(ref.fSpecialSubTrigger),
 154    fRemovePileUp(ref.fRemovePileUp),
 155    fRejectExtraSignals(ref.fRejectExtraSignals),
 156    fOfflineTriggerMask(ref.fOfflineTriggerMask),
 157    fHasV0AND(ref.fHasV0AND),
 158    fIsSDDFired(ref.fIsSDDFired),
 159    fRandom(ref.fRandom),
 160    fnHeaders(ref.fnHeaders),
 161    fNotRejectedStart(NULL),
 162    fNotRejectedEnd(NULL),
 163    fGeneratorNames(ref.fGeneratorNames),
 164    fCutString(NULL),
 165    fUtils(NULL),
 166    fEtaShift(ref.fEtaShift),
 167    fDoEtaShift(ref.fDoEtaShift),
 168    fDoReweightHistoMCPi0(ref.fDoReweightHistoMCPi0),
 169    fDoReweightHistoMCEta(ref.fDoReweightHistoMCEta),
 170    fDoReweightHistoMCK0s(ref.fDoReweightHistoMCK0s),
 171    fPathTrFReweighting(ref.fPathTrFReweighting),
 172    fNameHistoReweightingPi0(ref.fNameHistoReweightingPi0),
 173    fNameHistoReweightingEta(ref.fNameHistoReweightingEta),
 174    fNameHistoReweightingK0s(ref.fNameHistoReweightingK0s),
 175    fNameFitDataPi0(ref.fNameFitDataPi0),
 176    fNameFitDataEta(ref.fNameFitDataEta),
 177    fNameFitDataK0s(ref.fNameFitDataK0s),
 178    fHistoEventCuts(NULL),
 179    hCentrality(NULL),
 180    hCentralityVsNumberOfPrimaryTracks(NULL),
 181    hVertexZ(NULL),
 182    hTriggerClass(NULL),
 183    hTriggerClassSelected(NULL),
 184    hReweightMCHistPi0(ref.hReweightMCHistPi0),
 185    hReweightMCHistEta(ref.hReweightMCHistEta),
 186    hReweightMCHistK0s(ref.hReweightMCHistK0s),
 187    fFitDataPi0(ref.fFitDataPi0),
 188    fFitDataEta(ref.fFitDataEta),
 189    fFitDataK0s(ref.fFitDataK0s),
 190    fAddedSignalPDGCode(ref.fAddedSignalPDGCode),
 191    fPreSelCut(ref.fPreSelCut),
 192    fTriggerSelectedManually(ref.fTriggerSelectedManually),
 193    fSpecialTriggerName(ref.fSpecialTriggerName),
 194    fSpecialSubTriggerName(ref.fSpecialSubTriggerName),
 195    fNSpecialSubTriggerOptions(ref.fNSpecialSubTriggerOptions),
 196    fV0ReaderName(ref.fV0ReaderName)
 197 {
 198    // Copy Constructor
 199    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
 200    fCutString=new TObjString((GetCutNumber()).Data());
 201    fUtils = new AliAnalysisUtils();
 202    // dont copy histograms (if you like histograms, call InitCutHistograms())
 203
 204 }
 205
 206
 207 //________________________________________________________________________
 208 AliConvEventCuts::~AliConvEventCuts() {
 209    // Destructor
 210    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
 211    // if(fHistograms)
 212    //    delete fHistograms;
 213    // fHistograms = NULL;
 214    if(fCutString != NULL){
 215       delete fCutString;
 216       fCutString = NULL;
 217    }
 218    if(fNotRejectedStart){
 219       delete[] fNotRejectedStart;
 220       fNotRejectedStart = NULL;
 221    }
 222    if(fNotRejectedEnd){
 223       delete[] fNotRejectedEnd;
 224       fNotRejectedEnd = NULL;
 225    }
 226    if(fGeneratorNames){
 227       delete[] fGeneratorNames;
 228       fGeneratorNames = NULL;
 229    }
 230    if(fUtils){
 231      delete fUtils;
 232      fUtils = NULL;
 233    }
 234
 235 }
 236
 237 //________________________________________________________________________
 238 void AliConvEventCuts::InitCutHistograms(TString name, Bool_t preCut){
 239
 240    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
 241    TH1::AddDirectory(kFALSE);
 242
 243    if(fHistograms != NULL){
 244       delete fHistograms;
 245       fHistograms=NULL;
 246    }
 247    if(fHistograms==NULL){
 248       fHistograms=new TList();
 249       fHistograms->SetOwner(kTRUE);
 250       if(name=="")fHistograms->SetName(Form("ConvEventCuts_%s",GetCutNumber().Data()));
 251       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
 252    }
 253
 254    if (hReweightMCHistPi0){
 255       hReweightMCHistPi0->SetName("MCInputForWeightingPi0");
 256       fHistograms->Add(hReweightMCHistPi0);
 257    }
 258    if (hReweightMCHistEta){
 259       hReweightMCHistEta->SetName("MCInputForWeightingEta");
 260       fHistograms->Add(hReweightMCHistEta);
 261    }
 262    if (hReweightMCHistK0s){
 263       hReweightMCHistK0s->SetName("MCInputForWeightingK0s");
 264       fHistograms->Add(hReweightMCHistK0s);
 265    }
 266
 267    hCentrality=new TH1F(Form("Centrality %s",GetCutNumber().Data()),"Centrality",100,0,100);
 268    fHistograms->Add(hCentrality);
 269    hCentralityVsNumberOfPrimaryTracks=new TH2F(Form("Centrality vs Primary Tracks %s",GetCutNumber().Data()),"Centrality vs Primary Tracks ",100,0,100,4000,0,4000);
 270    fHistograms->Add(hCentralityVsNumberOfPrimaryTracks);
 271
 272    // Event Cuts and Info
 273    if(preCut){
 274       fHistoEventCuts=new TH1F(Form("ESD_EventCuts %s",GetCutNumber().Data()),"Event Cuts",7,-0.5,6.5);
 275       fHistoEventCuts->GetXaxis()->SetBinLabel(1,"in");
 276       fHistoEventCuts->GetXaxis()->SetBinLabel(2,"OfflineTrigger");
 277       fHistoEventCuts->GetXaxis()->SetBinLabel(3,"nvtxcontr");
 278       fHistoEventCuts->GetXaxis()->SetBinLabel(4,"VertexZ");
 279       fHistoEventCuts->GetXaxis()->SetBinLabel(5,"pileup");
 280       fHistoEventCuts->GetXaxis()->SetBinLabel(6,"centrsel");
 281       fHistoEventCuts->GetXaxis()->SetBinLabel(7,"out");
 282       fHistograms->Add(fHistoEventCuts);
 283
 284       hVertexZ=new TH1F(Form("VertexZ %s",GetCutNumber().Data()),"VertexZ",1000,-50,50);
 285       fHistograms->Add(hVertexZ);
 286
 287       hTriggerClass= new TH1F(Form("OfflineTrigger %s",GetCutNumber().Data()),"OfflineTrigger",35,-0.5,34.5);
 288       hTriggerClass->GetXaxis()->SetBinLabel( 1,"kMB");
 289       hTriggerClass->GetXaxis()->SetBinLabel( 2,"kINT7");
 290       hTriggerClass->GetXaxis()->SetBinLabel( 3,"kMUON");
 291       hTriggerClass->GetXaxis()->SetBinLabel( 4,"kHighMult");
 292       hTriggerClass->GetXaxis()->SetBinLabel( 5,"kKEMC1");
 293       hTriggerClass->GetXaxis()->SetBinLabel( 6,"kCINT5");
 294       hTriggerClass->GetXaxis()->SetBinLabel( 7,"kCMUS5/kMUSPB");
 295       hTriggerClass->GetXaxis()->SetBinLabel( 8,"kMUSH7/kMUSHPB");
 296       hTriggerClass->GetXaxis()->SetBinLabel( 9,"kMUL7/kMuonLikePB");
 297       hTriggerClass->GetXaxis()->SetBinLabel(10,"kMUU7/kMuonUnlikePB");
 298       hTriggerClass->GetXaxis()->SetBinLabel(11,"kEMC7/kEMC8");
 299       hTriggerClass->GetXaxis()->SetBinLabel(12,"kMUS7");
 300       hTriggerClass->GetXaxis()->SetBinLabel(13,"kPHI1");
 301       hTriggerClass->GetXaxis()->SetBinLabel(14,"kPHI7/kPHI8/kPHOSPb");
 302       hTriggerClass->GetXaxis()->SetBinLabel(15,"kEMCEJE");
 303       hTriggerClass->GetXaxis()->SetBinLabel(16,"kEMCEGA");
 304       hTriggerClass->GetXaxis()->SetBinLabel(17,"kCentral");
 305       hTriggerClass->GetXaxis()->SetBinLabel(18,"kSemiCentral");
 306       hTriggerClass->GetXaxis()->SetBinLabel(19,"kDG5");
 307       hTriggerClass->GetXaxis()->SetBinLabel(20,"kZED");
 308       hTriggerClass->GetXaxis()->SetBinLabel(21,"kSPI7/kSPI");
 309       hTriggerClass->GetXaxis()->SetBinLabel(22,"kINT8");
 310       hTriggerClass->GetXaxis()->SetBinLabel(23,"kMuonSingleLowPt8");
 311       hTriggerClass->GetXaxis()->SetBinLabel(24,"kMuonSingleHighPt8");
 312       hTriggerClass->GetXaxis()->SetBinLabel(25,"kMuonLikeLowPt8");
 313       hTriggerClass->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
 314       hTriggerClass->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
 315       hTriggerClass->GetXaxis()->SetBinLabel(28,"kUserDefined");
 316       hTriggerClass->GetXaxis()->SetBinLabel(29,"kTRD");
 317       hTriggerClass->GetXaxis()->SetBinLabel(30,"kFastOnly");
 318       hTriggerClass->GetXaxis()->SetBinLabel(31,"kAnyINT");
 319       hTriggerClass->GetXaxis()->SetBinLabel(32,"kAny");
 320       hTriggerClass->GetXaxis()->SetBinLabel(33,"V0AND");
 321       hTriggerClass->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
 322       hTriggerClass->GetXaxis()->SetBinLabel(35,"failed Physics Selection");
 323       fHistograms->Add(hTriggerClass);
 324    }
 325    if(!preCut){
 326       hTriggerClassSelected= new TH1F(Form("OfflineTriggerSelected %s",GetCutNumber().Data()),"OfflineTriggerSelected",34,-0.5,33.5);
 327       hTriggerClassSelected->GetXaxis()->SetBinLabel( 1,"kMB");
 328       hTriggerClassSelected->GetXaxis()->SetBinLabel( 2,"kINT7");
 329       hTriggerClassSelected->GetXaxis()->SetBinLabel( 3,"kMUON");
 330       hTriggerClassSelected->GetXaxis()->SetBinLabel( 4,"kHighMult");
 331       hTriggerClassSelected->GetXaxis()->SetBinLabel( 5,"kKEMC1");
 332       hTriggerClassSelected->GetXaxis()->SetBinLabel( 6,"kCINT5");
 333       hTriggerClassSelected->GetXaxis()->SetBinLabel( 7,"kCMUS5/kMUSPB");
 334       hTriggerClassSelected->GetXaxis()->SetBinLabel( 8,"kMUSH7/kMUSHPB");
 335       hTriggerClassSelected->GetXaxis()->SetBinLabel( 9,"kMUL7/kMuonLikePB");
 336       hTriggerClassSelected->GetXaxis()->SetBinLabel(10,"kMUU7/kMuonUnlikePB");
 337       hTriggerClassSelected->GetXaxis()->SetBinLabel(11,"kEMC7/kEMC8");
 338       hTriggerClassSelected->GetXaxis()->SetBinLabel(12,"kMUS7");
 339       hTriggerClassSelected->GetXaxis()->SetBinLabel(13,"kPHI1");
 340       hTriggerClassSelected->GetXaxis()->SetBinLabel(14,"kPHI7/kPHI8/kPHOSPb");
 341       hTriggerClassSelected->GetXaxis()->SetBinLabel(15,"kEMCEJE");
 342       hTriggerClassSelected->GetXaxis()->SetBinLabel(16,"kEMCEGA");
 343       hTriggerClassSelected->GetXaxis()->SetBinLabel(17,"kCentral");
 344       hTriggerClassSelected->GetXaxis()->SetBinLabel(18,"kSemiCentral");
 345       hTriggerClassSelected->GetXaxis()->SetBinLabel(19,"kDG5");
 346       hTriggerClassSelected->GetXaxis()->SetBinLabel(20,"kZED");
 347       hTriggerClassSelected->GetXaxis()->SetBinLabel(21,"kSPI7/kSPI");
 348       hTriggerClassSelected->GetXaxis()->SetBinLabel(22,"kINT8");
 349       hTriggerClassSelected->GetXaxis()->SetBinLabel(23,"kMuonSingleLowPt8");
 350       hTriggerClassSelected->GetXaxis()->SetBinLabel(24,"kMuonSingleHighPt8");
 351       hTriggerClassSelected->GetXaxis()->SetBinLabel(25,"kMuonLikeLowPt8");
 352       hTriggerClassSelected->GetXaxis()->SetBinLabel(26,"kMuonUnlikeLowPt8");
 353       hTriggerClassSelected->GetXaxis()->SetBinLabel(27,"kMuonUnlikeLowPt0");
 354       hTriggerClassSelected->GetXaxis()->SetBinLabel(28,"kUserDefined");
 355       hTriggerClassSelected->GetXaxis()->SetBinLabel(29,"kTRD");
 356       hTriggerClassSelected->GetXaxis()->SetBinLabel(30,"kFastOnly");
 357       hTriggerClassSelected->GetXaxis()->SetBinLabel(31,"kAnyINT");
 358       hTriggerClassSelected->GetXaxis()->SetBinLabel(32,"kAny");
 359       hTriggerClassSelected->GetXaxis()->SetBinLabel(33,"V0AND");
 360       hTriggerClassSelected->GetXaxis()->SetBinLabel(34,"NOT kFastOnly");
 361       fHistograms->Add(hTriggerClassSelected);
 362
 363    }
 364    TH1::AddDirectory(kTRUE);
 365 }
 366
 367 ///________________________________________________________________________
 368 Bool_t AliConvEventCuts::EventIsSelected(AliVEvent *fInputEvent, AliVEvent *fMCEvent){
 369    // Process Event Selection
 370
 371    Int_t cutindex=0;
 372    if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 373    cutindex++;
 374
 375    // Check for MC event
 376    if(fMCEvent && fInputEvent->IsA()==AliESDEvent::Class()){
 377       // Check if MC event is correctly loaded
 378       AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 379       if (!mcHandler){
 380          fEventQuality = 2;
 381          return kFALSE;
 382       }
 383       if (!mcHandler->InitOk() ){
 384          fEventQuality = 2;
 385          return kFALSE;
 386       }
 387       if (!mcHandler->TreeK() ){
 388          fEventQuality = 2;
 389          return kFALSE;
 390       }
 391       if (!mcHandler->TreeTR() ) {
 392          fEventQuality = 2;
 393          return kFALSE;
 394       }
 395    }
 396
 397    // Event Trigger
 398 //    cout << "before event trigger" << endl;
 399    if(!IsTriggerSelected(fInputEvent)){
 400       if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 401       fEventQuality = 3;
 402       return kFALSE;
 403    }
 404    cutindex++;
 405
 406    if(fInputEvent->IsA()==AliESDEvent::Class()){
 407       AliTriggerAnalysis fTriggerAnalysis;// = new AliTriggerAnalysis;
 408       fHasV0AND = fTriggerAnalysis.IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
 409       if(fHasV0AND&&hTriggerClass)hTriggerClass->Fill(32);
 410    }
 411 //   cout << "event number " << ((AliESDEvent*)fInputEvent)->GetEventNumberInFile() << " entered"<< endl;
 412
 413
 414    // Number of Contributors Cut
 415    if(GetNumberOfContributorsVtx(fInputEvent)<=0) {
 416       if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 417       fEventQuality = 5;
 418       return kFALSE;
 419    }
 420    cutindex++;
 421
 422    // Z Vertex Position Cut
 423    if(!VertexZCut(fInputEvent)){
 424       if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 425       fEventQuality = 4;
 426       return kFALSE;
 427    }
 428    cutindex++;
 429
 430    // Pile Up Rejection
 431
 432    if(fRemovePileUp){
 433       if(fInputEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
 434          if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 435          fEventQuality = 6;
 436          return kFALSE;
 437       }
 438    }
 439    cutindex++;
 440
 441    // Centrality Selection
 442    if(!IsCentralitySelected(fInputEvent,fMCEvent)){
 443       if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 444       fEventQuality = 1;
 445       return kFALSE;
 446    }
 447    cutindex++;
 448
 449    // Fill Event Histograms
 450    if(fHistoEventCuts)fHistoEventCuts->Fill(cutindex);
 451    if(hVertexZ)hVertexZ->Fill(fInputEvent->GetPrimaryVertex()->GetZ());
 452    if(hCentrality)hCentrality->Fill(GetCentrality(fInputEvent));
 453    if(hCentralityVsNumberOfPrimaryTracks)
 454       hCentralityVsNumberOfPrimaryTracks->Fill(GetCentrality(fInputEvent),
 455                                                ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
 456                                                 ->GetTask(fV0ReaderName.Data()))->GetNumberOfPrimaryTracks());
 457    fEventQuality = 0;
 458    return kTRUE;
 459 }
 460
 461 ///________________________________________________________________________
 462 Bool_t AliConvEventCuts::UpdateCutString() {
 463         ///Update the cut string (if it has been created yet)
 464
 465         if(fCutString && fCutString->GetString().Length() == kNCuts) {
 466                 fCutString->SetString(GetCutNumber());
 467         } else {
 468                 return kFALSE;
 469         }
 470         return kTRUE;
 471 }
 472
 473 ///________________________________________________________________________
 474 void AliConvEventCuts::LoadReweightingHistosMCFromFile() {
 475
 476         AliInfo("Entering loading of histograms for weighting");
 477         TFile *f = TFile::Open(fPathTrFReweighting.Data());
 478         if(!f){
 479                 AliError(Form("file for weighting %s not found",fPathTrFReweighting.Data()));
 480                 return;
 481         }
 482         if (fNameHistoReweightingPi0.CompareTo("") != 0 && fDoReweightHistoMCPi0 ){
 483                 cout << "I have to find: " <<  fNameHistoReweightingPi0.Data() << endl;
 484                 TH1D *hReweightMCHistPi0temp = (TH1D*)f->Get(fNameHistoReweightingPi0.Data());
 485                 hReweightMCHistPi0 = new TH1D(*hReweightMCHistPi0temp);
 486                 if (hReweightMCHistPi0) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingPi0.Data(),fPathTrFReweighting.Data() ));
 487                 else AliWarning(Form("%s not found in %s", fNameHistoReweightingPi0.Data() ,fPathTrFReweighting.Data()));
 488                 hReweightMCHistPi0->SetDirectory(0);
 489         }
 490         if (fNameFitDataPi0.CompareTo("") != 0 && fDoReweightHistoMCPi0 ){
 491                 cout << "I have to find: " <<  fNameFitDataPi0.Data() << endl;
 492                 TF1 *fFitDataPi0temp = (TF1*)f->Get(fNameFitDataPi0.Data());
 493                 fFitDataPi0 = new TF1(*fFitDataPi0temp);
 494                 if (fFitDataPi0) AliInfo(Form("%s has been loaded from %s", fNameFitDataPi0.Data(),fPathTrFReweighting.Data() ));
 495                 else AliWarning(Form("%s not found in %s",fPathTrFReweighting.Data(), fNameFitDataPi0.Data() ));
 496         }
 497
 498         if (fNameHistoReweightingEta.CompareTo("") != 0 && fDoReweightHistoMCEta){
 499                 cout << "I have to find: " <<  fNameHistoReweightingEta.Data() << endl;
 500                 TH1D *hReweightMCHistEtatemp = (TH1D*)f->Get(fNameHistoReweightingEta.Data());
 501                 hReweightMCHistEta = new TH1D(*hReweightMCHistEtatemp);
 502                 if (hReweightMCHistEta) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingEta.Data(),fPathTrFReweighting.Data() ));
 503                 else AliWarning(Form("%s not found in %s", fNameHistoReweightingEta.Data(),fPathTrFReweighting.Data() ));
 504                 hReweightMCHistEta->SetDirectory(0);
 505         }
 506
 507         if (fNameFitDataEta.CompareTo("") != 0 && fDoReweightHistoMCEta){
 508                 cout << "I have to find: " <<  fNameFitDataEta.Data() << endl;
 509                 TF1 *fFitDataEtatemp = (TF1*)f->Get(fNameFitDataEta.Data());
 510                 fFitDataEta = new TF1(*fFitDataEtatemp);
 511                 if (fFitDataEta) AliInfo(Form("%s has been loaded from %s", fNameFitDataEta.Data(),fPathTrFReweighting.Data() ));
 512                 else AliWarning(Form("%s not found in %s", fNameFitDataEta.Data(),fPathTrFReweighting.Data() ));
 513
 514         }
 515         if (fNameHistoReweightingK0s.CompareTo("") != 0 && fDoReweightHistoMCK0s){
 516                 cout << "I have to find: " <<  fNameHistoReweightingK0s.Data() << endl;
 517                 TH1D *hReweightMCHistK0stemp = (TH1D*)f->Get(fNameHistoReweightingK0s.Data());
 518                 hReweightMCHistK0s = new TH1D(*hReweightMCHistK0stemp);
 519                 if (hReweightMCHistK0s) AliInfo(Form("%s has been loaded from %s", fNameHistoReweightingK0s.Data(),fPathTrFReweighting.Data() ));
 520                 else AliWarning(Form("%s not found in %s", fNameHistoReweightingK0s.Data(),fPathTrFReweighting.Data() ));
 521                 hReweightMCHistK0s->SetDirectory(0);
 522         }
 523
 524         if (fNameFitDataK0s.CompareTo("") != 0 && fDoReweightHistoMCK0s){
 525                 cout << "I have to find: " <<  fNameFitDataK0s.Data() << endl;
 526                 TF1 *fFitDataK0stemp = (TF1*)f->Get(fNameFitDataK0s.Data());
 527                 fFitDataK0s = new TF1(*fFitDataK0stemp);
 528                 if (fFitDataK0s) AliInfo(Form("%s has been loaded from %s", fNameFitDataK0s.Data(),fPathTrFReweighting.Data() ));
 529                 else AliWarning(Form("%s not found in %s", fNameFitDataK0s.Data(),fPathTrFReweighting.Data() ));
 530         }
 531         f->Close();
 532         delete f;
 533 }
 534
 535
 536 ///________________________________________________________________________
 537 Bool_t AliConvEventCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 538         // Initialize Cuts from a given Cut string
 539         if(fDoReweightHistoMCPi0 || fDoReweightHistoMCEta || fDoReweightHistoMCK0s) {
 540                 AliInfo("Weighting was enabled");
 541                 LoadReweightingHistosMCFromFile();
 542         }
 543
 544         AliInfo(Form("Set Event Cut Number: %s",analysisCutSelection.Data()));
 545         if(analysisCutSelection.Length()!=kNCuts) {
 546                 AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 547                 return kFALSE;
 548         }
 549         if(!analysisCutSelection.IsDigit()){
 550                 AliError("Cut selection contains characters");
 551                 return kFALSE;
 552         }
 553
 554         if (fV0ReaderName.CompareTo("") == 0){
 555                 fV0ReaderName = "V0ReaderV1";
 556         }
 557         const char *cutSelection = analysisCutSelection.Data();
 558         #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 559         for(Int_t ii=0;ii<kNCuts;ii++){
 560                 ASSIGNARRAY(ii);
 561         }
 562
 563         // Set Individual Cuts
 564         for(Int_t ii=0;ii<kNCuts;ii++){
 565                 if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 566         }
 567
 568         PrintCutsWithValues();
 569
 570         return kTRUE;
 571 }
 572
 573 ///________________________________________________________________________
 574 Bool_t AliConvEventCuts::SetCut(cutIds cutID, const Int_t value) {
 575         ///Set individual cut ID
 576
 577         switch (cutID) {
 578         case kremovePileUp:
 579                 if( SetRemovePileUp(value)) {
 580                         fCuts[kremovePileUp] = value;
 581                         UpdateCutString();
 582                         return kTRUE;
 583                 } else return kFALSE;
 584         case kSelectSpecialTriggerAlias:
 585                 if( SetSelectSpecialTrigger(value)) {
 586                         fCuts[kSelectSpecialTriggerAlias] = value;
 587                         UpdateCutString();
 588                         return kTRUE;
 589                 } else return kFALSE;
 590         case kSelectSubTriggerClass:
 591                 if( SetSelectSubTriggerClass(value)) {
 592                         fCuts[kSelectSubTriggerClass] = value;
 593                         UpdateCutString();
 594                         return kTRUE;
 595                 } else return kFALSE;
 596         case kisHeavyIon:
 597                 if( SetIsHeavyIon(value)) {
 598                         fCuts[kisHeavyIon] = value;
 599                         UpdateCutString();
 600                         return kTRUE;
 601                 } else return kFALSE;
 602         case kCentralityMin:
 603                 if( SetCentralityMin(value)) {
 604                         fCuts[kCentralityMin] = value;
 605                         UpdateCutString();
 606                         return kTRUE;
 607                 } else return kFALSE;
 608         case kCentralityMax:
 609                 if( SetCentralityMax(value)) {
 610                         fCuts[kCentralityMax] = value;
 611                         UpdateCutString();
 612                         return kTRUE;
 613                 } else return kFALSE;
 614         case kExtraSignals:
 615                 if( SetRejectExtraSignalsCut(value)) {
 616                         fCuts[kExtraSignals] = value;
 617                         UpdateCutString();
 618                         return kTRUE;
 619                 } else return kFALSE;
 620         case kNCuts:
 621                 AliError("Cut id out of range");
 622                 return kFALSE;
 623         }
 624
 625         AliError("Cut id %d not recognized");
 626         return kFALSE;
 627 }
 628
 629 ///________________________________________________________________________
 630 void AliConvEventCuts::PrintCuts() {
 631         // Print out current Cut Selection
 632         for(Int_t ic = 0; ic < kNCuts; ic++) {
 633                 printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 634         }
 635 }
 636
 637 void AliConvEventCuts::PrintCutsWithValues() {
 638    // Print out current Cut Selection with value
 639         printf("\nEvent cutnumber \n");
 640         for(Int_t ic = 0; ic < kNCuts; ic++) {
 641                 printf("%d",fCuts[ic]);
 642         }
 643         printf("\n\n");
 644
 645         if (fIsHeavyIon == 0) {
 646                 printf("Running in pp mode \n");
 647                 if (fSpecialTrigger == 0){
 648                         printf("\t only events triggered by V0OR will be analysed \n");
 649                 } else if (fSpecialTrigger == 1){
 650                         printf("\t only events triggered by V0AND will be analysed \n");
 651                 } else if (fSpecialTrigger == 2){
 652                         printf("\t only events where SDD was present will be analysed \n");
 653                 } else if (fSpecialTrigger == 3){
 654                         printf("\t only events where SDD was present will be analysed and triggered by VOAND\n");
 655                 } else if (fSpecialTrigger > 3){
 656                         printf("\t only events triggered by %s %s\n", fSpecialTriggerName.Data(), fSpecialSubTriggerName.Data());
 657                 }
 658         } else if (fIsHeavyIon == 1){
 659                 printf("Running in PbPb mode \n");
 660                 if (fDetectorCentrality == 0){
 661                         printf("\t centrality selection based on V0M \n");
 662                 } else if (fDetectorCentrality == 1){
 663                         printf("\t centrality selection based on Cl1 \n");
 664                 }
 665                 if (fModCentralityClass == 0){
 666                         printf("\t %d - %d \n", fCentralityMin*10, fCentralityMax*10);
 667                 } else if ( fModCentralityClass == 1){
 668                         printf("\t %d - %d \n", fCentralityMin*5, fCentralityMax*5);
 669                 } else if ( fModCentralityClass == 2){
 670                         printf("\t %d - %d \n", fCentralityMin*5+45, fCentralityMax*5+45);
 671                 } else if (fModCentralityClass == 3){
 672                         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*10, fCentralityMax*10);
 673                 } else if ( fModCentralityClass == 4){
 674                         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*5, fCentralityMax*5);
 675                 } else if ( fModCentralityClass == 5){
 676                         printf("\t %d - %d, with Track mult in MC as data \n", fCentralityMin*5+45, fCentralityMax*5+45);
 677                 }
 678                 if (fSpecialTrigger == 0){
 679                         printf("\t only events triggered by kMB, kCentral, kSemiCentral will be analysed \n");
 680                 } else if (fSpecialTrigger > 4){
 681                         printf("\t only events triggered by %s %s\n", fSpecialTriggerName.Data(), fSpecialSubTriggerName.Data());
 682                 }
 683         } else if (fIsHeavyIon == 2){
 684                 printf("Running in pPb mode \n");
 685                 if (fDetectorCentrality == 0){
 686                         printf("\t centrality selection based on V0A \n");
 687                 } else if (fDetectorCentrality == 1){
 688                         printf("\t centrality selection based on Cl1 \n");
 689                 }
 690                 if (fModCentralityClass == 0){
 691                         printf("\t %d - %d \n", fCentralityMin*10, fCentralityMax*10);
 692                 }
 693                 if (fSpecialTrigger == 0){
 694                         printf("\t only events triggered by kINT7 will be analysed \n");
 695                 } else if (fSpecialTrigger > 4){
 696                         printf("\t only events triggered by %s %s\n", fSpecialTriggerName.Data(), fSpecialSubTriggerName.Data());
 697                 }
 698         }
 699         printf("MC event cuts: \n");
 700         if (fRejectExtraSignals == 0) printf("\t no rejection was applied \n");
 701                 else if (fRejectExtraSignals == 1) printf("\t only MB header will be inspected \n");
 702                 else if (fRejectExtraSignals > 1) printf("\t special header have been selected \n");
 703 }
 704
 705 ///________________________________________________________________________
 706 Bool_t AliConvEventCuts::SetIsHeavyIon(Int_t isHeavyIon)
 707 {   // Set Cut
 708         switch(isHeavyIon){
 709         case 0:
 710                 fIsHeavyIon=0;
 711                 break;
 712         case 1:
 713                 fIsHeavyIon=1;
 714                 fDetectorCentrality=0;
 715                 break;
 716         case 2:
 717                 fIsHeavyIon=1;
 718                 fDetectorCentrality=1;
 719                 break;
 720         case 3: //allows to select centrality 0-45% in steps of 5% for V0 Multiplicity
 721                 fIsHeavyIon=1;
 722                 fDetectorCentrality=0;
 723                 fModCentralityClass=1;
 724                 break;
 725         case 4: //allows to select centrality 45-90% in steps of 5% for V0 Multiplicity
 726                 fIsHeavyIon=1;
 727                 fDetectorCentrality=0;
 728                 fModCentralityClass=2;
 729                 break;
 730         case 5: //strict cut on v0 tracks for MC
 731                 fIsHeavyIon=1;
 732                 fDetectorCentrality=0;
 733                 fModCentralityClass=3;
 734                 break;
 735         case 6: //allows to select centrality 0-45% in steps of 5% for track mult
 736                 //strict cut on v0 tracks for MC
 737                 fIsHeavyIon=1;
 738                 fDetectorCentrality=0;
 739                 fModCentralityClass=4;
 740                 break;
 741         case 7: //allows to select centrality 45-90% in steps of 5% for V0 Multiplicity
 742                 //strict cut on v0 tracks for MC
 743                 fIsHeavyIon=1;
 744                 fDetectorCentrality=0;
 745                 fModCentralityClass=5;
 746                 break;
 747         case 8:
 748                 fIsHeavyIon=2;
 749                 fDetectorCentrality=0;
 750                 break;
 751         case 9:
 752                 fIsHeavyIon=2;
 753                 fDetectorCentrality=1;
 754                 break;
 755         default:
 756                 AliError(Form("SetHeavyIon not defined %d",isHeavyIon));
 757                 return kFALSE;
 758         }
 759         return kTRUE;
 760 }
 761
 762 //___________________________________________________________________
 763 Bool_t AliConvEventCuts::SetCentralityMin(Int_t minCentrality)
 764 {
 765         // Set Cut
 766         if(minCentrality<0||minCentrality>9){
 767                 AliError(Form("minCentrality not defined %d",minCentrality));
 768                 return kFALSE;
 769         }
 770
 771         fCentralityMin=minCentrality;
 772         return kTRUE;
 773 }
 774
 775 //___________________________________________________________________
 776 Bool_t AliConvEventCuts::SetCentralityMax(Int_t maxCentrality)
 777 {
 778         // Set Cut
 779         if(maxCentrality<0||maxCentrality>9){
 780                 AliError(Form("maxCentrality not defined %d",maxCentrality));
 781                 return kFALSE;
 782         }
 783         fCentralityMax=maxCentrality;
 784         return kTRUE;
 785 }
 786
 787 ///________________________________________________________________________
 788 Bool_t AliConvEventCuts::SetSelectSpecialTrigger(Int_t selectSpecialTrigger)
 789 {// Set Cut
 790
 791         switch(selectSpecialTrigger){
 792         case 0:
 793                 fSpecialTrigger=0; // dont care
 794                 break;
 795         case 1:
 796                 fSpecialTrigger=1; // V0AND
 797                 break;
 798         case 2:
 799                 fSpecialTrigger=2; // with SDD requested
 800                 break;
 801         case 3:
 802                 fSpecialTrigger=3; // V0AND plus with SDD requested
 803                 break;
 804         // allows to run MB & 6 other different trigger classes in parallel with the same photon cut
 805         case 4:
 806                 fSpecialTrigger=4; // trigger alias kTRD
 807                 fOfflineTriggerMask=AliVEvent::kTRD;
 808                 fTriggerSelectedManually = kTRUE;
 809                 fSpecialTriggerName="AliVEvent::kTRD";
 810                 break;
 811         case 5:
 812                 fSpecialTrigger=5; // trigger alias kEMC
 813                 fOfflineTriggerMask=AliVEvent::kEMC7 | AliVEvent::kEMC8 | AliVEvent::kEMC1 ;
 814                 fTriggerSelectedManually = kTRUE;
 815                 fSpecialTriggerName="AliVEvent::kEMC7/kEMC8/kEMC1";
 816                 break;
 817         case 6:
 818                 fSpecialTrigger=6; // trigger alias kPHI
 819                 fOfflineTriggerMask=AliVEvent::kPHI7 | AliVEvent::kPHI1 | AliVEvent::kPHI8 | AliVEvent::kPHOSPb;
 820                 fTriggerSelectedManually = kTRUE;
 821                 fSpecialTriggerName="AliVEvent::kPHI7/kPHI1/kPHI8/kPHOSPb";
 822                 break;
 823         case 7:
 824                 fSpecialTrigger=7; // trigger alias kHighMult
 825                 fOfflineTriggerMask=AliVEvent::kHighMult;
 826                 fTriggerSelectedManually = kTRUE;
 827                 fSpecialTriggerName="AliVEvent::kHighMult";
 828                 break;
 829                 case 8:
 830                 fSpecialTrigger=8; // trigger alias kEMCEGA
 831                 fOfflineTriggerMask=AliVEvent::kEMCEGA;
 832                 fTriggerSelectedManually = kTRUE;
 833                 fSpecialTriggerName="AliVEvent::kEMCEGA";
 834                 break;
 835                 case 9:
 836                 fSpecialTrigger=9; // trigger alias kEMCEJE
 837                 fOfflineTriggerMask=AliVEvent::kEMCEJE;
 838                 fTriggerSelectedManually = kTRUE;
 839                 fSpecialTriggerName="AliVEvent::kEMCEJE";
 840                 break;
 841         default:
 842                 AliError("Warning: Special Trigger Not known");
 843                 return 0;
 844         }
 845         return 1;
 846 }
 847
 848 ///________________________________________________________________________
 849 Bool_t AliConvEventCuts::SetSelectSubTriggerClass(Int_t selectSpecialSubTriggerClass)
 850 {// Set Cut
 851
 852         if (fSpecialTrigger == 1){ //V0AND with different detectors
 853                 switch(selectSpecialSubTriggerClass){
 854                 case 0: //with VZERO
 855                         fSpecialTrigger=1;
 856                         fSpecialSubTrigger=0;
 857 //                      AliInfo("Info: Nothing to be done");
 858                         break;
 859                 case 1: //with TZERO
 860                         fSpecialTrigger=0;
 861                         fSpecialSubTrigger=0;
 862                         fOfflineTriggerMask=AliVEvent::kINT8;
 863                         fTriggerSelectedManually = kTRUE;
 864                         fSpecialTriggerName="AliVEvent::kINT8";
 865                         break;
 866                 default:
 867                         AliError("Warning: Special Subtrigger Class Not known");
 868                         return 0;
 869                 }
 870
 871         } else if (fSpecialTrigger == 4){ // Subdivision of TRD trigger classes
 872                 switch(selectSpecialSubTriggerClass){
 873                 case 0: // all together
 874                         fSpecialSubTrigger=0;
 875                         fSpecialSubTriggerName="";
 876 //                      AliInfo("Info: Nothing to be done");
 877                         break;
 878                 case 1: // 7WUHSH - V0AND with single electron in TRD & EMCAL
 879                         fSpecialSubTrigger=1;
 880                         fNSpecialSubTriggerOptions=1;
 881                         fSpecialSubTriggerName="7WUHEE";
 882                         break;
 883                 case 2: // 8WUHSH - T0AND with single electron in TRD & EMCAL
 884                         fSpecialSubTrigger=1;
 885                         fNSpecialSubTriggerOptions=1;
 886                         fSpecialSubTriggerName="8WUHEE";
 887                         break;
 888                 case 3: // 7WUHSE - V0AND with single high pt electron in TRD
 889                         fSpecialSubTrigger=1;
 890                         fNSpecialSubTriggerOptions=1;
 891                         fSpecialSubTriggerName="7WUHSE";
 892                         break;
 893                 case 4: // 8WUHSE - T0AND with single high pt electron in TRD
 894                         fSpecialSubTrigger=1;
 895                         fNSpecialSubTriggerOptions=1;
 896                         fSpecialSubTriggerName="8WUHSE";
 897                         break;
 898                 case 5: // 7WUHJE - V0AND with jet in TRD
 899                         fSpecialSubTrigger=1;
 900                         fNSpecialSubTriggerOptions=1;
 901                         fSpecialSubTriggerName="7WUHJT";
 902                         break;
 903                 case 6: // 8WUHJE - T0AND with jet in TRD
 904                         fSpecialSubTrigger=1;
 905                         fNSpecialSubTriggerOptions=1;
 906                         fSpecialSubTriggerName="8WUHJT";
 907                         break;
 908                 case 7: // 7WUHQU - V0AND with dielectron pair in TRD
 909                         fSpecialSubTrigger=1;
 910                         fNSpecialSubTriggerOptions=1;
 911                         fSpecialSubTriggerName="7WUHQU";
 912                         break;
 913                 case 8: // 8WUHQU - T0AND with dielectron pair in TRD
 914                         fSpecialSubTrigger=1;
 915                         fNSpecialSubTriggerOptions=1;
 916                         fSpecialSubTriggerName="8WUHQU";
 917                         break;
 918                 default:
 919                         AliError("Warning: Special Subtrigger Class Not known");
 920                         return 0;
 921                 }
 922         } else if (fSpecialTrigger == 5){ // Subdivision of kEMC trigger classes
 923                 switch(selectSpecialSubTriggerClass){
 924                 case 0: // all together
 925                         fSpecialSubTrigger=0;
 926                         fSpecialSubTriggerName="";
 927 //                      AliInfo("Info: Nothing to be done");
 928                         break;
 929                 case 1: // CEMC1 - V0OR and EMCAL fired
 930                         fOfflineTriggerMask=AliVEvent::kEMC1;
 931                         fSpecialTriggerName="AliVEvent::kEMC1";
 932                         fSpecialSubTrigger=1;
 933                         fNSpecialSubTriggerOptions=1;
 934                         fSpecialSubTriggerName="CEMC1";
 935                         break;
 936                 case 2: // CEMC7 - V0AND and EMCAL fired
 937                         fSpecialSubTrigger=1;
 938                         fOfflineTriggerMask=AliVEvent::kEMC7;
 939                         fSpecialTriggerName="AliVEvent::kEMC7";
 940                         fNSpecialSubTriggerOptions=1;
 941                         fSpecialSubTriggerName="CEMC7";
 942                         break;
 943                 case 3: // CEMC8  - T0OR and EMCAL fired
 944                         fOfflineTriggerMask=AliVEvent::kEMC8;
 945                         fSpecialTriggerName="AliVEvent::kEMC8";
 946                         fSpecialSubTrigger=1;
 947                         fNSpecialSubTriggerOptions=1;
 948                         fSpecialSubTriggerName="CEMC8";
 949                         break;
 950                 default:
 951                         AliError("Warning: Special Subtrigger Class Not known");
 952                         return 0;
 953                 }
 954         }  else if (fSpecialTrigger == 6){ // Subdivision of kPHI trigger classes
 955                 switch(selectSpecialSubTriggerClass){
 956                 case 0: // all together
 957                         fSpecialSubTrigger=0;
 958                         fSpecialSubTriggerName="";
 959 //                      AliInfo("Info: Nothing to be done");
 960                         break;
 961                 case 1: // CEMC1 - V0OR and EMCAL fired
 962                         fOfflineTriggerMask=AliVEvent::kPHI1;
 963                         fSpecialTriggerName="AliVEvent::kPHI1";
 964                         fSpecialSubTrigger=1;
 965                         fNSpecialSubTriggerOptions=1;
 966                         fSpecialSubTriggerName="CPHI1";
 967                         break;
 968                 case 2: // CEMC7 - V0AND and EMCAL fired
 969                         fSpecialSubTrigger=1;
 970                         fOfflineTriggerMask=AliVEvent::kPHI7;
 971                         fSpecialTriggerName="AliVEvent::kPHI7";
 972                         fNSpecialSubTriggerOptions=1;
 973                         fSpecialSubTriggerName="CPHI7";
 974                         break;
 975                 case 3: // CEMC8  - T0OR and EMCAL fired
 976                         fOfflineTriggerMask=AliVEvent::kPHI8;
 977                         fSpecialTriggerName="AliVEvent::kPHI8";
 978                         fSpecialSubTrigger=1;
 979                         fNSpecialSubTriggerOptions=1;
 980                         fSpecialSubTriggerName="CPHI8";
 981                         break;
 982                 default:
 983                         AliError("Warning: Special Subtrigger Class Not known");
 984                         return 0;
 985                 }
 986         } else if (fSpecialTrigger == 7){ // Subdivision of kHighMult trigger classes
 987                 switch(selectSpecialSubTriggerClass){
 988                 case 0: // all together
 989                         fSpecialSubTrigger=0;
 990                         fSpecialSubTriggerName="";
 991 //                      AliInfo("Info: Nothing to be done");
 992                         break;
 993                 case 1: // CSHM1 - V0OR and high mult fired
 994                         fSpecialSubTrigger=1;
 995                         fNSpecialSubTriggerOptions=1;
 996                         fSpecialSubTriggerName="CSHM1";
 997                         break;
 998                 case 2: // CSHM7 - V0AND and high mult fired
 999                         fSpecialSubTrigger=1;
1000                         fNSpecialSubTriggerOptions=1;
1001                         fSpecialSubTriggerName="CSHM7";
1002                         break;
1003                 case 3: // CSHM8  - T0OR and high mult fired
1004                         fSpecialSubTrigger=1;
1005                         fNSpecialSubTriggerOptions=1;
1006                         fSpecialSubTriggerName="CSHM8";
1007                         break;
1008                 default:
1009                         AliError("Warning: Special Subtrigger Class Not known");
1010                         return 0;
1011                 }
1012         }  else if (fSpecialTrigger == 8){ // Subdivision of kEMCEGA trigger classes
1013                 switch(selectSpecialSubTriggerClass){
1014                 case 0: // all together
1015                         fSpecialSubTrigger=0;
1016                         fSpecialSubTriggerName="";
1017 //                      AliInfo("Info: Nothing to be done");
1018                         break;
1019                 case 1: // 7EGA - CINT7 EGA
1020                         fSpecialSubTrigger=1;
1021                         fNSpecialSubTriggerOptions=1;
1022                         fSpecialSubTriggerName="7EGA";
1023                         break;
1024                 case 2: // 8EGA - CINT8 EGA
1025                         fSpecialSubTrigger=1;
1026                         fNSpecialSubTriggerOptions=1;
1027                         fSpecialSubTriggerName="8EGA";
1028                         break;
1029                 case 3: // 7EG1 - CINT7 EG1
1030                         fSpecialSubTrigger=1;
1031                         fNSpecialSubTriggerOptions=1;
1032                         fSpecialSubTriggerName="7EG1";
1033                         break;
1034                 case 4: // 8EG1 - CINT8 EG1
1035                         fSpecialSubTrigger=1;
1036                         fNSpecialSubTriggerOptions=1;
1037                         fSpecialSubTriggerName="8EG1";
1038                         break;
1039                 case 5: // 7EG2 - CINT7 EG2
1040                         fSpecialSubTrigger=1;
1041                         fNSpecialSubTriggerOptions=1;
1042                         fSpecialSubTriggerName="7EG2";
1043                         break;
1044                 case 6: // 8EG2 - CINT8 EG2
1045                         fSpecialSubTrigger=1;
1046                         fNSpecialSubTriggerOptions=1;
1047                         fSpecialSubTriggerName="8EG2";
1048                         break;
1049                 default:
1050                         AliError("Warning: Special Subtrigger Class Not known");
1051                         return 0;
1052                 }
1053         } else if (fSpecialTrigger == 9){ // Subdivision of kEMCEGA trigger classes
1054                 switch(selectSpecialSubTriggerClass){
1055                 case 0: // all together
1056                         fSpecialSubTrigger=0;
1057                         fSpecialSubTriggerName="";
1058 //                      AliInfo("Info: Nothing to be done");
1059                         break;
1060                 case 1: // 7EJE - CINT7 EJE
1061                         fSpecialSubTrigger=1;
1062                         fNSpecialSubTriggerOptions=1;
1063                         fSpecialSubTriggerName="7EJE";
1064                         break;
1065                 case 2: // 8EJE - CINT8 EJE
1066                         fSpecialSubTrigger=1;
1067                         fNSpecialSubTriggerOptions=1;
1068                         fSpecialSubTriggerName="8EJE";
1069                         break;
1070                 case 3: // 7EJ1 - CINT7 EJ1
1071                         fSpecialSubTrigger=1;
1072                         fNSpecialSubTriggerOptions=1;
1073                         fSpecialSubTriggerName="7EJ1";
1074                         break;
1075                 case 4: // 8EJ1 - CINT8 EJ1
1076                         fSpecialSubTrigger=1;
1077                         fNSpecialSubTriggerOptions=1;
1078                         fSpecialSubTriggerName="8EJ1";
1079                         break;
1080                 case 5: // 7EJ2 - CINT7 EJ2
1081                         fSpecialSubTrigger=1;
1082                         fNSpecialSubTriggerOptions=1;
1083                         fSpecialSubTriggerName="7EJ2";
1084                         break;
1085                 case 6: // 8EJ2 - CINT8 EJ2
1086                         fSpecialSubTrigger=1;
1087                         fNSpecialSubTriggerOptions=1;
1088                         fSpecialSubTriggerName="8EJ2";
1089                         break;
1090                 default:
1091                         AliError("Warning: Special Subtrigger Class Not known");
1092                         return 0;
1093                 }
1094         }
1095         return 1;
1096 }
1097
1098 ///________________________________________________________________________
1099 Bool_t AliConvEventCuts::SetMultiplicityMethod(Int_t multiplicityMethod)
1100 {
1101         // Set Cut
1102         fMultiplicityMethod=multiplicityMethod;
1103
1104         // 0 Photon Multiplicity
1105         // 1 TPC Track multiplicity
1106         // 2 V0 Mult
1107         // 3 SPD Mult
1108
1109         return kTRUE;
1110 }
1111
1112 ///________________________________________________________________________
1113 Bool_t AliConvEventCuts::SetRemovePileUp(Int_t removePileUp)
1114 {// Set Cut
1115         switch(removePileUp){
1116         case 0:
1117                 fRemovePileUp=kFALSE;
1118                 break;
1119         case 1:
1120                 fRemovePileUp=kTRUE;
1121                 break;
1122         default:
1123                 AliError("RemovePileUpCut not defined");
1124                 return kFALSE;
1125         }
1126         return kTRUE;
1127 }
1128
1129 ///________________________________________________________________________
1130 Bool_t AliConvEventCuts::SetRejectExtraSignalsCut(Int_t extraSignal) {
1131
1132         switch(extraSignal){
1133         case 0:
1134                 fRejectExtraSignals = 0;
1135                 break; // No Rejection
1136         case 1:
1137                 fRejectExtraSignals = 1;
1138                 break; // MinBias Header
1139         case 2:
1140                 fRejectExtraSignals = 2;
1141                 break; // User String Array
1142         case 3:
1143                 fRejectExtraSignals = 3;
1144                 break; // Rejection for Gamma Correction only
1145         default:
1146                 AliError(Form("Extra Signal Rejection not defined %d",extraSignal));
1147                 return kFALSE;
1148         }
1149         return kTRUE;
1150 }
1151
1152 //-------------------------------------------------------------
1153 Double_t AliConvEventCuts::GetCentrality(AliVEvent *event)
1154 {   // Get Event Centrality
1155
1156         AliESDEvent *esdEvent=dynamic_cast<AliESDEvent*>(event);
1157         if(esdEvent){
1158                 AliCentrality *fESDCentrality=(AliCentrality*)esdEvent->GetCentrality();
1159                 if(fDetectorCentrality==0){
1160                         if (fIsHeavyIon==2){
1161                                 return fESDCentrality->GetCentralityPercentile("V0A"); // default for pPb
1162                         } else{
1163                                 return fESDCentrality->GetCentralityPercentile("V0M"); // default
1164                         }
1165                 }
1166                 if(fDetectorCentrality==1){
1167                         return fESDCentrality->GetCentralityPercentile("CL1");
1168                 }
1169         }
1170
1171         AliAODEvent *aodEvent=dynamic_cast<AliAODEvent*>(event);
1172         if(aodEvent){
1173                 if(aodEvent->GetHeader()){return aodEvent->GetHeader()->GetCentrality();}
1174         }
1175
1176         return -1;
1177 }
1178
1179 //_____________________________________________________________________________________
1180 Bool_t AliConvEventCuts::IsCentralitySelected(AliVEvent *event, AliVEvent *fMCEvent)
1181 {   // Centrality Selection
1182         if(!fIsHeavyIon)return kTRUE;
1183
1184         if(fCentralityMin == fCentralityMax ) return kTRUE;//0-100%
1185         else if(fCentralityMax==0) fCentralityMax=10; //CentralityRange = fCentralityMin-100%
1186
1187         Double_t centrality=GetCentrality(event);
1188         if(centrality<0)return kFALSE;
1189
1190         Int_t centralityC=0;
1191         if (fModCentralityClass == 0){
1192                 centralityC= Int_t(centrality/10);
1193                 if(centralityC >= fCentralityMin && centralityC < fCentralityMax)
1194                         return kTRUE;
1195                 else return kFALSE;
1196         }
1197         else if (fModCentralityClass ==1){
1198                 centralityC= Int_t(centrality);
1199                 if(centralityC >= fCentralityMin*5 && centralityC < fCentralityMax*5){
1200                         return kTRUE;
1201                 } else return kFALSE;
1202         }
1203         else if (fModCentralityClass ==2){
1204                 centralityC= Int_t(centrality);
1205                 if(centralityC >= ((fCentralityMin*5)+45) && centralityC < ((fCentralityMax*5)+45))
1206                         return kTRUE;
1207                 else return kFALSE;
1208         }
1209
1210         Int_t nprimaryTracks = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetNumberOfPrimaryTracks();
1211         Int_t PrimaryTracks10[10][2] =
1212                 {
1213                         {9999,9999}, //  0
1214                         {1210, 928}, // 10
1215                         { 817, 658}, // 20
1216                         { 536, 435}, // 30
1217                         { 337, 276}, // 40
1218                         { 197, 162}, // 50
1219                         { 106, 100}, // 60
1220                         {  51,  44}, // 70
1221                         {  21,  18}, // 80
1222                         {   0,   0}  // 90
1223                 };
1224         Int_t PrimaryTracks5a[10][2] =
1225                 {
1226                         {9999,9999}, // 0
1227                         {1485,1168}, // 5
1228                         {1210, 928}, // 10
1229                         { 995, 795}, // 15
1230                         { 817, 658}, // 20
1231                         { 666, 538}, // 25
1232                         { 536, 435}, // 30
1233                         { 428, 350}, // 35
1234                         { 337, 276}, // 40
1235                         { 260, 214}  // 45
1236                 };
1237         Int_t PrimaryTracks5b[10][2] =
1238                 {
1239                         { 260, 214}, // 45
1240                         { 197, 162}, // 50
1241                         { 147, 125}, // 55
1242                         { 106, 100}, // 60
1243                         {  75,  63}, // 65
1244                         {  51,  44}, // 70
1245                         {  34,  29}, // 75
1246                         {  21,  18}, // 80
1247                         {  13,  11}, // 85
1248                         {   0,   0}  // 90
1249                 };
1250         Int_t column = 0;
1251         if(event->IsA()==AliESDEvent::Class()) column = 0;
1252         if(event->IsA()==AliAODEvent::Class()) column = 1;
1253
1254         if (fModCentralityClass == 3){
1255                 if(fMCEvent){
1256                         if(nprimaryTracks > PrimaryTracks10[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks10[fCentralityMin][column])
1257                                 return kTRUE;
1258                         else return kFALSE;
1259                 }
1260                 else{
1261                         centralityC= Int_t(centrality/10);
1262                         if(centralityC >= fCentralityMin && centralityC < fCentralityMax)
1263                                 return kTRUE;
1264                         else return kFALSE;
1265                 }
1266         }
1267         else if (fModCentralityClass ==4){
1268                 if(fMCEvent){
1269                         if(nprimaryTracks > PrimaryTracks5a[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks5a[fCentralityMin][column])
1270                                 return kTRUE;
1271                         else return kFALSE;
1272                 }
1273                 else{
1274                         centralityC= Int_t(centrality);
1275                         if(centralityC >= fCentralityMin*5 && centralityC < fCentralityMax*5){
1276                                 return kTRUE;
1277                         } else return kFALSE;
1278                 }
1279         }
1280         else if (fModCentralityClass ==5){
1281                 if(fMCEvent){
1282                         if(nprimaryTracks > PrimaryTracks5b[fCentralityMax][column] && nprimaryTracks <= PrimaryTracks5b[fCentralityMin][column])
1283                                 return kTRUE;
1284                         else return kFALSE;
1285                 }
1286                 else{
1287                         centralityC= Int_t(centrality);
1288                         if(centralityC >= ((fCentralityMin*5)+45) && centralityC < ((fCentralityMax*5)+45))
1289                                 return kTRUE;
1290                         else return kFALSE;
1291                 }
1292         }
1293
1294         return kFALSE;
1295 }
1296
1297 ///________________________________________________________________________
1298 Bool_t AliConvEventCuts::VertexZCut(AliVEvent *event){
1299         // Cut on z position of primary vertex
1300         Double_t fVertexZ=event->GetPrimaryVertex()->GetZ();
1301         Double_t fVertexZSPD = 0;
1302         AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(event);
1303         if(fESDEvent){
1304                 fVertexZSPD = fESDEvent->GetPrimaryVertexSPD()->GetZ();
1305         }
1306         AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(event);
1307         if(fAODEvent){
1308                 fVertexZSPD = fAODEvent->GetPrimaryVertexSPD()->GetZ();
1309         }
1310
1311         if(abs(fVertexZ)>fMaxVertexZ)return kFALSE;
1312
1313         TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
1314                                                                                                         ->GetTask(fV0ReaderName.Data()))->GetPeriodName();
1315         if (periodName.CompareTo("LHC11h")==0){
1316                 if (abs(fVertexZ-fVertexZSPD) > 0.1) return kFALSE;
1317         }
1318         if (fIsHeavyIon == 2){
1319                 if(fUtils->IsFirstEventInChunk(event)) return kFALSE;
1320                 if(!fUtils->IsVertexSelected2013pA(event)) return kFALSE;
1321                 if(fUtils->IsPileUpEvent(event)) return kFALSE;
1322         }
1323
1324         return kTRUE;
1325 }
1326
1327 ///________________________________________________________________________
1328 Int_t AliConvEventCuts::GetNumberOfContributorsVtx(AliVEvent *event){
1329         // returns number of contributors to the vertex
1330
1331         AliESDEvent *fESDEvent=dynamic_cast<AliESDEvent*>(event);
1332         if(fESDEvent){
1333                 if (fESDEvent->GetPrimaryVertex() != NULL){
1334                         if(fESDEvent->GetPrimaryVertex()->GetNContributors()>0) {
1335         //     cout << "accepted global" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertex()->GetNContributors() << endl;
1336                                 return fESDEvent->GetPrimaryVertex()->GetNContributors();
1337                         }
1338                 }
1339
1340                 if(fESDEvent->GetPrimaryVertexSPD() !=NULL){
1341                         if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
1342         //     cout << "accepted SPD" << fESDEvent->GetEventNumberInFile() << " with NCont: " << fESDEvent->GetPrimaryVertexSPD()->GetNContributors() << endl;
1343                                 return fESDEvent->GetPrimaryVertexSPD()->GetNContributors();
1344                         }  else {
1345                                 AliWarning(Form("Number of contributors from bad vertex type:: %s",fESDEvent->GetPrimaryVertex()->GetName()));
1346         //            cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
1347                                 return 0;
1348                         }
1349                 }
1350         }
1351
1352         AliAODEvent *fAODEvent=dynamic_cast<AliAODEvent*>(event);
1353         if(fAODEvent){
1354                 if (fAODEvent->GetPrimaryVertex() != NULL){
1355                         if(fAODEvent->GetPrimaryVertex()->GetNContributors()>0) {
1356                                 return fAODEvent->GetPrimaryVertex()->GetNContributors();
1357                         }
1358                 }
1359                 if(fAODEvent->GetPrimaryVertexSPD() !=NULL){
1360                         if(fAODEvent->GetPrimaryVertexSPD()->GetNContributors()>0) {
1361                                 return fAODEvent->GetPrimaryVertexSPD()->GetNContributors();
1362                         } else {
1363                                 AliWarning(Form("Number of contributors from bad vertex type:: %s",fAODEvent->GetPrimaryVertex()->GetName()));
1364                                 return 0;
1365                         }
1366                 }
1367         }
1368         // cout << "rejected " << fESDEvent->GetEventNumberInFile() << endl;
1369         return 0;
1370 }
1371
1372
1373 ///________________________________________________________________________
1374 Bool_t AliConvEventCuts::IsTriggerSelected(AliVEvent *fInputEvent)
1375 {
1376
1377         AliInputEventHandler *fInputHandler=(AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
1378
1379         UInt_t isSelected = AliVEvent::kAny;
1380         TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
1381         //    cout <<   periodName.Data() << endl;
1382
1383         if (fInputHandler==NULL) return kFALSE;
1384         if( fInputHandler->GetEventSelection() || fInputEvent->IsA()==AliAODEvent::Class()) {
1385
1386                 TString firedTrigClass = fInputEvent->GetFiredTriggerClasses();
1387                 if (!fTriggerSelectedManually){
1388                         if (fPreSelCut) fOfflineTriggerMask = AliVEvent::kAny;
1389                         else {
1390                                 if (fIsHeavyIon == 1) fOfflineTriggerMask = AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral;
1391                                 else if (fIsHeavyIon == 2) fOfflineTriggerMask = AliVEvent::kINT7;
1392                                 else if (periodName.CompareTo("LHC11c") == 0 || periodName.CompareTo("LHC11d") == 0 || periodName.CompareTo("LHC11e") == 0 || periodName.CompareTo("LHC11f") == 0 || periodName.CompareTo("LHC11g") == 0  || periodName.CompareTo("LHC12a") == 0 || periodName.CompareTo("LHC12b") == 0 || periodName.CompareTo("LHC12c") == 0 || periodName.CompareTo("LHC12d") == 0 || periodName.CompareTo("LHC12f") == 0  || periodName.CompareTo("LHC12g") == 0  || periodName.CompareTo("LHC12h") == 0  || periodName.CompareTo("LHC12i") == 0  ||periodName.CompareTo("LHC13g") == 0 ) {
1393                                         fOfflineTriggerMask = AliVEvent::kINT7;
1394         //                              cout << "will take kINT7 as trigger mask" << endl;
1395                                 }
1396                                 else fOfflineTriggerMask = AliVEvent::kMB;
1397                         }
1398                 }
1399                 // Get the actual offline trigger mask for the event and AND it with the
1400                 // requested mask. If no mask requested select by default the event.
1401         //       if (fPreSelCut) cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask for Precut: " << fOfflineTriggerMask << endl;
1402         //       else cout << "Trigger selected from outside: "<< fTriggerSelectedManually <<"\t Offline Trigger mask: " << fOfflineTriggerMask << endl;
1403
1404                 if (fOfflineTriggerMask){
1405                         isSelected = fOfflineTriggerMask & fInputHandler->IsEventSelected();
1406                         if (isSelected && !fPreSelCut){
1407                                 if (fSpecialSubTrigger>0){
1408                                         if (!firedTrigClass.Contains(fSpecialSubTriggerName.Data())) isSelected = 0;
1409                                 }
1410                         }
1411                 }
1412         }
1413         fIsSDDFired = !(fInputHandler->IsEventSelected() & AliVEvent::kFastOnly);
1414
1415         // Fill Histogram
1416         if(hTriggerClass){
1417                 if (fIsSDDFired) hTriggerClass->Fill(33);
1418                 if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClass->Fill(0);
1419                 if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClass->Fill(1);
1420                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClass->Fill(2);
1421                 if (fInputHandler->IsEventSelected() & AliVEvent::kHighMult)hTriggerClass->Fill(3);
1422                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMC1)hTriggerClass->Fill(4);
1423                 if (fInputHandler->IsEventSelected() & AliVEvent::kCINT5)hTriggerClass->Fill(5);
1424                 if (fInputHandler->IsEventSelected() & AliVEvent::kCMUS5)hTriggerClass->Fill(6);
1425         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSPB)hTriggerClass->Fill(6);
1426                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUSH7)hTriggerClass->Fill(7);
1427         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSHPB)hTriggerClass->Fill(7);
1428                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUL7)hTriggerClass->Fill(8);
1429         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikePB)hTriggerClass->Fill(8);
1430                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUU7)hTriggerClass->Fill(9);
1431         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikePB)hTriggerClass->Fill(9);
1432                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMC7)hTriggerClass->Fill(10);
1433         //       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC8)hTriggerClass->Fill(10);
1434                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUS7)hTriggerClass->Fill(11);
1435                 if (fInputHandler->IsEventSelected() & AliVEvent::kPHI1)hTriggerClass->Fill(12);
1436                 if (fInputHandler->IsEventSelected() & AliVEvent::kPHI7)hTriggerClass->Fill(13);
1437         //       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI8)hTriggerClass->Fill(13);
1438         //       if (fInputHandler->IsEventSelected() & AliVEvent::kPHOSPb)hTriggerClass->Fill(13);
1439                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)hTriggerClass->Fill(14);
1440                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA)hTriggerClass->Fill(15);
1441                 if (fInputHandler->IsEventSelected() & AliVEvent::kCentral)hTriggerClass->Fill(16);
1442                 if (fInputHandler->IsEventSelected() & AliVEvent::kSemiCentral)hTriggerClass->Fill(17);
1443                 if (fInputHandler->IsEventSelected() & AliVEvent::kDG5)hTriggerClass->Fill(18);
1444                 if (fInputHandler->IsEventSelected() & AliVEvent::kZED)hTriggerClass->Fill(19);
1445                 if (fInputHandler->IsEventSelected() & AliVEvent::kSPI7)hTriggerClass->Fill(20);
1446         //       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI)hTriggerClass->Fill(20);
1447                 if (fInputHandler->IsEventSelected() & AliVEvent::kINT8)hTriggerClass->Fill(21);
1448                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleLowPt8)hTriggerClass->Fill(22);
1449                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleHighPt8)hTriggerClass->Fill(23);
1450                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikeLowPt8)hTriggerClass->Fill(24);
1451                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClass->Fill(25);
1452                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClass->Fill(26);
1453                 if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClass->Fill(27);
1454                 if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClass->Fill(28);
1455                 if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClass->Fill(29);
1456                 if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClass->Fill(30);
1457                 if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClass->Fill(31);
1458                 if (!fInputHandler->IsEventSelected()) hTriggerClass->Fill(34);
1459         }
1460
1461         if(hTriggerClassSelected && isSelected){
1462                 if (!fIsSDDFired) hTriggerClassSelected->Fill(33);
1463                 if (fInputHandler->IsEventSelected() & AliVEvent::kMB)hTriggerClassSelected->Fill(0);
1464                 if (fInputHandler->IsEventSelected() & AliVEvent::kINT7)hTriggerClassSelected->Fill(1);
1465                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUON)hTriggerClassSelected->Fill(2);
1466                 if (fInputHandler->IsEventSelected() & AliVEvent::kHighMult)hTriggerClassSelected->Fill(3);
1467                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMC1)hTriggerClassSelected->Fill(4);
1468                 if (fInputHandler->IsEventSelected() & AliVEvent::kCINT5)hTriggerClassSelected->Fill(5);
1469                 if (fInputHandler->IsEventSelected() & AliVEvent::kCMUS5)hTriggerClassSelected->Fill(6);
1470         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSPB)hTriggerClassSelected->Fill(6);
1471                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUSH7)hTriggerClassSelected->Fill(7);
1472         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMUSHPB)hTriggerClassSelected->Fill(7);
1473                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUL7)hTriggerClassSelected->Fill(8);
1474         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikePB)hTriggerClassSelected->Fill(8);
1475                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUU7)hTriggerClassSelected->Fill(9);
1476         //       if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikePB)hTriggerClassSelected->Fill(9);
1477                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMC7)hTriggerClassSelected->Fill(10);
1478         //       if (fInputHandler->IsEventSelected() & AliVEvent::kEMC8)hTriggerClassSelected->Fill(10);
1479                 if (fInputHandler->IsEventSelected() & AliVEvent::kMUS7)hTriggerClassSelected->Fill(11);
1480                 if (fInputHandler->IsEventSelected() & AliVEvent::kPHI1)hTriggerClassSelected->Fill(12);
1481                 if (fInputHandler->IsEventSelected() & AliVEvent::kPHI7)hTriggerClassSelected->Fill(13);
1482         //       if (fInputHandler->IsEventSelected() & AliVEvent::kPHI8)hTriggerClassSelected->Fill(13);
1483         //       if (fInputHandler->IsEventSelected() & AliVEvent::kPHOSPb)hTriggerClassSelected->Fill(13);
1484                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE)hTriggerClassSelected->Fill(14);
1485                 if (fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA)hTriggerClassSelected->Fill(15);
1486                 if (fInputHandler->IsEventSelected() & AliVEvent::kCentral)hTriggerClassSelected->Fill(16);
1487                 if (fInputHandler->IsEventSelected() & AliVEvent::kSemiCentral)hTriggerClassSelected->Fill(17);
1488                 if (fInputHandler->IsEventSelected() & AliVEvent::kDG5)hTriggerClassSelected->Fill(18);
1489                 if (fInputHandler->IsEventSelected() & AliVEvent::kZED)hTriggerClassSelected->Fill(19);
1490                 if (fInputHandler->IsEventSelected() & AliVEvent::kSPI7)hTriggerClassSelected->Fill(20);
1491         //       if (fInputHandler->IsEventSelected() & AliVEvent::kSPI)hTriggerClassSelected->Fill(20);
1492                 if (fInputHandler->IsEventSelected() & AliVEvent::kINT8)hTriggerClassSelected->Fill(21);
1493                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleLowPt8)hTriggerClassSelected->Fill(22);
1494                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonSingleHighPt8)hTriggerClassSelected->Fill(23);
1495                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonLikeLowPt8)hTriggerClassSelected->Fill(24);
1496                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt8)hTriggerClassSelected->Fill(25);
1497                 if (fInputHandler->IsEventSelected() & AliVEvent::kMuonUnlikeLowPt0)hTriggerClassSelected->Fill(26);
1498                 if (fInputHandler->IsEventSelected() & AliVEvent::kUserDefined)hTriggerClassSelected->Fill(27);
1499                 if (fInputHandler->IsEventSelected() & AliVEvent::kTRD)hTriggerClassSelected->Fill(28);
1500                 if (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly)hTriggerClassSelected->Fill(29);
1501                 if (fInputHandler->IsEventSelected() & AliVEvent::kAnyINT)hTriggerClassSelected->Fill(30);
1502                 if (fInputHandler->IsEventSelected() & AliVEvent::kAny)hTriggerClassSelected->Fill(31);
1503         }
1504
1505         if(!isSelected)return kFALSE;
1506         return kTRUE;
1507
1508 }
1509
1510 ///________________________________________________________________________
1511 TString AliConvEventCuts::GetCutNumber(){
1512    // returns TString with current cut number
1513    TString a(kNCuts);
1514    for(Int_t ii=0;ii<kNCuts;ii++){
1515       a.Append(Form("%d",fCuts[ii]));
1516    }
1517    return a;
1518 }
1519
1520 ///________________________________________________________________________
1521 void AliConvEventCuts::GetNotRejectedParticles(Int_t rejection, TList *HeaderList, AliVEvent *MCEvent){
1522
1523         TString periodName = ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data()))->GetPeriodName();
1524
1525
1526         if(fNotRejectedStart){
1527                 delete[] fNotRejectedStart;
1528                 fNotRejectedStart = NULL;
1529         }
1530         if(fNotRejectedEnd){
1531                 delete[] fNotRejectedEnd;
1532                 fNotRejectedEnd = NULL;
1533         }
1534         if(fGeneratorNames){
1535                 delete[] fGeneratorNames;
1536                 fGeneratorNames = NULL;
1537         }
1538
1539         if(rejection == 0) return; // No Rejection
1540
1541         AliGenCocktailEventHeader *cHeader = 0x0;
1542         AliAODMCHeader *cHeaderAOD = 0x0;
1543         Bool_t headerFound = kFALSE;
1544         AliStack *fMCStack = 0x0;
1545         TClonesArray *fMCStackAOD = 0x0;
1546         if(MCEvent->IsA()==AliMCEvent::Class()){
1547                 cHeader = dynamic_cast<AliGenCocktailEventHeader*>(dynamic_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
1548                 if(cHeader) headerFound = kTRUE;
1549                 if(dynamic_cast<AliMCEvent*>(MCEvent))fMCStack = dynamic_cast<AliStack*>(dynamic_cast<AliMCEvent*>(MCEvent)->Stack());
1550         }
1551         if(MCEvent->IsA()==AliAODEvent::Class()){ // MCEvent is a AODEvent in case of AOD
1552                 cHeaderAOD = dynamic_cast<AliAODMCHeader*>(MCEvent->FindListObject(AliAODMCHeader::StdBranchName()));
1553                 fMCStackAOD = dynamic_cast<TClonesArray*>(MCEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1554
1555
1556                 if(cHeaderAOD) headerFound = kTRUE;
1557         }
1558
1559         if(headerFound){
1560                 TList *genHeaders = 0x0;
1561                 if(cHeader) genHeaders = cHeader->GetHeaders();
1562                 if(cHeaderAOD){
1563                         genHeaders = cHeaderAOD->GetCocktailHeaders();
1564                         if(genHeaders->GetEntries()==1){
1565                                 SetRejectExtraSignalsCut(0);
1566                                 return;
1567                         }
1568                 }
1569                 AliGenEventHeader* gh = 0;
1570                 fnHeaders = 0;
1571                 Int_t firstindexA = 0;
1572                 Int_t lastindexA =  -1;
1573                 if(rejection == 1 || rejection == 3) fnHeaders = 1; // MinBiasHeader
1574                 if(rejection == 2){ // TList of Headers Names
1575                         for(Int_t i = 0; i<genHeaders->GetEntries();i++){
1576                                 gh = (AliGenEventHeader*)genHeaders->At(i);
1577                                 TString GeneratorName = gh->GetName();
1578                                 lastindexA = lastindexA + gh->NProduced();
1579 //                              cout << i << "\t" << GeneratorName.Data() << endl;
1580                                 for(Int_t j = 0; j<HeaderList->GetEntries();j++){
1581                                         TString GeneratorInList = ((TObjString*)HeaderList->At(j))->GetString();
1582 //                                      cout << GeneratorInList.Data() << endl;
1583                                         if(GeneratorName.CompareTo(GeneratorInList) == 0){
1584 //                                              cout << "accepted" << endl;
1585                                                 if (GeneratorInList.CompareTo("PARAM") == 0 || GeneratorInList.CompareTo("BOX") == 0 ){
1586                                                         if(fMCStack){
1587                                                                 if (fMCStack->Particle(firstindexA)->GetPdgCode() == fAddedSignalPDGCode ) {
1588                                                                         if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1c")==0 ){
1589                                                                                 if (gh->NProduced() > 10 && fMCStack->Particle(firstindexA+10)->GetPdgCode() == fAddedSignalPDGCode ){
1590 //                                                                                      cout << "cond 1: "<< fnHeaders << endl;
1591                                                                                         fnHeaders++;
1592                                                                                         continue;
1593                                                                                 }
1594                                                                                 continue;
1595                                                                         } else {
1596 //                                                                              cout << "cond 2: " << fnHeaders << endl;
1597                                                                                 fnHeaders++;
1598                                                                                 continue;
1599                                                                         }
1600                                                                 }
1601                                                         }
1602                                                         if ( fMCStackAOD){
1603                                                                 AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindexA));
1604                                                                 if (  aodMCParticle->GetPdgCode() == fAddedSignalPDGCode ){
1605                                                                         if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1c")==0 ){
1606                                                                                 if (gh->NProduced() > 10){
1607                                                                                         AliAODMCParticle *aodMCParticle2 = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindexA+10));
1608                                                                                         if (  aodMCParticle2->GetPdgCode() == fAddedSignalPDGCode ){
1609 //                                                                                              cout << "cond 1: " << fnHeaders << endl;
1610                                                                                                 fnHeaders++;
1611                                                                                                 continue;
1612                                                                                         }
1613                                                                                 }
1614                                                                                 continue;
1615                                                                         } else {
1616 //                                                                              cout << "cond 2: " << fnHeaders << endl;
1617                                                                                 fnHeaders++;
1618                                                                                 continue;
1619                                                                         }
1620                                                                 }
1621                                                         }
1622                                                         continue;
1623                                                 }
1624 //                                              cout << "cond 3: "<< fnHeaders << endl;
1625                                                 fnHeaders++;
1626                                                 continue;
1627                                         }
1628                                 }
1629                                 firstindexA = firstindexA + gh->NProduced();
1630                         }
1631                 }
1632 //              cout << "number of headers: " <<fnHeaders << endl;
1633
1634                 fNotRejectedStart = new Int_t[fnHeaders];
1635                 fNotRejectedEnd = new Int_t[fnHeaders];
1636                 fGeneratorNames = new TString[fnHeaders];
1637
1638                 if(rejection == 1 || rejection == 3){
1639                         fNotRejectedStart[0] = 0;
1640                         fNotRejectedEnd[0] = ((AliGenEventHeader*)genHeaders->At(0))->NProduced()-1;
1641                         fGeneratorNames[0] = ((AliGenEventHeader*)genHeaders->At(0))->GetName();
1642                         return;
1643                 }
1644
1645                 Int_t firstindex = 0;
1646                 Int_t lastindex =  -1;
1647                 Int_t number = 0;
1648
1649                 for(Int_t i = 0; i<genHeaders->GetEntries();i++){
1650                         gh = (AliGenEventHeader*)genHeaders->At(i);
1651                         TString GeneratorName = gh->GetName();
1652                         lastindex = lastindex + gh->NProduced();
1653                         for(Int_t j = 0; j<HeaderList->GetEntries();j++){
1654                                 TString GeneratorInList = ((TObjString*)HeaderList->At(j))->GetString();
1655 //                              cout << i << "\t" << GeneratorName.Data() << endl;
1656                                 if(GeneratorName.CompareTo(GeneratorInList) == 0){
1657                                         if (GeneratorInList.CompareTo("PARAM") == 0 || GeneratorInList.CompareTo("BOX") == 0 ){
1658                                                 if(fMCStack){
1659                                                         if (fMCStack->Particle(firstindex)->GetPdgCode() == fAddedSignalPDGCode ) {
1660                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1c")==0 ){
1661 //                                                                      cout << "produced " << gh->NProduced() << " with box generator" << endl;
1662                                                                         if (gh->NProduced() > 10 && fMCStack->Particle(firstindex+10)->GetPdgCode() == fAddedSignalPDGCode){
1663 //                                                                              cout << "one of them was a pi0 or eta" <<  endl;
1664                                                                                 fNotRejectedStart[number] = firstindex;
1665                                                                                 fNotRejectedEnd[number] = lastindex;
1666                                                                                 fGeneratorNames[number] = GeneratorName;
1667                                                                                 number++;
1668 //                                                                              cout << "Number of particles produced for: " << i << "\t" << GeneratorName.Data() << "\t" << lastindex-firstindex+1 << endl;
1669                                                                                 continue;
1670                                                                         }
1671                                                                 } else {
1672                                                                         fNotRejectedStart[number] = firstindex;
1673                                                                         fNotRejectedEnd[number] = lastindex;
1674                                                                         fGeneratorNames[number] = GeneratorName;
1675                                                                         number++;
1676                                                                         continue;
1677                                                                 }
1678                                                         }
1679                                                 }
1680                                                 if ( fMCStackAOD){
1681                                                         AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindex));
1682                                                         if (  aodMCParticle->GetPdgCode() == fAddedSignalPDGCode ){
1683                                                                 if (periodName.CompareTo("LHC14a1b")==0 || periodName.CompareTo("LHC14a1c")==0 ){
1684                                                                         if (gh->NProduced() > 10) {
1685                                                                                 AliAODMCParticle *aodMCParticle2 = static_cast<AliAODMCParticle*>(fMCStackAOD->At(firstindex+10));
1686                                                                                 if ( aodMCParticle2->GetPdgCode() == fAddedSignalPDGCode ){
1687                                                                                         fNotRejectedEnd[number] = lastindex;
1688                                                                                         fNotRejectedStart[number] = firstindex;
1689                                                                                         fGeneratorNames[number] = GeneratorName;
1690                                                                                         number++;
1691                                                                         }
1692                                                                                         continue;
1693                                                                                 }
1694                                                                 } else {
1695                                                                         fNotRejectedStart[number] = firstindex;
1696                                                                         fNotRejectedEnd[number] = lastindex;
1697                                                                         fGeneratorNames[number] = GeneratorName;
1698                                                                         number++;
1699                                                                         continue;
1700                                                                 }
1701                                                         }
1702                                                 }
1703                                                 continue;
1704                                         } else {
1705                                                 fNotRejectedStart[number] = firstindex;
1706                                                 fNotRejectedEnd[number] = lastindex;
1707                                                 fGeneratorNames[number] = GeneratorName;
1708 //                                              cout << "Number of particles produced for: " << i << "\t" << GeneratorName.Data() << "\t" << lastindex-firstindex+1 << endl;
1709                                                 number++;
1710                                                 continue;
1711                                         }
1712
1713                                 }
1714                         }
1715                         firstindex = firstindex + gh->NProduced();
1716                 }
1717 //              for (Int_t i = 0; i < number; i++){
1718 //                      cout << i << "\t" <<fGeneratorNames[i] << "\t" << fNotRejectedStart[i] << "\t" <<fNotRejectedEnd[i] << endl;
1719 //              }
1720
1721         } else { // No Cocktail Header Found
1722                 fNotRejectedStart = new Int_t[1];
1723                 fNotRejectedEnd = new Int_t[1];
1724
1725                 fnHeaders = 1;
1726                 fNotRejectedStart[0] = 0;
1727                 fNotRejectedEnd[0] = static_cast<AliMCEvent*>(MCEvent)->Stack()->GetNprimary()-1;
1728                 fGeneratorNames = new TString[1];
1729                 fGeneratorNames[0] = "NoCocktailGeneratorFound";
1730
1731                 AliGenPythiaEventHeader *mcHeaderPythia = dynamic_cast<AliGenPythiaEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
1732                 if (mcHeaderPythia) fGeneratorNames[0] = "NoCocktailGeneratorFound_Pythia";
1733                 AliGenDPMjetEventHeader *mcHeaderPhojet = dynamic_cast<AliGenDPMjetEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
1734                 if (mcHeaderPhojet) fGeneratorNames[0] = "NoCocktailGeneratorFound_Phojet";
1735                 AliGenHijingEventHeader *mcHeaderHijing = dynamic_cast<AliGenHijingEventHeader*>(static_cast<AliMCEvent*>(MCEvent)->GenEventHeader());
1736                 if (mcHeaderHijing) fGeneratorNames[0] = "NoCocktailGeneratorFound_Hijing";
1737
1738                 SetRejectExtraSignalsCut(0);
1739         }
1740
1741 }
1742
1743 //_________________________________________________________________________
1744 Int_t AliConvEventCuts::IsParticleFromBGEvent(Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
1745
1746         // Not Accepted == kFALSE == 0
1747         //     Accepted ==  kTRUE == 1
1748         //  FirstHeader ==  kTRUE == 3
1749         if(index < 0) return 0; // No Particle
1750
1751 //      if (index == 100){
1752 //              cout << "possible headers" << endl;
1753 //              for(Int_t i = 0;i<fnHeaders;i++){
1754 //                      cout << i << "\t" <<fGeneratorNames[i] << "\t" << fNotRejectedStart[i] << "\t" <<fNotRejectedEnd[i] << endl;
1755 //              }
1756 //      }
1757         Int_t accepted = 0;
1758         if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
1759                 if( index >= MCStack->GetNprimary()){ // Secondary Particle
1760                         if( ((TParticle*)MCStack->Particle(index))->GetMother(0) < 0) return 1; // Secondary Particle without Mother??
1761                         return IsParticleFromBGEvent(((TParticle*)MCStack->Particle(index))->GetMother(0),MCStack,InputEvent);
1762                 }
1763                 for(Int_t i = 0;i<fnHeaders;i++){
1764                         if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
1765                                 accepted = 1;
1766                                 if(i == 0) accepted = 2; // MB Header
1767                         }
1768                 }
1769         }
1770         else if(InputEvent->IsA()==AliAODEvent::Class()){
1771                 TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1772                 if (AODMCTrackArray){
1773                         AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
1774                         if(!aodMCParticle) return 1; // Photon Without a Mother ? --> Accepted
1775                         if(!aodMCParticle->IsPrimary()){
1776                                 if( aodMCParticle->GetMother() < 0) return 1;// Secondary Particle without Mother??
1777                                 return IsParticleFromBGEvent(aodMCParticle->GetMother(),MCStack,InputEvent);
1778                         }
1779                         index = abs(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index))->GetLabel());
1780                         for(Int_t i = 0;i<fnHeaders;i++){
1781                                 if(index >= fNotRejectedStart[i] && index <= fNotRejectedEnd[i]){
1782                                         accepted = 1;
1783                                         if(i == 0) accepted = 2; // MB Header
1784                                 }
1785                         }
1786                 }
1787         }
1788
1789         return accepted;
1790 }
1791
1792 //_________________________________________________________________________
1793 Int_t AliConvEventCuts::IsEventAcceptedByCut(AliConvEventCuts *ReaderCuts, AliVEvent *InputEvent, AliMCEvent *MCEvent, Int_t isHeavyIon){
1794
1795         if ( !IsTriggerSelected(InputEvent) )
1796                 return 3;
1797
1798         if(isHeavyIon != 0 && !(IsCentralitySelected(InputEvent,MCEvent)))
1799                 return 1; // Check Centrality --> Not Accepted => eventQuality = 1
1800
1801         if(isHeavyIon == 0 && GetIsFromPileup()){
1802                 if(InputEvent->IsPileupFromSPD(3,0.8,3.,2.,5.)){
1803                         return 6; // Check Pileup --> Not Accepted => eventQuality = 6
1804                 }
1805         }
1806
1807         Bool_t hasV0And = ReaderCuts->HasV0AND();
1808         Bool_t isSDDFired = ReaderCuts->IsSDDFired();
1809         if( (IsSpecialTrigger() == 2 || IsSpecialTrigger() == 3) && !isSDDFired && !MCEvent)
1810                 return 7; // With SDD requested but no fired
1811
1812         if( (IsSpecialTrigger() == 1 || IsSpecialTrigger() == 3) && !hasV0And)
1813                 return 8; // V0AND requested but no fired
1814
1815         if(hCentrality)hCentrality->Fill(GetCentrality(InputEvent));
1816         if(hCentralityVsNumberOfPrimaryTracks)
1817                 hCentralityVsNumberOfPrimaryTracks->Fill(GetCentrality(InputEvent),
1818                                                                                                 ((AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()
1819                                                                                                 ->GetTask(fV0ReaderName.Data()))->GetNumberOfPrimaryTracks());
1820
1821         return 0;
1822 }
1823
1824 //_________________________________________________________________________
1825 Float_t AliConvEventCuts::GetWeightForMeson(TString period, Int_t index, AliStack *MCStack, AliVEvent *InputEvent){
1826         if (!(period.CompareTo("LHC12f1a") == 0 || period.CompareTo("LHC12f1b") == 0  || period.CompareTo("LHC12i3") == 0 || period.CompareTo("LHC11a10a") == 0 || period.CompareTo("LHC11a10b") == 0 || period.CompareTo("LHC11a10b_bis") == 0 || period.CompareTo("LHC11a10a_bis") == 0 || period.CompareTo("LHC11a10b_plus") == 0 || period.Contains("LHC13d2")|| period.Contains("LHC14a1") ||
1827         period.CompareTo("LHC13e7") == 0 || period.Contains("LHC13b2_efix") || period.CompareTo("LHC14b2") == 0 )) return 1.;
1828
1829         Int_t kCaseGen = 0;
1830         for (Int_t i = 0; i < fnHeaders; i++){
1831                 if (index >= fNotRejectedStart[i] && index < fNotRejectedEnd[i]+1){
1832                         if (fGeneratorNames[i].CompareTo("Pythia") == 0){
1833                                 kCaseGen = 1;
1834                         } else if (fGeneratorNames[i].CompareTo("DPMJET") == 0){
1835                                 kCaseGen = 2;
1836                         } else if (fGeneratorNames[i].CompareTo("HIJING") == 0 ||
1837                                                 fGeneratorNames[i].CompareTo("Hijing") == 0 ||
1838                                                 fGeneratorNames[i].Contains("hijing")){
1839                                 kCaseGen = 3;
1840                         } else if (fGeneratorNames[i].CompareTo("BOX") == 0){
1841                                 kCaseGen = 4;
1842                         } else if (fGeneratorNames[i].CompareTo("PARAM") == 0){
1843                                 kCaseGen = 5;
1844                         } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound") == 0){
1845                                 kCaseGen = 6;
1846                         } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Pythia") == 0){
1847                                 kCaseGen = 1;
1848                         } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Phojet") == 0){
1849                                 kCaseGen = 2;
1850                         } else if (fGeneratorNames[i].CompareTo("NoCocktailGeneratorFound_Hijing") == 0){
1851                                 kCaseGen = 3;
1852                         }
1853                         if (period.Contains("LHC13d2") || period.CompareTo("LHC13e7") == 0 || period.Contains("LHC13b2_efix")  || period.Contains("LHC14a1") || period.CompareTo("LHC14b2") == 0 ){
1854                                 kCaseGen = 3;
1855                         }
1856                 }
1857         }
1858         if (kCaseGen == 0) return 1;
1859
1860
1861         Double_t mesonPt = 0;
1862         Double_t mesonMass = 0;
1863         Int_t PDGCode = 0;
1864         if(!InputEvent || InputEvent->IsA()==AliESDEvent::Class()){
1865                 mesonPt = ((TParticle*)MCStack->Particle(index))->Pt();
1866                 mesonMass = ((TParticle*)MCStack->Particle(index))->GetCalcMass();
1867                 PDGCode = ((TParticle*)MCStack->Particle(index))->GetPdgCode();
1868         }
1869         else if(InputEvent->IsA()==AliAODEvent::Class()){
1870                 TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(InputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1871                 if (AODMCTrackArray){
1872                         AliAODMCParticle *aodMCParticle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(index));
1873                         mesonPt = aodMCParticle->Pt();
1874                         mesonMass = aodMCParticle->GetCalcMass();
1875                         PDGCode = aodMCParticle->GetPdgCode();
1876                 } else {
1877                         return 1;
1878                 }
1879         }
1880
1881         Float_t functionResultMC = 1.;
1882         if (kCaseGen == 1){ // Pythia 6
1883                 Float_t dNdyMC = 2.1462;
1884                 Float_t nMC = 7.06055;
1885                 Float_t tMC = 0.12533;
1886                 if ( PDGCode ==  111){
1887                         dNdyMC = 2.1462;
1888                         nMC = 7.06055;
1889                         tMC = 0.12533;
1890                 } else if ( PDGCode ==  221){
1891                         dNdyMC = 0.2357;
1892                         nMC = 5.9105;
1893                         tMC = 0.1525;
1894                 }
1895                 functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
1896         } else if (kCaseGen == 2){ // Phojet
1897                 Float_t dNdyMC = 2.35978;
1898                 Float_t nMC = 6.81795;
1899                 Float_t tMC = 0.11492;
1900                 if ( PDGCode ==  111){
1901                         dNdyMC = 2.35978;
1902                         nMC = 6.81795;
1903                         tMC = 0.11492;
1904                 } else if ( PDGCode ==  221){
1905                         dNdyMC = 0.3690;
1906                         nMC = 5.55809;
1907                         tMC = 0.13387;
1908                 }
1909                 functionResultMC = dNdyMC / ( 2 * TMath::Pi())*(nMC-1.)*(nMC-2.) / (nMC*tMC*(nMC*tMC+mesonMass*(nMC-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nMC*tMC), -nMC);
1910         } else if (kCaseGen == 4){ // BOX generators pp
1911         //       functionResultMC = 1./sqrt(1.-mesonMass*mesonMass/((mesonMass*mesonMass+mesonPt*mesonPt)*cosh(mesonY)*cosh(mesonY)));
1912                 Float_t a = 0.23437;
1913                 Float_t b = 5.6661;
1914                 Float_t c = -1430.5863;
1915                 Float_t d = -0.6966624;
1916                 Float_t e = 252.3742;
1917                 if ( PDGCode ==  111){
1918                         a = 0.23437;
1919                         b = 5.6661;
1920                         c = -1430.5863;
1921                         d = -0.6966624;
1922                         e = 252.3742;
1923                 } else if ( PDGCode ==  221){
1924                         a = 0.10399;
1925                         b = 4.35311;
1926                         c = -12.17723;
1927                         d = -0.01172;
1928                         e =1.85140;
1929                 }
1930                 functionResultMC = a*TMath::Power(mesonPt,-1.*(b+c/(TMath::Power(mesonPt,d)+e)))*1./mesonPt *1./1.6 *1./(2.* TMath::Pi());
1931         //       cout << functionResultMC << endl;
1932         } else if (kCaseGen == 3 ){ // HIJING
1933                 if ( PDGCode ==  111 && fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0){
1934                         functionResultMC = hReweightMCHistPi0->Interpolate(mesonPt);
1935                 }
1936                 if ( PDGCode ==  221 && fDoReweightHistoMCEta && hReweightMCHistEta!= 0x0){
1937                         functionResultMC = hReweightMCHistEta->Interpolate(mesonPt);
1938                 }
1939                 if ( PDGCode ==  310 && fDoReweightHistoMCK0s && hReweightMCHistK0s!= 0x0){
1940                         functionResultMC = hReweightMCHistK0s->Interpolate(mesonPt);
1941                 }
1942         }
1943
1944         Float_t functionResultData = 1;
1945         if (kCaseGen == 1 || kCaseGen == 2 || kCaseGen == 4 ){
1946                 Float_t dNdyData = 2.2328;
1947                 Float_t nData = 7.1473;
1948                 Float_t tData = 0.1346;
1949                 if ( PDGCode ==  111){
1950                         dNdyData = 2.2328;
1951                         nData = 7.1473;
1952                         tData = 0.1346;
1953                 } else if ( PDGCode ==  221){
1954                         dNdyData = 0.38992; //be careful this fit is not optimal, eta in data still has problems
1955                         nData = 5.72778;
1956                         tData = 0.13835;
1957                 }
1958                 functionResultData = dNdyData / ( 2 * TMath::Pi())*(nData-1.)*(nData-2.) / (nData*tData*(nData*tData+mesonMass*(nData-2.)))  * TMath::Power(1.+(TMath::Sqrt(mesonPt*mesonPt+mesonMass*mesonMass)-mesonMass)/(nData*tData), -nData);
1959         //       cout << functionResultData << endl;
1960         } else {
1961                 if ( PDGCode ==  111 && fDoReweightHistoMCPi0 && fFitDataPi0!= 0x0){
1962                         functionResultData = fFitDataPi0->Eval(mesonPt);
1963                 }
1964                 if ( PDGCode ==  221 && fDoReweightHistoMCEta && fFitDataEta!= 0x0){
1965                         functionResultData = fFitDataEta->Eval(mesonPt);
1966                 }
1967                 if ( PDGCode ==  310 && fDoReweightHistoMCK0s && fFitDataK0s!= 0x0){
1968                         functionResultData = fFitDataK0s->Eval(mesonPt);
1969                 }
1970
1971         }
1972
1973         Double_t weight = 1;
1974         if (PDGCode ==  111 || PDGCode ==  221){
1975                 if (functionResultData != 0. && functionResultMC != 0. && isfinite(functionResultData) && isfinite(functionResultMC)){
1976                         weight = functionResultData/functionResultMC;
1977                         if ( kCaseGen == 3){
1978                                 if (PDGCode ==  111){
1979                                 if (!(fDoReweightHistoMCPi0 && hReweightMCHistPi0!= 0x0 && PDGCode ==  111)){
1980                                         weight = 1.;
1981                                 }
1982                                 }
1983                                 if (PDGCode ==  221){
1984                                 if (!(fDoReweightHistoMCEta && hReweightMCHistEta!= 0x0 && PDGCode ==  221)){
1985                                         weight = 1.;
1986                                 }
1987                                 }
1988                         }
1989                         if (!isfinite(functionResultData)) weight = 1.;
1990                         if (!isfinite(weight)) weight = 1.;
1991                 }
1992         } else if (PDGCode ==  310 && functionResultMC != 0 && isfinite(functionResultMC)){
1993                         weight = functionResultMC;
1994         }
1995
1996         //    if (fModCentralityClass == 0 && fCentralityMin == 4 && fCentralityMax == 6 && PDGCode ==  111){
1997         //        cout << period.Data() << "\t" << kCaseGen << "\t" <<fModCentralityClass<< "\t" <<fCentralityMin<< "\t" <<fCentralityMax << "\t" << mesonPt << "\t" <<mesonMass<< "\t"<<functionResultData << "\t"<< functionResultMC << "\t" << weight <<endl;
1998         //    }
1999         return weight;
2000 }
2001
2002
2003 ///________________________________________________________________________
2004 void AliConvEventCuts::GetCorrectEtaShiftFromPeriod(TString periodName){
2005
2006    if(periodName.CompareTo("LHC12g") == 0 || //pilot run 2012
2007       periodName.CompareTo("LHC13b") == 0 || //mainly minimum bias
2008       periodName.CompareTo("LHC13c") == 0 || //mainly minimum bias
2009       periodName.CompareTo("LHC13d") == 0 || //mainly triggered
2010       periodName.CompareTo("LHC13e") == 0 || //mainly triggered
2011       periodName.CompareTo("LHC13c3") == 0 || //MC Starlight, anchor LHC13d+e
2012       periodName.CompareTo("LHC13c2") == 0 || //MC Starlight, coherent J/Psi, UPC muon anchor LHC13d+e
2013       periodName.CompareTo("LHC13b4") == 0 || //MC Pythia 6 (Jet-Jet), anchor LHC13b
2014       periodName.CompareTo("LHC13b2_fix_1") == 0 || //MC DPMJET, anchr LHC13b+c
2015       periodName.CompareTo("LHC13b2_efix_p1") == 0 || //MC DPMJET, anchr LHC13b+c
2016       periodName.CompareTo("LHC13b2_efix_p2") == 0 || //MC DPMJET, anchr LHC13b+c
2017       periodName.CompareTo("LHC13b2_efix_p3") == 0 || //MC DPMJET, anchr LHC13b+c
2018       periodName.CompareTo("LHC13b2_efix_p4") == 0 || //MC DPMJET, anchr LHC13b+c
2019       periodName.CompareTo("LHC13e7") == 0 || //MC DPMJET, anchr LHC13b+c
2020       periodName.CompareTo("LHC13b3") == 0 || //MC HIJING, weighted to number of events per run, anchor LHC13b
2021       periodName.CompareTo("LHC13b2") == 0 ||  // MC DPMJET, wrong energy, anchor LHC13b
2022       periodName.CompareTo("LHC13b2_plus") == 0 || // MC DPMJET, weighted to number event per run, anchor LHC13b
2023       periodName.CompareTo("LHC13c1_bis") == 0 || // MC AMPT fast generation, pT hardbin, anchor ?
2024       periodName.CompareTo("LHC13c1") == 0 || // MC AMPT fast generation, anchor ?
2025       periodName.CompareTo("LHC13b1") == 0 || // MC DPMJET, fragments, with fixed label 0, anchor LHC12g
2026       periodName.CompareTo("LHC12g4b_fix") == 0 || // MC DPMJET, with fixed label 0, anchor LHC12g
2027       periodName.CompareTo("LHC12g1_fix") == 0 || // MC ?, with fixed label 0, anchor LHC12g
2028       periodName.CompareTo("LHC12g4c") == 0 || // MC DPMJET, shifted vertex runs, anchor LHC12g
2029       periodName.CompareTo("LHC12h6") == 0 || // MC muon cocktail, anchor LHC12g
2030       periodName.CompareTo("LHC12g4b") == 0 || // MC DPMJET 3rd iteration, anchor LHC12g
2031       periodName.CompareTo("LHC12g4a") == 0 || // MC DPMJET improved, anchor LHC12g
2032       periodName.CompareTo("LHC12g4") == 0 || // MC DPMJET, anchor LHC12g
2033       periodName.CompareTo("LHC12g5") == 0 || // MC PHOJET, anchor LHC12g
2034       periodName.CompareTo("LHC12g2") == 0 || // MC Starlight background, anchor LHC12g
2035       periodName.CompareTo("LHC12g1") == 0 ) // MC ?, anchor LHC12g
2036       {
2037          printf(" Gamma Conversion Cuts %s :: pPb Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),-0.465);
2038          SetEtaShift(-0.465);
2039       }
2040    else if(periodName.CompareTo("LHC13f") == 0 ||
2041            periodName.CompareTo("LHC13c6b") == 0 ||// MC Jpsi -> mumu, anchor LHC13f
2042            periodName.CompareTo("LHC13c5") == 0 || //MC Starlight, gamma gamma UPC muon, anchor LHC13f
2043            periodName.CompareTo("LHC13c4") == 0 )//MC Starlight, coherent JPsi, UPC muon, anchor LHC13f
2044       {
2045          printf(" Gamma Conversion Cuts %s :: Pbp Run doing Eta Shift of %f \n\n",(GetCutNumber()).Data(),0.465);
2046          SetEtaShift(+0.465);
2047       }
2048    else printf(" Gamma Conversion Cuts %s :: Automatic Eta Shift requested but Period is not known -> No Shift \n\n",(GetCutNumber()).Data());
2049 }
2050