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