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 //      cout << "\n"<< fMCMother->GetPdgCode() << "\n" << endl;
 486         for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 487
 488                 TParticle* temp = (TParticle*)fMCStack->Particle( index );
 489 //              cout << temp->GetPdgCode() << endl;
 490                 switch( temp->GetPdgCode() ) {
 491                 case 211:
 492                         posPion      =  temp;
 493                         labelPosPion = index;
 494                         break;
 495                 case -211:
 496                         negPion      =  temp;
 497                         labelNegPion = index;
 498                         break;
 499                 case 111:
 500                         neutPion         =  temp;
 501                         labelNeutPion    = index;
 502                         break;
 503                 }
 504         }
 505
 506         if( posPion && negPion && neutPion ) return kTRUE;
 507         return kFALSE;
 508 }
 509
 510
 511 //________________________________________________________________________
 512 Bool_t AliConversionMesonCuts::MesonIsSelectedMCChiC(TParticle *fMCMother,AliStack *fMCStack,Int_t & labelelectronChiC, Int_t & labelpositronChiC, Int_t & labelgammaChiC, Double_t fRapidityShift){
 513    // Returns true for all ChiC within acceptance cuts for decay into JPsi + gamma -> e+ + e- + gamma
 514    // If bMCDaughtersInAcceptance is selected, it requires in addition that both daughter photons are within acceptance cuts
 515
 516    if(!fMCStack)return kFALSE;
 517    // if(fMCMother->GetPdgCode()==20443 ){
 518    //    return kFALSE;
 519    // }
 520    if(fMCMother->GetPdgCode()==10441 || fMCMother->GetPdgCode()==10443 || fMCMother->GetPdgCode()==445 ){
 521       if(fMCMother->R()>fMaxR)  return kFALSE; // cuts on distance from collision point
 522
 523       Double_t rapidity = 10.;
 524       if(fMCMother->Energy() - fMCMother->Pz() == 0 || fMCMother->Energy() + fMCMother->Pz() == 0){
 525          rapidity=8.-fRapidityShift;
 526       }
 527       else{
 528          rapidity = 0.5*(TMath::Log((fMCMother->Energy()+fMCMother->Pz()) / (fMCMother->Energy()-fMCMother->Pz())))-fRapidityShift;
 529       }
 530
 531       // Rapidity Cut
 532       if(abs(rapidity)>fRapidityCutMeson)return kFALSE;
 533
 534       // Select only -> ChiC radiative (JPsi+gamma) decay channel
 535       if(fMCMother->GetNDaughters()!=2)return kFALSE;
 536
 537       TParticle *jpsi   = 0x0;
 538       TParticle *gamma  = 0x0;
 539       TParticle *positron = 0x0;
 540       TParticle *electron = 0x0;
 541
 542       Int_t labeljpsiChiC = -1;
 543
 544       for(Int_t index= fMCMother->GetFirstDaughter();index<= fMCMother->GetLastDaughter();index++){
 545
 546          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 547
 548          switch( temp->GetPdgCode() ) {
 549          case 443:
 550             jpsi =  temp;
 551             labeljpsiChiC = index;
 552             break;
 553          case 22:
 554             gamma    =  temp;
 555             labelgammaChiC = index;
 556             break;
 557          }
 558       }
 559
 560       if ( !jpsi || ! gamma) return kFALSE;
 561       if(jpsi->GetNDaughters()!=2)return kFALSE;
 562
 563
 564       for(Int_t index= jpsi->GetFirstDaughter();index<= jpsi->GetLastDaughter();index++){
 565          TParticle* temp = (TParticle*)fMCStack->Particle( index );
 566          switch( temp->GetPdgCode() ) {
 567          case -11:
 568             electron =  temp;
 569             labelelectronChiC = index;
 570             break;
 571          case 11:
 572             positron =  temp;
 573             labelpositronChiC = index;
 574             break;
 575          }
 576       }
 577       if( !electron || !positron) return kFALSE;
 578       if( positron && electron && gamma) return kTRUE;
 579    }
 580    return kFALSE;
 581 }
 582
 583 ///________________________________________________________________________
 584 Bool_t AliConversionMesonCuts::MesonIsSelected(AliAODConversionMother *pi0,Bool_t IsSignal, Double_t fRapidityShift)
 585 {
 586
 587    // Selection of reconstructed Meson candidates
 588    // Use flag IsSignal in order to fill Fill different
 589    // histograms for Signal and Background
 590    TH1 *hist=0x0;
 591
 592    if(IsSignal){hist=hMesonCuts;}
 593    else{hist=hMesonBGCuts;}
 594
 595    Int_t cutIndex=0;
 596
 597    if(hist)hist->Fill(cutIndex);
 598    cutIndex++;
 599
 600    // Undefined Rapidity -> Floating Point exception
 601    if((pi0->E()+pi0->Pz())/(pi0->E()-pi0->Pz())<=0){
 602       if(hist)hist->Fill(cutIndex);
 603       cutIndex++;
 604 //        cout << "undefined rapidity" << endl;
 605       return kFALSE;
 606    }
 607    else{
 608       // PseudoRapidity Cut --> But we cut on Rapidity !!!
 609       cutIndex++;
 610       if(abs(pi0->Rapidity()-fRapidityShift)>fRapidityCutMeson){
 611          if(hist)hist->Fill(cutIndex);
 612 //               cout << abs(pi0->Rapidity()-fRapidityShift) << ">" << fRapidityCutMeson << endl;
 613          return kFALSE;
 614       }
 615    }
 616    cutIndex++;
 617
 618    // Opening Angle Cut
 619    //fOpeningAngle=2*TMath::ATan(0.134/pi0->P());// physical minimum opening angle
 620    if( pi0->GetOpeningAngle() < fOpeningAngle){
 621 //        cout << pi0->GetOpeningAngle() << "<" << fOpeningAngle << endl;
 622       if(hist)hist->Fill(cutIndex);
 623       return kFALSE;
 624    }
 625    cutIndex++;
 626
 627    if ( fAlphaPtDepCut == kTRUE ) {
 628
 629         fAlphaCutMeson = fFAlphaCut->Eval( pi0->Pt() );
 630    }
 631
 632
 633    // Alpha Max Cut
 634    if(pi0->GetAlpha()>fAlphaCutMeson){
 635 //         cout << pi0->GetAlpha() << ">" << fAlphaCutMeson << endl;
 636       if(hist)hist->Fill(cutIndex);
 637       return kFALSE;
 638    }
 639    cutIndex++;
 640
 641    // Alpha Min Cut
 642    if(pi0->GetAlpha()<fAlphaMinCutMeson){
 643 //        cout << pi0->GetAlpha() << "<" << fAlphaMinCutMeson << endl;
 644       if(hist)hist->Fill(cutIndex);
 645       return kFALSE;
 646    }
 647    cutIndex++;
 648
 649    if (hDCAGammaGammaMesonBefore)hDCAGammaGammaMesonBefore->Fill(pi0->GetDCABetweenPhotons());
 650    if (hDCARMesonPrimVtxBefore)hDCARMesonPrimVtxBefore->Fill(pi0->GetDCARMotherPrimVtx());
 651
 652    if (pi0->GetDCABetweenPhotons() > fDCAGammaGammaCut){
 653 //        cout << pi0->GetDCABetweenPhotons() << ">" << fDCAGammaGammaCut << endl;
 654       if(hist)hist->Fill(cutIndex);
 655       return kFALSE;
 656    }
 657    cutIndex++;
 658
 659    if (pi0->GetDCARMotherPrimVtx() > fDCARMesonPrimVtxCut){
 660 //         cout << pi0->GetDCARMotherPrimVtx() << ">" << fDCARMesonPrimVtxCut << endl;
 661       if(hist)hist->Fill(cutIndex);
 662       return kFALSE;
 663    }
 664    cutIndex++;
 665
 666
 667    if (hDCAZMesonPrimVtxBefore)hDCAZMesonPrimVtxBefore->Fill(pi0->GetDCAZMotherPrimVtx());
 668
 669    if (abs(pi0->GetDCAZMotherPrimVtx()) > fDCAZMesonPrimVtxCut){
 670 //        cout << pi0->GetDCAZMotherPrimVtx() << ">" << fDCAZMesonPrimVtxCut << endl;
 671       if(hist)hist->Fill(cutIndex);
 672       return kFALSE;
 673    }
 674    cutIndex++;
 675
 676
 677    if (hDCAGammaGammaMesonAfter)hDCAGammaGammaMesonAfter->Fill(pi0->GetDCABetweenPhotons());
 678    if (hDCARMesonPrimVtxAfter)hDCARMesonPrimVtxAfter->Fill(pi0->GetDCARMotherPrimVtx());
 679    if (hDCAZMesonPrimVtxAfter)hDCAZMesonPrimVtxAfter->Fill(pi0->M(),pi0->GetDCAZMotherPrimVtx());
 680
 681    if(hist)hist->Fill(cutIndex);
 682    return kTRUE;
 683 }
 684
 685
 686
 687 ///________________________________________________________________________
 688 ///________________________________________________________________________
 689 Bool_t AliConversionMesonCuts::UpdateCutString() {
 690    ///Update the cut string (if it has been created yet)
 691
 692    if(fCutString && fCutString->GetString().Length() == kNCuts) {
 693       fCutString->SetString(GetCutNumber());
 694    } else {
 695       return kFALSE;
 696    }
 697    return kTRUE;
 698 }
 699
 700 ///________________________________________________________________________
 701 Bool_t AliConversionMesonCuts::InitializeCutsFromCutString(const TString analysisCutSelection ) {
 702    // Initialize Cuts from a given Cut string
 703    AliInfo(Form("Set Meson Cutnumber: %s",analysisCutSelection.Data()));
 704    if(analysisCutSelection.Length()!=kNCuts) {
 705       AliError(Form("Cut selection has the wrong length! size is %d, number of cuts is %d", analysisCutSelection.Length(), kNCuts));
 706       return kFALSE;
 707    }
 708    if(!analysisCutSelection.IsDigit()){
 709       AliError("Cut selection contains characters");
 710       return kFALSE;
 711    }
 712
 713    const char *cutSelection = analysisCutSelection.Data();
 714 #define ASSIGNARRAY(i)  fCuts[i] = cutSelection[i] - '0'
 715    for(Int_t ii=0;ii<kNCuts;ii++){
 716       ASSIGNARRAY(ii);
 717    }
 718
 719    // Set Individual Cuts
 720    for(Int_t ii=0;ii<kNCuts;ii++){
 721       if(!SetCut(cutIds(ii),fCuts[ii]))return kFALSE;
 722    }
 723
 724    PrintCutsWithValues();
 725    return kTRUE;
 726 }
 727 ///________________________________________________________________________
 728 Bool_t AliConversionMesonCuts::SetCut(cutIds cutID, const Int_t value) {
 729    ///Set individual cut ID
 730
 731    //cout << "Updating cut  " << fgkCutNames[cutID] << " (" << cutID << ") to " << value << endl;
 732    switch (cutID) {
 733    case kMesonKind:
 734       if( SetMesonKind(value)) {
 735          fCuts[kMesonKind] = value;
 736          UpdateCutString();
 737          return kTRUE;
 738       } else return kFALSE;
 739    case kSelectionCut:
 740       if( SetSelectionWindowCut(value)) {
 741          fCuts[kSelectionCut] = value;
 742          UpdateCutString();
 743          return kTRUE;
 744       } else return kFALSE;
 745    case kalphaMesonCut:
 746       if( SetAlphaMesonCut(value)) {
 747          fCuts[kalphaMesonCut] = value;
 748          UpdateCutString();
 749          return kTRUE;
 750       } else return kFALSE;
 751    case kRCut:
 752       if( SetRCut(value)) {
 753          fCuts[kRCut] = value;
 754          UpdateCutString();
 755          return kTRUE;
 756       } else return kFALSE;
 757
 758    case kRapidityMesonCut:
 759       if( SetRapidityMesonCut(value)) {
 760          fCuts[kRapidityMesonCut] = value;
 761          UpdateCutString();
 762          return kTRUE;
 763       } else return kFALSE;
 764
 765    case kBackgroundScheme:
 766       if( SetBackgroundScheme(value)) {
 767          fCuts[kBackgroundScheme] = value;
 768          UpdateCutString();
 769          return kTRUE;
 770       } else return kFALSE;
 771
 772    case kDegreesForRotationMethod:
 773       if( SetNDegreesForRotationMethod(value)) {
 774          fCuts[kDegreesForRotationMethod] = value;
 775          UpdateCutString();
 776          return kTRUE;
 777       } else return kFALSE;
 778
 779    case kNumberOfBGEvents:
 780       if( SetNumberOfBGEvents(value)) {
 781          fCuts[kNumberOfBGEvents] = value;
 782          UpdateCutString();
 783          return kTRUE;
 784       } else return kFALSE;
 785
 786    case kuseMCPSmearing:
 787       if( SetMCPSmearing(value)) {
 788          fCuts[kuseMCPSmearing] = value;
 789          UpdateCutString();
 790          return kTRUE;
 791       } else return kFALSE;
 792    case kElecShare:
 793       if( SetSharedElectronCut(value)) {
 794          fCuts[kElecShare] = value;
 795          UpdateCutString();
 796          return kTRUE;
 797       } else return kFALSE;
 798    case kToCloseV0s:
 799       if( SetToCloseV0sCut(value)) {
 800          fCuts[kToCloseV0s] = value;
 801          UpdateCutString();
 802          return kTRUE;
 803       } else return kFALSE;
 804    case kDcaGammaGamma:
 805       if( SetDCAGammaGammaCut(value)) {
 806          fCuts[kDcaGammaGamma] = value;
 807          UpdateCutString();
 808          return kTRUE;
 809       } else return kFALSE;
 810    case kDcaZPrimVtx:
 811       if( SetDCAZMesonPrimVtxCut(value)) {
 812          fCuts[kDcaZPrimVtx] = value;
 813          UpdateCutString();
 814          return kTRUE;
 815       } else return kFALSE;
 816    case kDcaRPrimVtx:
 817       if( SetDCARMesonPrimVtxCut(value)) {
 818          fCuts[kDcaRPrimVtx] = value;
 819          UpdateCutString();
 820          return kTRUE;
 821       } else return kFALSE;
 822
 823    case kNCuts:
 824       cout << "Error:: Cut id out of range"<< endl;
 825       return kFALSE;
 826    }
 827
 828    cout << "Error:: Cut id " << cutID << " not recognized "<< endl;
 829    return kFALSE;
 830
 831 }
 832
 833
 834 ///________________________________________________________________________
 835 void AliConversionMesonCuts::PrintCuts() {
 836    // Print out current Cut Selection
 837    for(Int_t ic = 0; ic < kNCuts; ic++) {
 838       printf("%-30s : %d \n", fgkCutNames[ic], fCuts[ic]);
 839    }
 840 }
 841
 842 ///________________________________________________________________________
 843 void AliConversionMesonCuts::PrintCutsWithValues() {
 844    // Print out current Cut Selection with values
 845         printf("\nMeson cutnumber \n");
 846         for(Int_t ic = 0; ic < kNCuts; ic++) {
 847                 printf("%d",fCuts[ic]);
 848         }
 849         printf("\n\n");
 850
 851         printf("Meson cuts \n");
 852         printf("\t |y| < %3.2f \n", fRapidityCutMeson);
 853         printf("\t theta_{open} < %3.2f\n", fOpeningAngle);
 854         if (!fAlphaPtDepCut) printf("\t %3.2f < alpha < %3.2f\n", fAlphaMinCutMeson, fAlphaCutMeson);
 855         printf("\t dca_{gamma,gamma} > %3.2f\n", fDCAGammaGammaCut);
 856         printf("\t dca_{R, prim Vtx} > %3.2f\n", fDCARMesonPrimVtxCut);
 857         printf("\t dca_{Z, prim Vtx} > %3.2f\n\n", fDCAZMesonPrimVtxCut);
 858         printf("\t Meson selection window for further analysis %3.3f > M_{gamma,gamma} > %3.3f\n\n", fSelectionLow, fSelectionHigh);
 859
 860          printf("Meson BG settings \n");
 861          if (!fDoBG){
 862                 if (!fUseRotationMethodInBG  & !fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing V0 mult \n");
 863                 if (!fUseRotationMethodInBG  & fUseTrackMultiplicityForBG) printf("\t BG scheme: mixing track mult \n");
 864                 if (fUseRotationMethodInBG )printf("\t BG scheme: rotation \n");
 865                 if (fdoBGProbability) printf("\t -> use BG probability \n");
 866                 if (fBackgroundHandler) printf("\t -> use new BG handler \n");
 867                 printf("\t depth of pool: %d\n", fNumberOfBGEvents);
 868                 if (fUseRotationMethodInBG )printf("\t degree's for BG rotation: %d\n", fnDegreeRotationPMForBG);
 869          }
 870
 871 }
 872
 873
 874 ///________________________________________________________________________
 875 Bool_t AliConversionMesonCuts::SetMesonKind(Int_t mesonKind){
 876    // Set Cut
 877    switch(mesonKind){
 878    case 0:
 879       fMesonKind=0;
 880       break;
 881    case 1:
 882       fMesonKind=1;;
 883       break;
 884    default:
 885       cout<<"Warning: Meson kind not defined"<<mesonKind<<endl;
 886       return kFALSE;
 887    }
 888    return kTRUE;
 889 }
 890
 891 ///________________________________________________________________________
 892 Bool_t AliConversionMesonCuts::SetRCut(Int_t rCut){
 893    // Set Cut
 894    switch(rCut){
 895    case 0:
 896       fMaxR = 180.;
 897       break;
 898    case 1:
 899       fMaxR = 180.;
 900       break;
 901    case 2:
 902       fMaxR = 180.;
 903       break;
 904    case 3:
 905       fMaxR = 70.;
 906       break;
 907    case 4:
 908       fMaxR = 70.;
 909       break;
 910    case 5:
 911       fMaxR = 180.;
 912       break;
 913       // High purity cuts for PbPb
 914    case 9:
 915       fMaxR = 180.;
 916       break;
 917    default:
 918       cout<<"Warning: rCut not defined"<<rCut<<endl;
 919       return kFALSE;
 920    }
 921    return kTRUE;
 922 }
 923
 924 ///________________________________________________________________________
 925 Bool_t AliConversionMesonCuts::SetSelectionWindowCut(Int_t selectionCut){
 926    // Set Cut
 927    switch(selectionCut){
 928    case 0:
 929       fSelectionLow = 0.08;
 930           fSelectionHigh = 0.145;
 931       break;
 932    case 1:
 933       fSelectionLow = 0.1;
 934           fSelectionHigh = 0.145;
 935       break;
 936    case 2:
 937           fSelectionLow = 0.11;
 938           fSelectionHigh = 0.145;
 939           break;
 940    case 3:
 941           fSelectionLow = 0.12;
 942           fSelectionHigh = 0.145;
 943       break;
 944    case 4:
 945           fSelectionLow = 0.1;
 946           fSelectionHigh = 0.15;
 947       break;
 948    case 5:
 949           fSelectionLow = 0.11;
 950           fSelectionHigh = 0.15;
 951       break;
 952    case 6:
 953           fSelectionLow = 0.12;
 954           fSelectionHigh = 0.15;
 955       break;
 956
 957    default:
 958       cout<<"Warning: SelectionCut not defined "<<selectionCut<<endl;
 959       return kFALSE;
 960    }
 961    return kTRUE;
 962 }
 963
 964
 965 ///________________________________________________________________________
 966 Bool_t AliConversionMesonCuts::SetAlphaMesonCut(Int_t alphaMesonCut)
 967 {   // Set Cut
 968    switch(alphaMesonCut){
 969    case 0:      // 0- 0.7
 970       fAlphaMinCutMeson  = 0.0;
 971       fAlphaCutMeson     = 0.7;
 972       fAlphaPtDepCut = kFALSE;
 973       break;
 974    case 1:      // Updated 31 October 2013 before 0.0 - 0.5
 975       if( fFAlphaCut ) delete fFAlphaCut;
 976       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 977       fFAlphaCut->SetParameter(0,0.7);
 978       fFAlphaCut->SetParameter(1,1.2);
 979       fAlphaMinCutMeson  =  0.0;
 980       fAlphaCutMeson     = -1.0;
 981       fAlphaPtDepCut = kTRUE;
 982       break;
 983    case 2:      // Updated 31 October 2013 before 0.5-1
 984       if( fFAlphaCut ) delete fFAlphaCut;
 985       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
 986       fFAlphaCut->SetParameter(0,0.8);
 987       fFAlphaCut->SetParameter(1,1.2);
 988       fAlphaMinCutMeson  =  0.0;
 989       fAlphaCutMeson     = -1.0;
 990       fAlphaPtDepCut = kTRUE;
 991       break;
 992    case 3:      // 0.0-1
 993       fAlphaMinCutMeson  = 0.0;
 994       fAlphaCutMeson     = 1.;
 995       fAlphaPtDepCut = kFALSE;
 996       break;
 997    case 4:      // 0-0.65
 998       fAlphaMinCutMeson  = 0.0;
 999       fAlphaCutMeson     = 0.65;
1000       fAlphaPtDepCut = kFALSE;
1001       break;
1002    case 5:      // 0-0.75
1003       fAlphaMinCutMeson  = 0.0;
1004       fAlphaCutMeson     = 0.75;
1005       fAlphaPtDepCut = kFALSE;
1006       break;
1007    case 6:      // 0-0.8
1008       fAlphaMinCutMeson  = 0.0;
1009       fAlphaCutMeson     = 0.8;
1010       fAlphaPtDepCut = kFALSE;
1011       break;
1012    case 7:      // 0.0-0.85
1013       fAlphaMinCutMeson  = 0.0;
1014       fAlphaCutMeson     = 0.85;
1015       fAlphaPtDepCut = kFALSE;
1016       break;
1017    case 8:      // 0.0-0.6
1018       fAlphaMinCutMeson  = 0.0;
1019       fAlphaCutMeson     = 0.6;
1020       fAlphaPtDepCut = kFALSE;
1021       break;
1022    case 9: // Updated 11 November 2013 before 0.0 - 0.3
1023       if( fFAlphaCut ) delete fFAlphaCut;
1024       fFAlphaCut= new TF1("fFAlphaCut","[0]*tanh([1]*x)",0.,100.);
1025       fFAlphaCut->SetParameter(0,0.65);
1026       fFAlphaCut->SetParameter(1,1.2);
1027       fAlphaMinCutMeson  =  0.0;
1028       fAlphaCutMeson     = -1.0;
1029       fAlphaPtDepCut = kTRUE;
1030       break;
1031    default:
1032       cout<<"Warning: AlphaMesonCut not defined "<<alphaMesonCut<<endl;
1033       return kFALSE;
1034    }
1035    return kTRUE;
1036 }
1037
1038 ///________________________________________________________________________
1039 Bool_t AliConversionMesonCuts::SetRapidityMesonCut(Int_t RapidityMesonCut){
1040    // Set Cut
1041    switch(RapidityMesonCut){
1042    case 0:  // changed from 0.9 to 1.35
1043       fRapidityCutMeson   = 1.35;
1044       break;
1045    case 1:  //
1046       fRapidityCutMeson   = 0.8;
1047       break;
1048    case 2:  //
1049       fRapidityCutMeson   = 0.7;
1050       break;
1051    case 3:  //
1052       fRapidityCutMeson   = 0.6;
1053       break;
1054    case 4:  //
1055       fRapidityCutMeson   = 0.5;
1056       break;
1057    case 5:  //
1058       fRapidityCutMeson   = 0.85;
1059       break;
1060    case 6:  //
1061       fRapidityCutMeson   = 0.75;
1062       break;
1063    case 7:  //
1064       fRapidityCutMeson   = 0.3;
1065       break;
1066    case 8:  //changed, before 0.35
1067       fRapidityCutMeson   = 0.25;
1068       break;
1069    case 9:  //
1070       fRapidityCutMeson   = 0.4;
1071       break;
1072    default:
1073       cout<<"Warning: RapidityMesonCut not defined "<<RapidityMesonCut<<endl;
1074       return kFALSE;
1075    }
1076    return kTRUE;
1077 }
1078
1079
1080 ///________________________________________________________________________
1081 Bool_t AliConversionMesonCuts::SetBackgroundScheme(Int_t BackgroundScheme){
1082    // Set Cut
1083    switch(BackgroundScheme){
1084    case 0: //Rotation
1085       fUseRotationMethodInBG=kTRUE;
1086       fdoBGProbability=kFALSE;
1087       break;
1088    case 1: // mixed event with V0 multiplicity
1089       fUseRotationMethodInBG=kFALSE;
1090       fUseTrackMultiplicityForBG=kFALSE;
1091       fdoBGProbability=kFALSE;
1092       break;
1093    case 2: // mixed event with track multiplicity
1094       fUseRotationMethodInBG=kFALSE;
1095       fUseTrackMultiplicityForBG=kTRUE;
1096       fdoBGProbability=kFALSE;
1097       break;
1098    case 3: //Rotation
1099       fUseRotationMethodInBG=kTRUE;
1100       fdoBGProbability=kTRUE;
1101       break;
1102    case 4: //No BG calculation
1103       cout << "no BG calculation should be done" << endl;
1104       fUseRotationMethodInBG=kFALSE;
1105       fdoBGProbability=kFALSE;
1106       fDoBG=kFALSE;
1107       break;
1108    case 5: //Rotation
1109       fUseRotationMethodInBG=kTRUE;
1110       fdoBGProbability=kFALSE;
1111       fBackgroundHandler = 1;
1112       break;
1113    case 6: // mixed event with V0 multiplicity
1114       fUseRotationMethodInBG=kFALSE;
1115       fUseTrackMultiplicityForBG=kFALSE;
1116       fdoBGProbability=kFALSE;
1117       fBackgroundHandler = 1;
1118       break;
1119    case 7: // mixed event with track multiplicity
1120       fUseRotationMethodInBG=kFALSE;
1121       fUseTrackMultiplicityForBG=kTRUE;
1122       fdoBGProbability=kFALSE;
1123       fBackgroundHandler = 1;
1124       break;
1125    case 8: //Rotation
1126       fUseRotationMethodInBG=kTRUE;
1127       fdoBGProbability=kTRUE;
1128       fBackgroundHandler = 1;
1129       break;
1130    default:
1131       cout<<"Warning: BackgroundScheme not defined "<<BackgroundScheme<<endl;
1132       return kFALSE;
1133    }
1134    return kTRUE;
1135 }
1136
1137
1138 ///________________________________________________________________________
1139 Bool_t AliConversionMesonCuts::SetNDegreesForRotationMethod(Int_t DegreesForRotationMethod){
1140    // Set Cut
1141    switch(DegreesForRotationMethod){
1142    case 0:
1143       fnDegreeRotationPMForBG = 5;
1144       break;
1145    case 1:
1146       fnDegreeRotationPMForBG = 10;
1147       break;
1148    case 2:
1149       fnDegreeRotationPMForBG = 15;
1150       break;
1151    case 3:
1152       fnDegreeRotationPMForBG = 20;
1153       break;
1154    default:
1155       cout<<"Warning: DegreesForRotationMethod not defined "<<DegreesForRotationMethod<<endl;
1156       return kFALSE;
1157    }
1158    fCuts[kDegreesForRotationMethod]=DegreesForRotationMethod;
1159    return kTRUE;
1160 }
1161
1162 ///________________________________________________________________________
1163 Bool_t AliConversionMesonCuts::SetNumberOfBGEvents(Int_t NumberOfBGEvents){
1164    // Set Cut
1165    switch(NumberOfBGEvents){
1166    case 0:
1167       fNumberOfBGEvents = 5;
1168       break;
1169    case 1:
1170       fNumberOfBGEvents = 10;
1171       break;
1172    case 2:
1173       fNumberOfBGEvents = 15;
1174       break;
1175    case 3:
1176       fNumberOfBGEvents = 20;
1177       break;
1178    case 4:
1179       fNumberOfBGEvents = 2;
1180       break;
1181    case 5:
1182       fNumberOfBGEvents = 50;
1183       break;
1184    case 6:
1185       fNumberOfBGEvents = 80;
1186       break;
1187    case 7:
1188       fNumberOfBGEvents = 100;
1189       break;
1190    default:
1191       cout<<"Warning: NumberOfBGEvents not defined "<<NumberOfBGEvents<<endl;
1192       return kFALSE;
1193    }
1194    return kTRUE;
1195 }
1196 ///________________________________________________________________________
1197 Bool_t AliConversionMesonCuts::SetSharedElectronCut(Int_t sharedElec) {
1198
1199    switch(sharedElec){
1200    case 0:
1201       fDoSharedElecCut = kFALSE;
1202       break;
1203    case 1:
1204       fDoSharedElecCut = kTRUE;
1205       break;
1206    default:
1207       cout<<"Warning: Shared Electron Cut not defined "<<sharedElec<<endl;
1208       return kFALSE;
1209    }
1210
1211    return kTRUE;
1212 }
1213
1214 ///________________________________________________________________________
1215 Bool_t AliConversionMesonCuts::SetToCloseV0sCut(Int_t toClose) {
1216
1217    switch(toClose){
1218    case 0:
1219       fDoToCloseV0sCut = kFALSE;
1220       fminV0Dist = 250;
1221       break;
1222    case 1:
1223       fDoToCloseV0sCut = kTRUE;
1224       fminV0Dist = 1;
1225       break;
1226    case 2:
1227       fDoToCloseV0sCut = kTRUE;
1228       fminV0Dist = 2;
1229       break;
1230    case 3:
1231       fDoToCloseV0sCut = kTRUE;
1232       fminV0Dist = 3;
1233       break;
1234    default:
1235       cout<<"Warning: Shared Electron Cut not defined "<<toClose<<endl;
1236       return kFALSE;
1237    }
1238    return kTRUE;
1239 }
1240
1241 ///________________________________________________________________________
1242 Bool_t AliConversionMesonCuts::SetMCPSmearing(Int_t useMCPSmearing)
1243 {// Set Cut
1244    switch(useMCPSmearing){
1245    case 0:
1246       fUseMCPSmearing=0;
1247       fPBremSmearing=1.;
1248       fPSigSmearing=0.;
1249       fPSigSmearingCte=0.;
1250       break;
1251    case 1:
1252       fUseMCPSmearing=1;
1253       fPBremSmearing=1.0e-14;
1254       fPSigSmearing=0.;
1255       fPSigSmearingCte=0.;
1256       break;
1257    case 2:
1258       fUseMCPSmearing=1;
1259       fPBremSmearing=1.0e-15;
1260       fPSigSmearing=0.0;
1261       fPSigSmearingCte=0.;
1262       break;
1263    case 3:
1264       fUseMCPSmearing=1;
1265       fPBremSmearing=1.;
1266       fPSigSmearing=0.003;
1267       fPSigSmearingCte=0.002;
1268       break;
1269    case 4:
1270       fUseMCPSmearing=1;
1271       fPBremSmearing=1.;
1272       fPSigSmearing=0.003;
1273       fPSigSmearingCte=0.007;
1274       break;
1275    case 5:
1276       fUseMCPSmearing=1;
1277       fPBremSmearing=1.;
1278       fPSigSmearing=0.003;
1279       fPSigSmearingCte=0.016;
1280       break;
1281    case 6:
1282       fUseMCPSmearing=1;
1283       fPBremSmearing=1.;
1284       fPSigSmearing=0.007;
1285       fPSigSmearingCte=0.016;
1286       break;
1287    case 7:
1288       fUseMCPSmearing=1;
1289       fPBremSmearing=1.0e-16;
1290       fPSigSmearing=0.0;
1291       fPSigSmearingCte=0.;
1292       break;
1293    case 8:
1294       fUseMCPSmearing=1;
1295       fPBremSmearing=1.;
1296       fPSigSmearing=0.007;
1297       fPSigSmearingCte=0.014;
1298       break;
1299    case 9:
1300       fUseMCPSmearing=1;
1301       fPBremSmearing=1.;
1302       fPSigSmearing=0.007;
1303       fPSigSmearingCte=0.011;
1304       break;
1305
1306    default:
1307       AliError("Warning: UseMCPSmearing not defined");
1308       return kFALSE;
1309    }
1310    return kTRUE;
1311 }
1312
1313
1314 ///________________________________________________________________________
1315 Bool_t AliConversionMesonCuts::SetDCAGammaGammaCut(Int_t DCAGammaGamma){
1316    // Set Cut
1317    switch(DCAGammaGamma){
1318    case 0:  //
1319       fDCAGammaGammaCut   = 1000;
1320       break;
1321    case 1:  //
1322       fDCAGammaGammaCut   = 10;
1323       break;
1324    case 2:  //
1325       fDCAGammaGammaCut   = 5;
1326       break;
1327    case 3:  //
1328       fDCAGammaGammaCut   = 4;
1329       break;
1330    case 4:  //
1331       fDCAGammaGammaCut   = 3;
1332       break;
1333    case 5:  //
1334       fDCAGammaGammaCut   = 2.5;
1335       break;
1336    case 6:  //
1337       fDCAGammaGammaCut   = 2;
1338       break;
1339    case 7:  //
1340       fDCAGammaGammaCut   = 1.5;
1341       break;
1342    case 8:  //
1343       fDCAGammaGammaCut   = 1;
1344       break;
1345    case 9:  //
1346       fDCAGammaGammaCut   = 0.5;
1347       break;
1348    default:
1349       cout<<"Warning: DCAGammaGamma not defined "<<DCAGammaGamma<<endl;
1350       return kFALSE;
1351    }
1352    return kTRUE;
1353 }
1354
1355 ///________________________________________________________________________
1356 Bool_t AliConversionMesonCuts::SetDCAZMesonPrimVtxCut(Int_t DCAZMesonPrimVtx){
1357    // Set Cut
1358    switch(DCAZMesonPrimVtx){
1359    case 0:  //
1360       fDCAZMesonPrimVtxCut   = 1000;
1361       break;
1362    case 1:  //
1363       fDCAZMesonPrimVtxCut   = 10;
1364       break;
1365    case 2:  //
1366       fDCAZMesonPrimVtxCut   = 5;
1367       break;
1368    case 3:  //
1369       fDCAZMesonPrimVtxCut   = 4;
1370       break;
1371    case 4:  //
1372       fDCAZMesonPrimVtxCut   = 3;
1373       break;
1374    case 5:  //
1375       fDCAZMesonPrimVtxCut   = 2.5;
1376       break;
1377    case 6:  //
1378       fDCAZMesonPrimVtxCut   = 2;
1379       break;
1380    case 7:  //
1381       fDCAZMesonPrimVtxCut   = 1.5;
1382       break;
1383    case 8:  //
1384       fDCAZMesonPrimVtxCut   = 1;
1385       break;
1386    case 9:  //
1387       fDCAZMesonPrimVtxCut   = 0.5;
1388       break;
1389    default:
1390       cout<<"Warning: DCAZMesonPrimVtx not defined "<<DCAZMesonPrimVtx<<endl;
1391       return kFALSE;
1392    }
1393    return kTRUE;
1394 }
1395
1396 ///________________________________________________________________________
1397 Bool_t AliConversionMesonCuts::SetDCARMesonPrimVtxCut(Int_t DCARMesonPrimVtx){
1398    // Set Cut
1399    switch(DCARMesonPrimVtx){
1400    case 0:  //
1401       fDCARMesonPrimVtxCut   = 1000;
1402       break;
1403    case 1:  //
1404       fDCARMesonPrimVtxCut   = 10;
1405       break;
1406    case 2:  //
1407       fDCARMesonPrimVtxCut   = 5;
1408       break;
1409    case 3:  //
1410       fDCARMesonPrimVtxCut   = 4;
1411       break;
1412    case 4:  //
1413       fDCARMesonPrimVtxCut   = 3;
1414       break;
1415    case 5:  //
1416       fDCARMesonPrimVtxCut   = 2.5;
1417       break;
1418    case 6:  //
1419       fDCARMesonPrimVtxCut   = 2;
1420       break;
1421    case 7:  //
1422       fDCARMesonPrimVtxCut   = 1.5;
1423       break;
1424    case 8:  //
1425       fDCARMesonPrimVtxCut   = 1;
1426       break;
1427    case 9:  //
1428       fDCARMesonPrimVtxCut   = 0.5;
1429       break;
1430    default:
1431       cout<<"Warning: DCARMesonPrimVtx not defined "<<DCARMesonPrimVtx<<endl;
1432       return kFALSE;
1433    }
1434    return kTRUE;
1435 }
1436
1437
1438 ///________________________________________________________________________
1439 TString AliConversionMesonCuts::GetCutNumber(){
1440    // returns TString with current cut number
1441    TString a(kNCuts);
1442    for(Int_t ii=0;ii<kNCuts;ii++){
1443       a.Append(Form("%d",fCuts[ii]));
1444    }
1445    return a;
1446 }
1447
1448 ///________________________________________________________________________
1449 void AliConversionMesonCuts::FillElectonLabelArray(AliAODConversionPhoton* photon, Int_t nV0){
1450
1451    Int_t posLabel = photon->GetTrackLabelPositive();
1452    Int_t negLabel = photon->GetTrackLabelNegative();
1453
1454    fElectronLabelArray[nV0*2] = posLabel;
1455    fElectronLabelArray[(nV0*2)+1] = negLabel;
1456 }
1457
1458 ///________________________________________________________________________
1459 Bool_t AliConversionMesonCuts::RejectSharedElectronV0s(AliAODConversionPhoton* photon, Int_t nV0, Int_t nV0s){
1460
1461    Int_t posLabel = photon->GetTrackLabelPositive();
1462    Int_t negLabel = photon->GetTrackLabelNegative();
1463
1464    for(Int_t i = 0; i<nV0s*2;i++){
1465       if(i==nV0*2)     continue;
1466       if(i==(nV0*2)+1) continue;
1467       if(fElectronLabelArray[i] == posLabel){
1468          return kFALSE;}
1469       if(fElectronLabelArray[i] == negLabel){
1470          return kFALSE;}
1471    }
1472
1473    return kTRUE;
1474 }
1475 ///________________________________________________________________________
1476 Bool_t AliConversionMesonCuts::RejectToCloseV0s(AliAODConversionPhoton* photon, TList *photons, Int_t nV0){
1477    Double_t posX = photon->GetConversionX();
1478    Double_t posY = photon->GetConversionY();
1479    Double_t posZ = photon->GetConversionZ();
1480
1481    for(Int_t i = 0;i<photons->GetEntries();i++){
1482       if(nV0 == i) continue;
1483       AliAODConversionPhoton *photonComp = (AliAODConversionPhoton*) photons->At(i);
1484       Double_t posCompX = photonComp->GetConversionX();
1485       Double_t posCompY = photonComp->GetConversionY();
1486       Double_t posCompZ = photonComp->GetConversionZ();
1487
1488       Double_t dist = pow((posX - posCompX),2)+pow((posY - posCompY),2)+pow((posZ - posCompZ),2);
1489
1490       if(dist < fminV0Dist*fminV0Dist){
1491          if(photon->GetChi2perNDF() < photonComp->GetChi2perNDF()) return kTRUE;
1492          else {
1493             return kFALSE;}
1494       }
1495
1496    }
1497    return kTRUE;
1498 }
1499
1500 ///________________________________________________________________________
1501 void AliConversionMesonCuts::SmearParticle(AliAODConversionPhoton* photon)
1502 {
1503
1504    if (photon==NULL) return;
1505    Double_t facPBrem = 1.;
1506    Double_t facPSig = 0.;
1507
1508    Double_t phi=0.;
1509    Double_t theta=0.;
1510    Double_t P=0.;
1511
1512
1513    P=photon->P();
1514    phi=photon->Phi();
1515    if( photon->P()!=0){
1516       theta=acos( photon->Pz()/ photon->P());
1517    }
1518
1519    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1520       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1521    }
1522
1523    if( fPBremSmearing != 1.){
1524       if(fBrem!=NULL){
1525          facPBrem = fBrem->GetRandom();
1526       }
1527    }
1528
1529    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1530    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1531    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1532    photon->SetE(photon->P());
1533 }
1534 ///________________________________________________________________________
1535 void AliConversionMesonCuts::SmearVirtualPhoton(AliAODConversionPhoton* photon)
1536 {
1537
1538    if (photon==NULL) return;
1539    Double_t facPBrem = 1.;
1540    Double_t facPSig = 0.;
1541
1542    Double_t phi=0.;
1543    Double_t theta=0.;
1544    Double_t P=0.;
1545
1546
1547    P=photon->P();
1548    phi=photon->Phi();
1549    if( photon->P()!=0){
1550       theta=acos( photon->Pz()/ photon->P());
1551    }
1552
1553    if( fPSigSmearing != 0. || fPSigSmearingCte!=0. ){
1554       facPSig = TMath::Sqrt(fPSigSmearingCte*fPSigSmearingCte+fPSigSmearing*fPSigSmearing*P*P)*fRandom.Gaus(0.,1.);
1555    }
1556
1557    if( fPBremSmearing != 1.){
1558       if(fBrem!=NULL){
1559          facPBrem = fBrem->GetRandom();
1560       }
1561    }
1562
1563    photon->SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1564    photon->SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1565    photon->SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1566
1567 }
1568 ///________________________________________________________________________
1569 TLorentzVector AliConversionMesonCuts::SmearElectron(TLorentzVector particle)
1570 {
1571
1572    //if (particle==0) return;
1573    Double_t facPBrem = 1.;
1574    Double_t facPSig = 0.;
1575
1576    Double_t phi=0.;
1577    Double_t theta=0.;
1578    Double_t P=0.;
1579
1580
1581    P=particle.P();
1582
1583
1584    phi=particle.Phi();
1585    if (phi < 0.) phi += 2. * TMath::Pi();
1586
1587    if( particle.P()!=0){
1588       theta=acos( particle.Pz()/ particle.P());
1589    }
1590
1591
1592    Double_t fPSigSmearingHalf    =  fPSigSmearing  / 2.0;  //The parameter was set for gammas with 2 particles and here we have just one electron
1593    Double_t sqrtfPSigSmearingCteHalf =  fPSigSmearingCte / 2.0 ;  //The parameter was set for gammas with 2 particles and here we have just one electron
1594
1595
1596
1597    if( fPSigSmearingHalf != 0. || sqrtfPSigSmearingCteHalf!=0. ){
1598       facPSig = TMath::Sqrt(sqrtfPSigSmearingCteHalf*sqrtfPSigSmearingCteHalf+fPSigSmearingHalf*fPSigSmearingHalf*P*P)*fRandom.Gaus(0.,1.);
1599    }
1600
1601    if( fPBremSmearing != 1.){
1602       if(fBrem!=NULL){
1603          facPBrem = fBrem->GetRandom();
1604       }
1605    }
1606
1607    TLorentzVector SmearedParticle;
1608
1609    SmearedParticle.SetXYZM( facPBrem* (1+facPSig)* P*sin(theta)*cos(phi) , facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)  ,
1610                           facPBrem* (1+facPSig)* P*cos(theta) , TDatabasePDG::Instance()->GetParticle(  ::kElectron   )->Mass()) ;
1611
1612    //particle.SetPx(facPBrem* (1+facPSig)* P*sin(theta)*cos(phi)) ;
1613    //particle.SetPy(facPBrem* (1+facPSig)* P*sin(theta)*sin(phi)) ;
1614    //particle.SetPz(facPBrem* (1+facPSig)* P*cos(theta)) ;
1615
1616    return SmearedParticle;
1617
1618 }
1619