PWGGA/GammaConv/AliConversionMesonCuts.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 "AliConversionMesonCuts.h"
  24
  25 #include "AliKFVertex.h"
  26 #include "AliAODTrack.h"
  27 #include "AliESDtrack.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliInputEventHandler.h"
  30 #include "AliMCEventHandler.h"
  31 #include "AliAODHandler.h"
  32 #include "AliPIDResponse.h"
  33 #include "TH1.h"
  34 #include "TH2.h"
  35 #include "AliStack.h"
  36 #include "AliAODConversionMother.h"
  37 #include "TObjString.h"
  38 #include "AliAODEvent.h"
  39 #include "AliESDEvent.h"
  40 #include "AliCentrality.h"
  41 #include "TList.h"
  42 #include "TPDGCode.h"
  43 #include "TDatabasePDG.h"
  44 #include "AliAODMCParticle.h"
  45
  46 class iostream;
  47
  48 using namespace std;
  49
  50 ClassImp(AliConversionMesonCuts)
  51
  52
  53 const char* AliConversionMesonCuts::fgkCutNames[AliConversionMesonCuts::kNCuts] = {
  54    "MesonKind", //0
  55    "BackgroundScheme", //1
  56    "NumberOfBGEvents", //2
  57    "DegreesForRotationMethod", //3
  58    "RapidityMesonCut", //4
  59    "RCut", //5
  60    "AlphaMesonCut", //6
  61    "SelectionWindow", //7
  62    "SharedElectronCuts", //8
  63    "RejectToCloseV0s", //9
  64    "UseMCPSmearing", //10
  65    "DcaGammaGamma", //11
  66    "DcaRPrimVtx", //12
  67    "DcaZPrimVtx" //13
  68 };
  69
  70
  71 //________________________________________________________________________
  72 AliConversionMesonCuts::AliConversionMesonCuts(const char *name,const char *title) :
  73    AliAnalysisCuts(name,title),
  74    fHistograms(NULL),
  75    fMesonKind(0),
  76    fMaxR(200),
  77    fSelectionLow(0.08),
  78    fSelectionHigh(0.145),
  79    fAlphaMinCutMeson(0),
  80    fAlphaCutMeson(1),
  81    fRapidityCutMeson(1),
  82    fUseRotationMethodInBG(kFALSE),
  83    fDoBG(kTRUE),
  84    fdoBGProbability(kFALSE),
  85    fUseTrackMultiplicityForBG(kFALSE),
  86    fnDegreeRotationPMForBG(0),
  87    fNumberOfBGEvents(0),
  88    fOpeningAngle(0.005),
  89    fDoToCloseV0sCut(kFALSE),
  90    fminV0Dist(200.),
  91    fDoSharedElecCut(kFALSE),
  92    fUseMCPSmearing(kFALSE),
  93    fPBremSmearing(0),
  94    fPSigSmearing(0),
  95    fPSigSmearingCte(0),
  96    fBrem(NULL),
  97    fRandom(0),
  98    fFAlphaCut(0),
  99    fAlphaPtDepCut(kFALSE),
 100    fElectronLabelArraySize(500),
 101    fElectronLabelArray(NULL),
 102    fDCAGammaGammaCut(1000),
 103    fDCAZMesonPrimVtxCut(1000),
 104    fDCARMesonPrimVtxCut(1000),
 105    fBackgroundHandler(0),
 106    fCutString(NULL),
 107    hMesonCuts(NULL),
 108    hMesonBGCuts(NULL),
 109    hDCAGammaGammaMesonBefore(NULL),
 110    hDCAZMesonPrimVtxBefore(NULL),
 111    hDCARMesonPrimVtxBefore(NULL),
 112    hDCAGammaGammaMesonAfter(NULL),
 113    hDCAZMesonPrimVtxAfter(NULL),
 114    hDCARMesonPrimVtxAfter(NULL)
 115
 116 {
 117    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=0;}
 118    fCutString=new TObjString((GetCutNumber()).Data());
 119    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
 120    if (fBrem == NULL){
 121       fBrem = new TF1("fBrem","pow(-log(x),[0]/log(2.0)-1.0)/TMath::Gamma([0]/log(2.0))",0.00001,0.999999999);
 122       // tests done with 1.0e-14
 123       fBrem->SetParameter(0,fPBremSmearing);
 124       fBrem->SetNpx(100000);
 125    }
 126
 127 }
 128
 129 //________________________________________________________________________
 130 AliConversionMesonCuts::AliConversionMesonCuts(const AliConversionMesonCuts &ref) :
 131    AliAnalysisCuts(ref),
 132    fHistograms(NULL),
 133    fMesonKind(ref.fMesonKind),
 134    fMaxR(ref.fMaxR),
 135    fSelectionLow(ref.fSelectionLow),
 136    fSelectionHigh(ref.fSelectionHigh),
 137    fAlphaMinCutMeson(ref.fAlphaMinCutMeson),
 138    fAlphaCutMeson(ref.fAlphaCutMeson),
 139    fRapidityCutMeson(ref.fRapidityCutMeson),
 140    fUseRotationMethodInBG(ref.fUseRotationMethodInBG),
 141    fDoBG(ref.fDoBG),
 142    fdoBGProbability(ref.fdoBGProbability),
 143    fUseTrackMultiplicityForBG(ref.fUseTrackMultiplicityForBG),
 144    fnDegreeRotationPMForBG(ref.fnDegreeRotationPMForBG),
 145    fNumberOfBGEvents(ref. fNumberOfBGEvents),
 146    fOpeningAngle(ref.fOpeningAngle),
 147    fDoToCloseV0sCut(ref.fDoToCloseV0sCut),
 148    fminV0Dist(ref.fminV0Dist),
 149    fDoSharedElecCut(ref.fDoSharedElecCut),
 150    fUseMCPSmearing(ref.fUseMCPSmearing),
 151    fPBremSmearing(ref.fPBremSmearing),
 152    fPSigSmearing(ref.fPSigSmearing),
 153    fPSigSmearingCte(ref.fPSigSmearingCte),
 154    fBrem(NULL),
 155    fRandom(ref.fRandom),
 156    fFAlphaCut(NULL),
 157    fAlphaPtDepCut(ref.fAlphaPtDepCut),
 158    fElectronLabelArraySize(ref.fElectronLabelArraySize),
 159    fElectronLabelArray(NULL),
 160    fDCAGammaGammaCut(ref.fDCAGammaGammaCut),
 161    fDCAZMesonPrimVtxCut(ref.fDCAZMesonPrimVtxCut),
 162    fDCARMesonPrimVtxCut(ref.fDCARMesonPrimVtxCut),
 163    fBackgroundHandler(ref.fBackgroundHandler),
 164    fCutString(NULL),
 165    hMesonCuts(NULL),
 166    hMesonBGCuts(NULL),
 167    hDCAGammaGammaMesonBefore(NULL),
 168    hDCAZMesonPrimVtxBefore(NULL),
 169    hDCARMesonPrimVtxBefore(NULL),
 170    hDCAGammaGammaMesonAfter(NULL),
 171    hDCAZMesonPrimVtxAfter(NULL),
 172    hDCARMesonPrimVtxAfter(NULL)
 173 {
 174    // Copy Constructor
 175    for(Int_t jj=0;jj<kNCuts;jj++){fCuts[jj]=ref.fCuts[jj];}
 176    fCutString=new TObjString((GetCutNumber()).Data());
 177    fElectronLabelArray = new Int_t[fElectronLabelArraySize];
 178    if (fBrem == NULL)fBrem = (TF1*)ref.fBrem->Clone("fBrem");
 179    // Histograms are not copied, if you need them, call InitCutHistograms
 180 }
 181
 182
 183 //________________________________________________________________________
 184 AliConversionMesonCuts::~AliConversionMesonCuts() {
 185    // Destructor
 186    //Deleting fHistograms leads to seg fault it it's added to output collection of a task
 187    // if(fHistograms)
 188    //   delete fHistograms;
 189    // fHistograms = NULL;
 190    if(fCutString != NULL){
 191       delete fCutString;
 192       fCutString = NULL;
 193    }
 194    if(fElectronLabelArray){
 195       delete fElectronLabelArray;
 196       fElectronLabelArray = NULL;
 197    }
 198
 199 }
 200
 201 //________________________________________________________________________
 202 void AliConversionMesonCuts::InitCutHistograms(TString name, Bool_t additionalHists){
 203
 204    // Initialize Cut Histograms for QA (only initialized and filled if function is called)
 205    TH1::AddDirectory(kFALSE);
 206
 207    if(fHistograms != NULL){
 208       delete fHistograms;
 209       fHistograms=NULL;
 210    }
 211    if(fHistograms==NULL){
 212       fHistograms=new TList();
 213       fHistograms->SetOwner(kTRUE);
 214       if(name=="")fHistograms->SetName(Form("ConvMesonCuts_%s",GetCutNumber().Data()));
 215       else fHistograms->SetName(Form("%s_%s",name.Data(),GetCutNumber().Data()));
 216    }
 217
 218    // Meson Cuts
 219    hMesonCuts=new TH1F(Form("MesonCuts %s",GetCutNumber().Data()),"MesonCuts",13,-0.5,12.5);
 220    hMesonCuts->GetXaxis()->SetBinLabel(1,"in");
 221    hMesonCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 222    hMesonCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 223    hMesonCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 224    hMesonCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 225    hMesonCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 226    hMesonCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
 227    hMesonCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
 228    hMesonCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
 229    hMesonCuts->GetXaxis()->SetBinLabel(10,"out");
 230    fHistograms->Add(hMesonCuts);
 231
 232    hMesonBGCuts=new TH1F(Form("MesonBGCuts %s",GetCutNumber().Data()),"MesonBGCuts",13,-0.5,12.5);
 233    hMesonBGCuts->GetXaxis()->SetBinLabel(1,"in");
 234    hMesonBGCuts->GetXaxis()->SetBinLabel(2,"undef rapidity");
 235    hMesonBGCuts->GetXaxis()->SetBinLabel(3,"rapidity cut");
 236    hMesonBGCuts->GetXaxis()->SetBinLabel(4,"opening angle");
 237    hMesonBGCuts->GetXaxis()->SetBinLabel(5,"alpha max");
 238    hMesonBGCuts->GetXaxis()->SetBinLabel(6,"alpha min");
 239    hMesonBGCuts->GetXaxis()->SetBinLabel(7,"dca gamma gamma");
 240    hMesonBGCuts->GetXaxis()->SetBinLabel(8,"dca R prim Vtx");
 241    hMesonBGCuts->GetXaxis()->SetBinLabel(9,"dca Z prim Vtx");
 242    hMesonBGCuts->GetXaxis()->SetBinLabel(10,"out");
 243
 244    fHistograms->Add(hMesonBGCuts);
 245
 246    if (additionalHists){
 247       hDCAGammaGammaMesonBefore=new TH1F(Form("DCAGammaGammaMeson Before %s",GetCutNumber().Data()),"DCAGammaGammaMeson Before",200,0,10);
 248       fHistograms->Add(hDCAGammaGammaMesonBefore);
 249
 250       hDCARMesonPrimVtxBefore=new TH1F(Form("DCARMesonPrimVtx Before %s",GetCutNumber().Data()),"DCARMesonPrimVtx Before",200,0,10);
 251       fHistograms->Add(hDCARMesonPrimVtxBefore);
 252
 253       hDCAZMesonPrimVtxBefore=new TH1F(Form("DCAZMesonPrimVtx Before %s",GetCutNumber().Data()),"DCAZMesonPrimVtx Before",401,-10,10);
 254       fHistograms->Add(hDCAZMesonPrimVtxBefore);
 255
 256    }
 257
 258    hDCAGammaGammaMesonAfter=new TH1F(Form("DCAGammaGammaMeson After %s",GetCutNumber().Data()),"DCAGammaGammaMeson After",200,0,10);
 259    fHistograms->Add(hDCAGammaGammaMesonAfter);
 260
 261    hDCAZMesonPrimVtxAfter=new TH2F(Form("InvMassDCAZMesonPrimVtx After %s",GetCutNumber().Data()),"InvMassDCAZMesonPrimVtx After",800,0,0.8,401,-10,10);
 262    fHistograms->Add(hDCAZMesonPrimVtxAfter);
 263
 264    hDCARMesonPrimVtxAfter=new TH1F(Form("DCARMesonPrimVtx After %s",GetCutNumber().Data()),"DCARMesonPrimVtx After",200,0,10);
 265    fHistograms->Add(hDCARMesonPrimVtxAfter);
 266
 267    TH1::AddDirectory(kTRUE);
 268 }
 269
 270 //________________________________________________________________________
 271 Bool_t AliConversionMesonCuts::MesonIsSelectedMC(TParticle *fMCMother,AliStack *fMCStack, Double_t fRapidityShift){
 272    // Returns true for all pions within acceptance cuts for decay into 2 photons
 273    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 274
 275    if(!fMCStack)return kFALSE;
 276
 277    if(fMCMother->GetPdgCode()==111 || fMCMother->GetPdgCode()==221){
 278       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 279
 280       Double_t rapidity = 10.;
 281       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 282          rapidity=8.-fRapidityShift;
 283       } else{
 284          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 285       }
 286
 287       // Rapidity Cut
 288       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 289
 290       // Select only -> 2y decay channel
 291       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 292
 293       for(Int_t i=0;i<2;i++){
 294          TParticle *MDaughter=fMCStack->Particle(fMCMother->GetDaughter(i));
 295          // Is Daughter a Photon?
 296          if(MDaughter->GetPdgCode()!=22)return kFALSE;
 297          // Is Photon in Acceptance?
 298          //   if(bMCDaughtersInAcceptance){
 299          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
 300          //   }
 301       }
 302       return kTRUE;
 303    }
 304    return kFALSE;
 305 }
 306 //________________________________________________________________________
 307 Bool_t AliConversionMesonCuts::MesonIsSelectedAODMC(AliAODMCParticle *MCMother,TClonesArray *AODMCArray, Double_t fRapidityShift){
 308    // Returns true for all pions within acceptance cuts for decay into 2 photons
 309    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 310
 311    if(!AODMCArray)return kFALSE;
 312
 313    if(MCMother->GetPdgCode()==111 || MCMother->GetPdgCode()==221){
 314       Double_t rMeson = sqrt( (MCMother->Xv()*MCMother->Xv()) + (MCMother->Yv()*MCMother->Yv()) ) ;
 315       if(rMeson>fMaxR)  return kFALSE; // cuts on distance from collision point
 316
 317       Double_t rapidity = 10.;
 318       if(MCMother->E() - MCMother->Pz() == 0 || MCMother->E() + MCMother->Pz() == 0){
 319          rapidity=8.-fRapidityShift;
 320       } else{
 321          rapidity = 0.5*(TMath::Log((MCMother->E()+MCMother->Pz()) / (MCMother->E()-MCMother->Pz())))-fRapidityShift;
 322       }
 323
 324       // Rapidity Cut
 325       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 326
 327       // Select only -> 2y decay channel
 328       if(MCMother->GetNDaughters()!=2)return kFALSE;
 329
 330       for(Int_t i=0;i<2;i++){
 331          AliAODMCParticle *MDaughter=static_cast<AliAODMCParticle*>(AODMCArray->At(MCMother->GetDaughter(i)));
 332          // Is Daughter a Photon?
 333          if(MDaughter->GetPdgCode()!=22)return kFALSE;
 334          // Is Photon in Acceptance?
 335          //   if(bMCDaughtersInAcceptance){
 336          //     if(!PhotonIsSelectedMC(MDaughter,fMCStack)){return kFALSE;}
 337          //   }
 338       }
 339       return kTRUE;
 340    }
 341    return kFALSE;
 342 }
 343
 344 //________________________________________________________________________
 345 Bool_t AliConversionMesonCuts::MesonIsSelectedMCDalitz(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelelectron, Int_t &labelpositron, Int_t &labelgamma, Double_t fRapidityShift){
 346
 347    // Returns true for all pions within acceptance cuts for decay into 2 photons
 348    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 349
 350    if( !fMCStack )return kFALSE;
 351
 352    if(  fMCMother->GetPdgCode() != 111 && fMCMother->GetPdgCode() != 221 ) return kFALSE;
 353
 354    if(  fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 355
 356    Double_t rapidity = 10.;
 357
 358    if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 359       rapidity=8.-fRapidityShift;
 360    }
 361    else{
 362       rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 363    }
 364
 365    // Rapidity Cut
 366    if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 367
 368    // Select only -> Dalitz decay channel
 369    if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 370
 371    TParticle *positron = 0x0;
 372    TParticle *electron = 0x0;
 373    TParticle    *gamma = 0x0;
 374
 375    for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 376
 377       TParticle* temp = (TParticle*)fMCStack->Particle( index );
 378
 379       switch( temp->GetPdgCode() ) {
 380       case ::kPositron:
 381          positron      =  temp;
 382          labelpositron = index;
 383          break;
 384       case ::kElectron:
 385          electron      =  temp;
 386          labelelectron = index;
 387          break;
 388       case ::kGamma:
 389          gamma         =  temp;
 390          labelgamma    = index;
 391          break;
 392       }
 393    }
 394
 395    if( positron && electron && gamma) return kTRUE;
 396    return kFALSE;
 397 }
 398
 399 //________________________________________________________________________
 400 Bool_t AliConversionMesonCuts::MesonIsSelectedMCEtaPiPlPiMiGamma(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelGamma, Double_t fRapidityShift){
 401
 402         // Returns true for all pions within acceptance cuts for decay into 2 photons
 403         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 404
 405         if( !fMCStack )return kFALSE;
 406
 407         if( fMCMother->GetPdgCode() != 221 ) return kFALSE;
 408
 409         if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 410
 411         Double_t rapidity = 10.;
 412
 413         if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 414                 rapidity=8.-fRapidityShift;
 415         }
 416         else{
 417                 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 418         }
 419
 420         // Rapidity Cut
 421         if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 422
 423         // Select only -> Dalitz decay channel
 424         if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 425
 426         TParticle *posPion = 0x0;
 427         TParticle *negPion = 0x0;
 428         TParticle    *gamma = 0x0;
 429
 430         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 431
 432                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
 433
 434                 switch( temp->GetPdgCode() ) {
 435                 case 211:
 436                         posPion      =  temp;
 437                         labelPosPion = index;
 438                         break;
 439                 case -211:
 440                         negPion      =  temp;
 441                         labelNegPion = index;
 442                         break;
 443                 case ::kGamma:
 444                         gamma         =  temp;
 445                         labelGamma    = index;
 446                         break;
 447                 }
 448         }
 449
 450         if( posPion && negPion && gamma) return kTRUE;
 451         return kFALSE;
 452 }
 453
 454 //________________________________________________________________________
 455 Bool_t AliConversionMesonCuts::MesonIsSelectedMCPiPlPiMiPiZero(TParticle *fMCMother,AliStack *fMCStack, Int_t &labelNegPion, Int_t &labelPosPion, Int_t &labelNeutPion, Double_t fRapidityShift){
 456
 457         // Returns true for all pions within acceptance cuts for decay into 2 photons
 458         // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 459
 460         if( !fMCStack )return kFALSE;
 461
 462         if( !(fMCMother->GetPdgCode() == 221 || fMCMother->GetPdgCode() == 223) ) return kFALSE;
 463
 464         if( fMCMother->R()>fMaxR ) return kFALSE; // cuts on distance from collision point
 465
 466         Double_t rapidity = 10.;
 467
 468         if( fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0 ){
 469                 rapidity=8.-fRapidityShift;
 470         }
 471         else{
 472                 rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 473         }
 474
 475         // Rapidity Cut
 476         if( abs(rapidity) > fRapidityCutMeson )return kFALSE;
 477
 478         // Select only -> pi+ pi- pi0
 479         if( fMCMother->GetNDaughters() != 3 )return kFALSE;
 480
 481         TParticle *posPion = 0x0;
 482         TParticle *negPion = 0x0;
 483         TParticle *neutPion = 0x0;
 484
 485         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 486
 487                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
 488                 switch( temp->GetPdgCode() ) {
 489                 case 211:
 490                         posPion      =  temp;
 491                         labelPosPion = index;
 492                         break;
 493                 case -211:
 494                         negPion      =  temp;
 495                         labelNegPion = index;
 496                         break;
 497                 case 111:
 498                         neutPion         =  temp;
 499                         labelNeutPion    = index;
 500                         break;
 501                 }
 502         }
 503
 504         if( posPion && negPion && neutPion ) return kTRUE;
 505         return kFALSE;
 506 }
 507
 508
 509 //________________________________________________________________________
 510 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
 511    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
 512    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 513
 514    if(!fMCStack)return kFALSE;
 515    // if(fMCMother->GetPdgCode()==20443 ){
 516    //    return kFALSE;
 517    // }
 518    if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
 519       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 520
 521       Double_t rapidity = 10.;
 522       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 523          rapidity=8.-fRapidityShift;
 524       }
 525       else{
 526          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 527       }
 528
 529       // Rapidity Cut
 530       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 531
 532       // Select only -> ChiC radiative (JPsi+gamma) decay channel
 533       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 534
 535       TParticle *jpsi   = 0x0;
 536       TParticle *gamma  = 0x0;
 537       TParticle *positron = 0x0;
 538       TParticle *electron = 0x0;
 539
 540       Int_t labeljpsiChiC = -1;
 541
 542       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 543
 544          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 545
 546          switch( temp->GetPdgCode() ) {
 547          case 443:
 548             jpsi =  temp;
 549             labeljpsiChiC = index;
 550             break;
 551          case 22:
 552             gamma    =  temp;
 553             labelgammaChiC = index;
 554             break;
 555          }
 556       }
 557
 558       if ( !jpsi || ! gamma) return kFALSE;
 559       if(jpsi->GetNDaughters()!=2)return kFALSE;
 560
 561
 562       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
 563          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 564          switch( temp->GetPdgCode() ) {
 565          case -11:
 566             electron =  temp;
 567             labelelectronChiC = index;
 568             break;
 569          case 11:
 570             positron =  temp;
 571             labelpositronChiC = index;
 572             break;
 573          }
 574       }
 575       if( !electron || !positron) return kFALSE;
 576       if( positron && electron && gamma) return kTRUE;
 577    }
 578    return kFALSE;
 579 }
 580
 581 ///________________________________________________________________________
 582 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
 583 {
 584
 585    // Selection of reconstructed Meson candidates
 586    // Use flag IsSignal in order to fill Fill different
 587    // histograms for Signal and Background
 588    TH1 *hist=0x0;
 589
 590    if(IsSignal){hist=hMesonCuts;}
 591    else{hist=hMesonBGCuts;}
 592
 593    Int_t cutIndex=0;
 594
 595    if(hist)hist->Fill(cutIndex);
 596    cutIndex++;
 597
 598    // Undefined Rapidity -> Floating Point exception
 599    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
 600       if(hist)hist->Fill(cutIndex);
 601       cutIndex++;
 602 //        cout << "undefined rapidity" << endl;
 603       return kFALSE;
 604    }
 605    else{
 606       // PseudoRapidity Cut --> But we cut on Rapidity !!!
 607       cutIndex++;
 608       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
 609          if(hist)hist->Fill(cutIndex);
 610 //               cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
 611          return kFALSE;
 612       }
 613    }
 614    cutIndex++;
 615
 616    // Opening Angle Cut
 617    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
 618    if( pi0->GetOpeningAngle() < fOpeningAngle){
 619 //        cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl;
 620       if(hist)hist->Fill(cutIndex);
 621       return kFALSE;
 622    }
 623    cutIndex++;
 624
 625    if ( fAlphaPtDepCut == kTRUE ) {
 626
 627         fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
 628    }
 629
 630
 631    // Alpha Max Cut
 632    if(pi0->GetAlpha()>fAlphaCutMeson){
 633 //         cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl;
 634       if(hist)hist->Fill(cutIndex);
 635       return kFALSE;
 636    }
 637    cutIndex++;
 638
 639    // Alpha Min Cut
 640    if(pi0->GetAlpha()<fAlphaMinCutMeson){
 641 //        cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl;
 642       if(hist)hist->Fill(cutIndex);
 643       return kFALSE;
 644    }
 645    cutIndex++;
 646
 647    if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
 648    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
 649
 650    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
 651 //        cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl;
 652       if(hist)hist->Fill(cutIndex);
 653       return kFALSE;
 654    }
 655    cutIndex++;
 656
 657    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
 658 //         cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl;
 659       if(hist)hist->Fill(cutIndex);
 660       return kFALSE;
 661    }
 662    cutIndex++;
 663
 664
 665    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
 666
 667    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
 668 //        cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl;
 669       if(hist)hist->Fill(cutIndex);
 670       return kFALSE;
 671    }
 672    cutIndex++;
 673
 674
 675    if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
 676    if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
 677    if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
 678
 679    if(hist)hist->Fill(cutIndex);
 680    return kTRUE;
 681 }
 682
 683
 684
 685 ///________________________________________________________________________
 686 ///________________________________________________________________________
 687 Bool_t AliConversionMesonCuts::UpdateCutString() {
 688    ///Update the cut string (if it has been created yet)
 689
 690    if(fCutString && fCutString->GetString().Length() == kNCuts) {
 691       fCutString->SetString(GetCutNumber());
 692    } else {
 693       return kFALSE;
 694    }
 695    return kTRUE;
 696 }
 697
 698 ///________________________________________________________________________
 699 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 700    // Initialize Cuts from a given Cut string
 701    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
 702    if(analysisCutSelection.Length()!=kNCuts) {
 703       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 704       return kFALSE;
 705    }
 706    if(!analysisCutSelection.IsDigit()){
 707       AliError("Cut selection contains characters");
 708       return kFALSE;
 709    }
 710
 711    const char *cutSelection = analysisCutSelection.Data();
 712 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 713    for(Int_t ii=0;ii<kNCuts;ii++){
 714       ASSIGNARRAY(ii);
 715    }
 716
 717    // Set Individual Cuts
 718    for(Int_t ii=0;ii<kNCuts;ii++){
 719       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 720    }
 721
 722    PrintCutsWithValues();
 723    return kTRUE;
 724 }
 725 ///________________________________________________________________________
 726 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
 727    ///Set individual cut ID
 728
 729    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 730    switch (cutID) {
 731    case kMesonKind:
 732       if( SetMesonKind(value)) {
 733          fCuts[kMesonKind] = value;
 734          UpdateCutString();
 735          return kTRUE;
 736       } else return kFALSE;
 737    case kSelectionCut:
 738       if( SetSelectionWindowCut(value)) {
 739          fCuts[kSelectionCut] = value;
 740          UpdateCutString();
 741          return kTRUE;
 742       } else return kFALSE;
 743    case kalphaMesonCut:
 744       if( SetAlphaMesonCut(value)) {
 745          fCuts[kalphaMesonCut] = value;
 746          UpdateCutString();
 747          return kTRUE;
 748       } else return kFALSE;
 749    case kRCut:
 750       if( SetRCut(value)) {
 751          fCuts[kRCut] = value;
 752          UpdateCutString();
 753          return kTRUE;
 754       } else return kFALSE;
 755
 756    case kRapidityMesonCut:
 757       if( SetRapidityMesonCut(value)) {
 758          fCuts[kRapidityMesonCut] = value;
 759          UpdateCutString();
 760          return kTRUE;
 761       } else return kFALSE;
 762
 763    case kBackgroundScheme:
 764       if( SetBackgroundScheme(value)) {
 765          fCuts[kBackgroundScheme] = value;
 766          UpdateCutString();
 767          return kTRUE;
 768       } else return kFALSE;
 769
 770    case kDegreesForRotationMethod:
 771       if( SetNDegreesForRotationMethod(value)) {
 772          fCuts[kDegreesForRotationMethod] = value;
 773          UpdateCutString();
 774          return kTRUE;
 775       } else return kFALSE;
 776
 777    case kNumberOfBGEvents:
 778       if( SetNumberOfBGEvents(value)) {
 779          fCuts[kNumberOfBGEvents] = value;
 780          UpdateCutString();
 781          return kTRUE;
 782       } else return kFALSE;
 783
 784    case kuseMCPSmearing:
 785       if( SetMCPSmearing(value)) {
 786          fCuts[kuseMCPSmearing] = value;
 787          UpdateCutString();
 788          return kTRUE;
 789       } else return kFALSE;
 790    case kElecShare:
 791       if( SetSharedElectronCut(value)) {
 792          fCuts[kElecShare] = value;
 793          UpdateCutString();
 794          return kTRUE;
 795       } else return kFALSE;
 796    case kToCloseV0s:
 797       if( SetToCloseV0sCut(value)) {
 798          fCuts[kToCloseV0s] = value;
 799          UpdateCutString();
 800          return kTRUE;
 801       } else return kFALSE;
 802    case kDcaGammaGamma:
 803       if( SetDCAGammaGammaCut(value)) {
 804          fCuts[kDcaGammaGamma] = value;
 805          UpdateCutString();
 806          return kTRUE;
 807       } else return kFALSE;
 808    case kDcaZPrimVtx:
 809       if( SetDCAZMesonPrimVtxCut(value)) {
 810          fCuts[kDcaZPrimVtx] = value;
 811          UpdateCutString();
 812          return kTRUE;
 813       } else return kFALSE;
 814    case kDcaRPrimVtx:
 815       if( SetDCARMesonPrimVtxCut(value)) {
 816          fCuts[kDcaRPrimVtx] = value;
 817          UpdateCutString();
 818          return kTRUE;
 819       } else return kFALSE;
 820
 821    case kNCuts:
 822       cout << "Error:: Cut id out of range"<< endl;
 823       return kFALSE;
 824    }
 825
 826    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 827    return kFALSE;
 828
 829 }
 830
 831
 832 ///________________________________________________________________________
 833 void AliConversionMesonCuts::PrintCuts() {
 834    // Print out current Cut Selection
 835    for(Int_t ic = 0; ic < kNCuts; ic++) {
 836       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 837    }
 838 }
 839
 840 ///________________________________________________________________________
 841 void AliConversionMesonCuts::PrintCutsWithValues() {
 842    // Print out current Cut Selection with values
 843         printf("\nMeson cutnumber \n");
 844         for(Int_t ic = 0; ic < kNCuts; ic++) {
 845                 printf("%d",fCuts[ic]);
 846         }
 847         printf("\n\n");
 848
 849         printf("Meson cuts \n");
 850         printf("\t |y| < %3.2f \n", fRapidityCutMeson);
 851         printf("\t theta_{open} < %3.2f\n", fOpeningAngle);
 852         if (!fAlphaPtDepCut) printf("\t %3.2f < alpha < %3.2f\n", fAlphaMinCutMeson, fAlphaCutMeson);
 853         printf("\t dca_{gamma,gamma} > %3.2f\n", fDCAGammaGammaCut);
 854         printf("\t dca_{R, prim Vtx} > %3.2f\n", fDCARMesonPrimVtxCut);
 855         printf("\t dca_{Z, prim Vtx} > %3.2f\n\n", fDCAZMesonPrimVtxCut);
 856         printf("\t Meson selection window for further analysis %3.3f > M_{gamma,gamma} > %3.3f\n\n", fSelectionLow, fSelectionHigh);
 857
 858          printf("Meson BG settings \n");
 859          if (!fDoBG){
 860                 if (!fUseRotationMethodInBG  & !fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing V0 mult \n");
 861                 if (!fUseRotationMethodInBG  & fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing track mult \n");
 862                 if (fUseRotationMethodInBG )printf("\t BG scheme: rotation \n");
 863                 if (fdoBGProbability) printf("\t -> use BG probability \n");
 864                 if (fBackgroundHandler) printf("\t -> use new BG handler \n");
 865                 printf("\t depth of pool: %d\n", fNumberOfBGEvents);
 866                 if (fUseRotationMethodInBG )printf("\t degree's for BG rotation: %d\n", fnDegreeRotationPMForBG);
 867          }
 868
 869 }
 870
 871
 872 ///________________________________________________________________________
 873 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
 874    // Set Cut
 875    switch(mesonKind){
 876    case 0:
 877       fMesonKind=0;
 878       break;
 879    case 1:
 880       fMesonKind=1;;
 881       break;
 882    default:
 883       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
 884       return kFALSE;
 885    }
 886    return kTRUE;
 887 }
 888
 889 ///________________________________________________________________________
 890 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
 891    // Set Cut
 892    switch(rCut){
 893    case 0:
 894       fMaxR = 180.;
 895       break;
 896    case 1:
 897       fMaxR = 180.;
 898       break;
 899    case 2:
 900       fMaxR = 180.;
 901       break;
 902    case 3:
 903       fMaxR = 70.;
 904       break;
 905    case 4:
 906       fMaxR = 70.;
 907       break;
 908    case 5:
 909       fMaxR = 180.;
 910       break;
 911       // High purity cuts for PbPb
 912    case 9:
 913       fMaxR = 180.;
 914       break;
 915    default:
 916       cout<<"Warning: rCut not defined"<<rCut<<endl;
 917       return kFALSE;
 918    }
 919    return kTRUE;
 920 }
 921
 922 ///________________________________________________________________________
 923 Bool_t AliConversionMesonCuts::SetSelectionWindowCut(Int_t selectionCut){
 924    // Set Cut
 925    switch(selectionCut){
 926    case 0:
 927       fSelectionLow = 0.08;
 928           fSelectionHigh = 0.145;
 929       break;
 930    case 1:
 931       fSelectionLow = 0.1;
 932           fSelectionHigh = 0.145;
 933       break;
 934    case 2:
 935           fSelectionLow = 0.11;
 936           fSelectionHigh = 0.145;
 937           break;
 938    case 3:
 939           fSelectionLow = 0.12;
 940           fSelectionHigh = 0.145;
 941       break;
 942    case 4:
 943           fSelectionLow = 0.1;
 944           fSelectionHigh = 0.15;
 945       break;
 946    case 5:
 947           fSelectionLow = 0.11;
 948           fSelectionHigh = 0.15;
 949       break;
 950    case 6:
 951           fSelectionLow = 0.12;
 952           fSelectionHigh = 0.15;
 953       break;
 954
 955    default:
 956       cout<<"Warning: SelectionCut not defined "<<selectionCut<<endl;
 957       return kFALSE;
 958    }
 959    return kTRUE;
 960 }
 961
 962
 963 ///________________________________________________________________________
 964 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 965 {   // Set Cut
 966    switch(alphaMesonCut){
 967    case 0:      // 0- 0.7
 968       fAlphaMinCutMeson  = 0.0;
 969       fAlphaCutMeson     = 0.7;
 970       fAlphaPtDepCut = kFALSE;
 971       break;
 972    case 1:      // Updated 31 October 2013 before 0.0 - 0.5
 973       if( fFAlphaCut ) delete fFAlphaCut;
 974       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 975       fFAlphaCut->SetParameter(0,0.7);
 976       fFAlphaCut->SetParameter(1,1.2);
 977       fAlphaMinCutMeson  =  0.0;
 978       fAlphaCutMeson     = -1.0;
 979       fAlphaPtDepCut = kTRUE;
 980       break;
 981    case 2:      // Updated 31 October 2013 before 0.5-1
 982       if( fFAlphaCut ) delete fFAlphaCut;
 983       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 984       fFAlphaCut->SetParameter(0,0.8);
 985       fFAlphaCut->SetParameter(1,1.2);
 986       fAlphaMinCutMeson  =  0.0;
 987       fAlphaCutMeson     = -1.0;
 988       fAlphaPtDepCut = kTRUE;
 989       break;
 990    case 3:      // 0.0-1
 991       fAlphaMinCutMeson  = 0.0;
 992       fAlphaCutMeson     = 1.;
 993       fAlphaPtDepCut = kFALSE;
 994       break;
 995    case 4:      // 0-0.65
 996       fAlphaMinCutMeson  = 0.0;
 997       fAlphaCutMeson     = 0.65;
 998       fAlphaPtDepCut = kFALSE;
 999       break;
1000    case 5:      // 0-0.75
1001       fAlphaMinCutMeson  = 0.0;
1002       fAlphaCutMeson     = 0.75;
1003       fAlphaPtDepCut = kFALSE;
1004       break;
1005    case 6:      // 0-0.8
1006       fAlphaMinCutMeson  = 0.0;
1007       fAlphaCutMeson     = 0.8;
1008       fAlphaPtDepCut = kFALSE;
1009       break;
1010    case 7:      // 0.0-0.85
1011       fAlphaMinCutMeson  = 0.0;
1012       fAlphaCutMeson     = 0.85;
1013       fAlphaPtDepCut = kFALSE;
1014       break;
1015    case 8:      // 0.0-0.6
1016       fAlphaMinCutMeson  = 0.0;
1017       fAlphaCutMeson     = 0.6;
1018       fAlphaPtDepCut = kFALSE;
1019       break;
1020    case 9: // Updated 11 November 2013 before 0.0 - 0.3
1021       if( fFAlphaCut ) delete fFAlphaCut;
1022       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
1023       fFAlphaCut->SetParameter(0,0.65);
1024       fFAlphaCut->SetParameter(1,1.2);
1025       fAlphaMinCutMeson  =  0.0;
1026       fAlphaCutMeson     = -1.0;
1027       fAlphaPtDepCut = kTRUE;
1028       break;
1029    default:
1030       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
1031       return kFALSE;
1032    }
1033    return kTRUE;
1034 }
1035
1036 ///________________________________________________________________________
1037 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
1038    // Set Cut
1039    switch(RapidityMesonCut){
1040    case 0:  // changed from 0.9 to 1.35
1041       fRapidityCutMeson   = 1.35;
1042       break;
1043    case 1:  //
1044       fRapidityCutMeson   = 0.8;
1045       break;
1046    case 2:  //
1047       fRapidityCutMeson   = 0.7;
1048       break;
1049    case 3:  //
1050       fRapidityCutMeson   = 0.6;
1051       break;
1052    case 4:  //
1053       fRapidityCutMeson   = 0.5;
1054       break;
1055    case 5:  //
1056       fRapidityCutMeson   = 0.85;
1057       break;
1058    case 6:  //
1059       fRapidityCutMeson   = 0.75;
1060       break;
1061    case 7:  //
1062       fRapidityCutMeson   = 0.3;
1063       break;
1064    case 8:  //changed, before 0.35
1065       fRapidityCutMeson   = 0.25;
1066       break;
1067    case 9:  //
1068       fRapidityCutMeson   = 0.4;
1069       break;
1070    default:
1071       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
1072       return kFALSE;
1073    }
1074    return kTRUE;
1075 }
1076
1077
1078 ///________________________________________________________________________
1079 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
1080    // Set Cut
1081    switch(BackgroundScheme){
1082    case 0: //Rotation
1083       fUseRotationMethodInBG=kTRUE;
1084       fdoBGProbability=kFALSE;
1085       break;
1086    case 1: // mixed event with V0 multiplicity
1087       fUseRotationMethodInBG=kFALSE;
1088       fUseTrackMultiplicityForBG=kFALSE;
1089       fdoBGProbability=kFALSE;
1090       break;
1091    case 2: // mixed event with track multiplicity
1092       fUseRotationMethodInBG=kFALSE;
1093       fUseTrackMultiplicityForBG=kTRUE;
1094       fdoBGProbability=kFALSE;
1095       break;
1096    case 3: //Rotation
1097       fUseRotationMethodInBG=kTRUE;
1098       fdoBGProbability=kTRUE;
1099       break;
1100    case 4: //No BG calculation
1101       cout << "no BG calculation should be done" << endl;
1102       fUseRotationMethodInBG=kFALSE;
1103       fdoBGProbability=kFALSE;
1104       fDoBG=kFALSE;
1105       break;
1106    case 5: //Rotation
1107       fUseRotationMethodInBG=kTRUE;
1108       fdoBGProbability=kFALSE;
1109       fBackgroundHandler = 1;
1110       break;
1111    case 6: // mixed event with V0 multiplicity
1112       fUseRotationMethodInBG=kFALSE;
1113       fUseTrackMultiplicityForBG=kFALSE;
1114       fdoBGProbability=kFALSE;
1115       fBackgroundHandler = 1;
1116       break;
1117    case 7: // mixed event with track multiplicity
1118       fUseRotationMethodInBG=kFALSE;
1119       fUseTrackMultiplicityForBG=kTRUE;
1120       fdoBGProbability=kFALSE;
1121       fBackgroundHandler = 1;
1122       break;
1123    case 8: //Rotation
1124       fUseRotationMethodInBG=kTRUE;
1125       fdoBGProbability=kTRUE;
1126       fBackgroundHandler = 1;
1127       break;
1128    default:
1129       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1130       return kFALSE;
1131    }
1132    return kTRUE;
1133 }
1134
1135
1136 ///________________________________________________________________________
1137 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
1138    // Set Cut
1139    switch(DegreesForRotationMethod){
1140    case 0:
1141       fnDegreeRotationPMForBG = 5;
1142       break;
1143    case 1:
1144       fnDegreeRotationPMForBG = 10;
1145       break;
1146    case 2:
1147       fnDegreeRotationPMForBG = 15;
1148       break;
1149    case 3:
1150       fnDegreeRotationPMForBG = 20;
1151       break;
1152    default:
1153       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1154       return kFALSE;
1155    }
1156    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1157    return kTRUE;
1158 }
1159
1160 ///________________________________________________________________________
1161 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
1162    // Set Cut
1163    switch(NumberOfBGEvents){
1164    case 0:
1165       fNumberOfBGEvents = 5;
1166       break;
1167    case 1:
1168       fNumberOfBGEvents = 10;
1169       break;
1170    case 2:
1171       fNumberOfBGEvents = 15;
1172       break;
1173    case 3:
1174       fNumberOfBGEvents = 20;
1175       break;
1176    case 4:
1177       fNumberOfBGEvents = 2;
1178       break;
1179    case 5:
1180       fNumberOfBGEvents = 50;
1181       break;
1182    case 6:
1183       fNumberOfBGEvents = 80;
1184       break;
1185    case 7:
1186       fNumberOfBGEvents = 100;
1187       break;
1188    default:
1189       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1190       return kFALSE;
1191    }
1192    return kTRUE;
1193 }
1194 ///________________________________________________________________________
1195 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1196
1197    switch(sharedElec){
1198    case 0:
1199       fDoSharedElecCut = kFALSE;
1200       break;
1201    case 1:
1202       fDoSharedElecCut = kTRUE;
1203       break;
1204    default:
1205       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1206       return kFALSE;
1207    }
1208
1209    return kTRUE;
1210 }
1211
1212 ///________________________________________________________________________
1213 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1214
1215    switch(toClose){
1216    case 0:
1217       fDoToCloseV0sCut = kFALSE;
1218       fminV0Dist = 250;
1219       break;
1220    case 1:
1221       fDoToCloseV0sCut = kTRUE;
1222       fminV0Dist = 1;
1223       break;
1224    case 2:
1225       fDoToCloseV0sCut = kTRUE;
1226       fminV0Dist = 2;
1227       break;
1228    case 3:
1229       fDoToCloseV0sCut = kTRUE;
1230       fminV0Dist = 3;
1231       break;
1232    default:
1233       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1234       return kFALSE;
1235    }
1236    return kTRUE;
1237 }
1238
1239 ///________________________________________________________________________
1240 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1241 {// Set Cut
1242    switch(useMCPSmearing){
1243    case 0:
1244       fUseMCPSmearing=0;
1245       fPBremSmearing=1.;
1246       fPSigSmearing=0.;
1247       fPSigSmearingCte=0.;
1248       break;
1249    case 1:
1250       fUseMCPSmearing=1;
1251       fPBremSmearing=1.0e-14;
1252       fPSigSmearing=0.;
1253       fPSigSmearingCte=0.;
1254       break;
1255    case 2:
1256       fUseMCPSmearing=1;
1257       fPBremSmearing=1.0e-15;
1258       fPSigSmearing=0.0;
1259       fPSigSmearingCte=0.;
1260       break;
1261    case 3:
1262       fUseMCPSmearing=1;
1263       fPBremSmearing=1.;
1264       fPSigSmearing=0.003;
1265       fPSigSmearingCte=0.002;
1266       break;
1267    case 4:
1268       fUseMCPSmearing=1;
1269       fPBremSmearing=1.;
1270       fPSigSmearing=0.003;
1271       fPSigSmearingCte=0.007;
1272       break;
1273    case 5:
1274       fUseMCPSmearing=1;
1275       fPBremSmearing=1.;
1276       fPSigSmearing=0.003;
1277       fPSigSmearingCte=0.016;
1278       break;
1279    case 6:
1280       fUseMCPSmearing=1;
1281       fPBremSmearing=1.;
1282       fPSigSmearing=0.007;
1283       fPSigSmearingCte=0.016;
1284       break;
1285    case 7:
1286       fUseMCPSmearing=1;
1287       fPBremSmearing=1.0e-16;
1288       fPSigSmearing=0.0;
1289       fPSigSmearingCte=0.;
1290       break;
1291    case 8:
1292       fUseMCPSmearing=1;
1293       fPBremSmearing=1.;
1294       fPSigSmearing=0.007;
1295       fPSigSmearingCte=0.014;
1296       break;
1297    case 9:
1298       fUseMCPSmearing=1;
1299       fPBremSmearing=1.;
1300       fPSigSmearing=0.007;
1301       fPSigSmearingCte=0.011;
1302       break;
1303
1304    default:
1305       AliError("Warning: UseMCPSmearing not defined");
1306       return kFALSE;
1307    }
1308    return kTRUE;
1309 }
1310
1311
1312 ///________________________________________________________________________
1313 Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1314    // Set Cut
1315    switch(DCAGammaGamma){
1316    case 0:  //
1317       fDCAGammaGammaCut   = 1000;
1318       break;
1319    case 1:  //
1320       fDCAGammaGammaCut   = 10;
1321       break;
1322    case 2:  //
1323       fDCAGammaGammaCut   = 5;
1324       break;
1325    case 3:  //
1326       fDCAGammaGammaCut   = 4;
1327       break;
1328    case 4:  //
1329       fDCAGammaGammaCut   = 3;
1330       break;
1331    case 5:  //
1332       fDCAGammaGammaCut   = 2.5;
1333       break;
1334    case 6:  //
1335       fDCAGammaGammaCut   = 2;
1336       break;
1337    case 7:  //
1338       fDCAGammaGammaCut   = 1.5;
1339       break;
1340    case 8:  //
1341       fDCAGammaGammaCut   = 1;
1342       break;
1343    case 9:  //
1344       fDCAGammaGammaCut   = 0.5;
1345       break;
1346    default:
1347       cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1348       return kFALSE;
1349    }
1350    return kTRUE;
1351 }
1352
1353 ///________________________________________________________________________
1354 Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1355    // Set Cut
1356    switch(DCAZMesonPrimVtx){
1357    case 0:  //
1358       fDCAZMesonPrimVtxCut   = 1000;
1359       break;
1360    case 1:  //
1361       fDCAZMesonPrimVtxCut   = 10;
1362       break;
1363    case 2:  //
1364       fDCAZMesonPrimVtxCut   = 5;
1365       break;
1366    case 3:  //
1367       fDCAZMesonPrimVtxCut   = 4;
1368       break;
1369    case 4:  //
1370       fDCAZMesonPrimVtxCut   = 3;
1371       break;
1372    case 5:  //
1373       fDCAZMesonPrimVtxCut   = 2.5;
1374       break;
1375    case 6:  //
1376       fDCAZMesonPrimVtxCut   = 2;
1377       break;
1378    case 7:  //
1379       fDCAZMesonPrimVtxCut   = 1.5;
1380       break;
1381    case 8:  //
1382       fDCAZMesonPrimVtxCut   = 1;
1383       break;
1384    case 9:  //
1385       fDCAZMesonPrimVtxCut   = 0.5;
1386       break;
1387    default:
1388       cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1389       return kFALSE;
1390    }
1391    return kTRUE;
1392 }
1393
1394 ///________________________________________________________________________
1395 Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1396    // Set Cut
1397    switch(DCARMesonPrimVtx){
1398    case 0:  //
1399       fDCARMesonPrimVtxCut   = 1000;
1400       break;
1401    case 1:  //
1402       fDCARMesonPrimVtxCut   = 10;
1403       break;
1404    case 2:  //
1405       fDCARMesonPrimVtxCut   = 5;
1406       break;
1407    case 3:  //
1408       fDCARMesonPrimVtxCut   = 4;
1409       break;
1410    case 4:  //
1411       fDCARMesonPrimVtxCut   = 3;
1412       break;
1413    case 5:  //
1414       fDCARMesonPrimVtxCut   = 2.5;
1415       break;
1416    case 6:  //
1417       fDCARMesonPrimVtxCut   = 2;
1418       break;
1419    case 7:  //
1420       fDCARMesonPrimVtxCut   = 1.5;
1421       break;
1422    case 8:  //
1423       fDCARMesonPrimVtxCut   = 1;
1424       break;
1425    case 9:  //
1426       fDCARMesonPrimVtxCut   = 0.5;
1427       break;
1428    default:
1429       cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1430       return kFALSE;
1431    }
1432    return kTRUE;
1433 }
1434
1435
1436 ///________________________________________________________________________
1437 TString AliConversionMesonCuts::GetCutNumber(){
1438    // returns TString with current cut number
1439    TString a(kNCuts);
1440    for(Int_t ii=0;ii<kNCuts;ii++){
1441       a.Append(Form("%d",fCuts[ii]));
1442    }
1443    return a;
1444 }
1445
1446 ///________________________________________________________________________
1447 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
1448
1449    Int_t posLabel = photon->GetTrackLabelPositive();
1450    Int_t negLabel = photon->GetTrackLabelNegative();
1451
1452    fElectronLabelArray[nV0*2] = posLabel;
1453    fElectronLabelArray[(nV0*2)+1] = negLabel;
1454 }
1455
1456 ///________________________________________________________________________
1457 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1458
1459    Int_t posLabel = photon->GetTrackLabelPositive();
1460    Int_t negLabel = photon->GetTrackLabelNegative();
1461
1462    for(Int_t i = 0; i<nV0s*2;i++){
1463       if(i==nV0*2)     continue;
1464       if(i==(nV0*2)+1) continue;
1465       if(fElectronLabelArray[i] == posLabel){
1466          return kFALSE;}
1467       if(fElectronLabelArray[i] == negLabel){
1468          return kFALSE;}
1469    }
1470
1471    return kTRUE;
1472 }
1473 ///________________________________________________________________________
1474 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
1475    Double_t posX = photon->GetConversionX();
1476    Double_t posY = photon->GetConversionY();
1477    Double_t posZ = photon->GetConversionZ();
1478
1479    for(Int_t i = 0;i<photons->GetEntries();i++){
1480       if(nV0 == i) continue;
1481       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1482       Double_t posCompX = photonComp->GetConversionX();
1483       Double_t posCompY = photonComp->GetConversionY();
1484       Double_t posCompZ = photonComp->GetConversionZ();
1485
1486       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
1487
1488       if(dist < fminV0Dist*fminV0Dist){
1489          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1490          else {
1491             return kFALSE;}
1492       }
1493
1494    }
1495    return kTRUE;
1496 }
1497
1498 ///________________________________________________________________________
1499 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1500 {
1501
1502    if (photon==NULL) return;
1503    Double_t facPBrem = 1.;
1504    Double_t facPSig = 0.;
1505
1506    Double_t phi=0.;
1507    Double_t theta=0.;
1508    Double_t P=0.;
1509
1510
1511    P=photon->P();
1512    phi=photon->Phi();
1513    if( photon->P()!=0){
1514       theta=acos( photon->Pz()/ photon->P());
1515    }
1516
1517    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1518       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1519    }
1520
1521    if( fPBremSmearing != 1.){
1522       if(fBrem!=NULL){
1523          facPBrem = fBrem->GetRandom();
1524       }
1525    }
1526
1527    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1528    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1529    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1530    photon->SetE(photon->P());
1531 }
1532 ///________________________________________________________________________
1533 void AliConversionMesonCuts::SmearVirtualPhoton(AliAODConversionPhoton* photon)
1534 {
1535
1536    if (photon==NULL) return;
1537    Double_t facPBrem = 1.;
1538    Double_t facPSig = 0.;
1539
1540    Double_t phi=0.;
1541    Double_t theta=0.;
1542    Double_t P=0.;
1543
1544
1545    P=photon->P();
1546    phi=photon->Phi();
1547    if( photon->P()!=0){
1548       theta=acos( photon->Pz()/ photon->P());
1549    }
1550
1551    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1552       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1553    }
1554
1555    if( fPBremSmearing != 1.){
1556       if(fBrem!=NULL){
1557          facPBrem = fBrem->GetRandom();
1558       }
1559    }
1560
1561    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1562    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1563    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1564
1565 }
1566 ///________________________________________________________________________
1567 TLorentzVector AliConversionMesonCuts::SmearElectron(TLorentzVector particle)
1568 {
1569
1570    //if (particle==0) return;
1571    Double_t facPBrem = 1.;
1572    Double_t facPSig = 0.;
1573
1574    Double_t phi=0.;
1575    Double_t theta=0.;
1576    Double_t P=0.;
1577
1578
1579    P=particle.P();
1580
1581
1582    phi=particle.Phi();
1583    if (phi < 0.) phi += 2. * TMath::Pi();
1584
1585    if( particle.P()!=0){
1586       theta=acos( particle.Pz()/ particle.P());
1587    }
1588
1589
1590    Double_t fPSigSmearingHalf    =  fPSigSmearing  / 2.0;  //The parameter was set for gammas with 2 particles and here we have just one electron
1591    Double_t sqrtfPSigSmearingCteHalf =  fPSigSmearingCte / 2.0 ;  //The parameter was set for gammas with 2 particles and here we have just one electron
1592
1593
1594
1595    if( fPSigSmearingHalf != 0. || sqrtfPSigSmearingCteHalf!=0. ){
1596       facPSig = TMath::Sqrt(sqrtfPSigSmearingCteHalf*sqrtfPSigSmearingCteHalf+fPSigSmearingHalf*fPSigSmearingHalf*P*P)*fRandom.Gaus(0.,1.);
1597    }
1598
1599    if( fPBremSmearing != 1.){
1600       if(fBrem!=NULL){
1601          facPBrem = fBrem->GetRandom();
1602       }
1603    }
1604
1605    TLorentzVector SmearedParticle;
1606
1607    SmearedParticle.SetXYZM( facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) , facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)  ,
1608                           facPBrem* (1+facPSig)* P*cos(theta) , TDatabasePDG::Instance()->GetParticle(  ::kElectron   )->Mass()) ;
1609
1610    //particle.SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1611    //particle.SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1612    //particle.SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1613
1614    return SmearedParticle;
1615
1616 }
1617