PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /////////////////////////////////////////////////////////////
  18 //
  19 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
  20 // and comparison with the MC truth and cut variables distributions.
  21 //
  22 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
  23 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
  24 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
  25 /////////////////////////////////////////////////////////////
  26
  27 #include <Riostream.h>
  28 #include <TClonesArray.h>
  29 #include <TCanvas.h>
  30 #include <TNtuple.h>
  31 #include <TList.h>
  32 #include <TH1F.h>
  33 #include <TH2F.h>
  34 #include <TDatabasePDG.h>
  35
  36 #include "AliAnalysisManager.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODMCHeader.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODRecoDecayHF2Prong.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSED0Mass.h"
  48
  49
  50 ClassImp(AliAnalysisTaskSED0Mass)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  55 AliAnalysisTaskSE(),
  56 fOutputPPR(0),
  57 fOutputmycuts(0),
  58 fDistr(0),
  59 fNentries(0),
  60 fVHFPPR(0),
  61 fVHFmycuts(0),
  62 fArray(0),
  63 fReadMC(0),
  64 fLsNormalization(1.)
  65
  66 {
  67   // Default constructor
  68   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  69 }
  70
  71 //________________________________________________________________________
  72 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
  73 AliAnalysisTaskSE(name),
  74 fOutputPPR(0),
  75 fOutputmycuts(0),
  76 fDistr(0),
  77 fNentries(0),
  78 fVHFPPR(0),
  79 fVHFmycuts(0),
  80 fArray(0),
  81 fReadMC(0),
  82 fLsNormalization(1.)
  83
  84 {
  85   // Default constructor
  86   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  87
  88   // Output slot #1 writes into a TList container
  89   DefineOutput(1,TList::Class());  //My private output
  90   // Output slot #2 writes into a TList container
  91   DefineOutput(2,TList::Class());  //My private output
  92   // Output slot #3 writes into a TH1F container
  93   DefineOutput(3,TH1F::Class());  //My private output
  94   // Output slot #4 writes into a TList container
  95   DefineOutput(4,TList::Class());  //My private output
  96 }
  97
  98 //________________________________________________________________________
  99 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
 100 {
 101    if (fOutputPPR) {
 102     delete fOutputPPR;
 103     fOutputPPR = 0;
 104   }
 105   if (fVHFPPR) {
 106     delete fVHFPPR;
 107     fVHFPPR = 0;
 108   }
 109   if (fOutputmycuts) {
 110     delete fOutputmycuts;
 111     fOutputmycuts = 0;
 112   }
 113
 114   if (fDistr) {
 115     delete fDistr;
 116     fDistr = 0;
 117   }
 118
 119   if (fVHFmycuts) {
 120     delete fVHFmycuts;
 121     fVHFmycuts = 0;
 122   }
 123   if (fNentries){
 124     delete fNentries;
 125     fNentries = 0;
 126   }
 127
 128
 129 }
 130
 131 //________________________________________________________________________
 132 void AliAnalysisTaskSED0Mass::Init()
 133 {
 134   // Initialization
 135
 136   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 137
 138   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
 139
 140   // 2 sets of dedidcated cuts -- defined in UserExec
 141   fVHFPPR = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 142   fVHFmycuts = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 143
 144   return;
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 149 {
 150
 151   // Create the output container
 152   //
 153   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 154
 155   // Several histograms are more conveniently managed in a TList
 156   fOutputPPR = new TList();
 157   fOutputPPR->SetOwner();
 158   fOutputPPR->SetName("listPPR");
 159
 160   fOutputmycuts = new TList();
 161   fOutputmycuts->SetOwner();
 162   fOutputmycuts->SetName("listloose");
 163
 164   fDistr = new TList();
 165   fDistr->SetOwner();
 166   fDistr->SetName("distributionslist");
 167
 168   const Int_t nhist=5;
 169
 170   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
 171
 172   for(Int_t i=0;i<nhist;i++){
 173
 174     nameMass="histMass_";
 175     nameMass+=i+1;
 176     nameSgn27="histSgn27_";
 177     nameSgn27+=i+1;
 178     nameSgn="histSgn_";
 179     nameSgn+=i+1;
 180     nameBkg="histBkg_";
 181     nameBkg+=i+1;
 182     nameRfl="histRfl_";
 183     nameRfl+=i+1;
 184     nameMassNocutsS="hMassS_";
 185     nameMassNocutsS+=i+1;
 186     nameMassNocutsB="hMassB_";
 187     nameMassNocutsB+=i+1;
 188
 189     //histograms of cut variable distributions
 190
 191     //  pT
 192     namedistr="hptpiS_";
 193     namedistr+=i+1;
 194     TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 195
 196     namedistr="hptKS_";
 197     namedistr+=i+1;
 198     TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 199
 200     namedistr="hptB_";
 201     namedistr+=i+1;
 202     TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 203
 204     //  pT no mass cut
 205     namedistr="hptpiSnoMcut_";
 206     namedistr+=i+1;
 207     TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 208
 209     namedistr="hptKSnoMcut_";
 210     namedistr+=i+1;
 211     TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 212
 213     namedistr="hptB1prongnoMcut_";
 214     namedistr+=i+1;
 215     TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 216
 217     namedistr="hptB2prongsnoMcut_";
 218     namedistr+=i+1;
 219     TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 220
 221
 222     //  dca
 223     namedistr="hdcaS_";
 224     namedistr+=i+1;
 225     TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 226     namedistr="hdcaB_";
 227     namedistr+=i+1;
 228     TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 229
 230     //  costhetastar
 231     namedistr="hcosthetastarS_";
 232     namedistr+=i+1;
 233     TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 234     namedistr="hcosthetastarB_";
 235     namedistr+=i+1;
 236     TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 237
 238     // impact parameter
 239     namedistr="hd0piS_";
 240     namedistr+=i+1;
 241     TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
 242
 243     namedistr="hd0KS_";
 244     namedistr+=i+1;
 245     TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
 246     namedistr="hd0B_";
 247     namedistr+=i+1;
 248     TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution;d0 [cm]",200,-0.1,0.1);
 249
 250     namedistr="hd0d0S_";
 251     namedistr+=i+1;
 252     TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 253     namedistr="hd0d0B_";
 254     namedistr+=i+1;
 255     TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 256
 257     //  costhetapoint
 258     namedistr="hcosthetapointS_";
 259     namedistr+=i+1;
 260     TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 261     namedistr="hcosthetapointB_";
 262     namedistr+=i+1;
 263     TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 264
 265     namedistr="hcosthpointd0d0S_";
 266     namedistr+=i+1;
 267     TH2F *hcosthpointd0d0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 268     namedistr="hcosthpointd0d0B_";
 269     namedistr+=i+1;
 270     TH2F *hcosthpointd0d0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 271
 272     fDistr->Add(hptpiS);
 273     fDistr->Add(hptKS);
 274     fDistr->Add(hptB);
 275
 276     fDistr->Add(hptpiSnoMcut);
 277     fDistr->Add(hptKSnoMcut);
 278     fDistr->Add(hptB1pnoMcut);
 279     fDistr->Add(hptB2pnoMcut);
 280
 281     fDistr->Add(hdcaS);
 282     fDistr->Add(hdcaB);
 283
 284     fDistr->Add(hd0piS);
 285     fDistr->Add(hd0KS);
 286     fDistr->Add(hd0B);
 287
 288     fDistr->Add(hd0d0S);
 289     fDistr->Add(hd0d0B);
 290
 291     fDistr->Add(hcosthetastarS);
 292     fDistr->Add(hcosthetastarB);
 293
 294     fDistr->Add(hcosthetapointS);
 295     fDistr->Add(hcosthetapointB);
 296
 297     fDistr->Add(hcosthpointd0d0S);
 298     fDistr->Add(hcosthpointd0d0B);
 299
 300
 301     //histograms of invariant mass distributions
 302
 303     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 304     TH1F *tmpMl=(TH1F*)tmpMt->Clone();
 305     tmpMt->Sumw2();
 306     tmpMl->Sumw2();
 307
 308     //to compare with AliAnalysisTaskCharmFraction
 309     TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.765,1.965);
 310     TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
 311     tmpS27t->Sumw2();
 312     tmpS27l->Sumw2();
 313
 314     //distribution w/o cuts
 315     //   TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
 316     TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
 317     TH1F *tmpMB=(TH1F*)tmpMS->Clone();
 318     tmpMB->SetName(nameMassNocutsB.Data());
 319     tmpMS->Sumw2();
 320     tmpMB->Sumw2();
 321
 322     //MC signal and background
 323     TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 324     TH1F *tmpSl=(TH1F*)tmpSt->Clone();
 325     tmpSt->Sumw2();
 326     tmpSl->Sumw2();
 327
 328     TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 329     TH1F *tmpBl=(TH1F*)tmpBt->Clone();
 330     tmpBt->Sumw2();
 331     tmpBl->Sumw2();
 332
 333     //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 334     TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 335     TH1F *tmpRl=(TH1F*)tmpRt->Clone();
 336     tmpRt->Sumw2();
 337     tmpRl->Sumw2();
 338   //  printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
 339
 340     fOutputPPR->Add(tmpMt);
 341     fOutputPPR->Add(tmpSt);
 342     fOutputPPR->Add(tmpS27t);
 343     fOutputPPR->Add(tmpBt);
 344     fOutputPPR->Add(tmpRt);
 345
 346     fOutputmycuts->Add(tmpMl);
 347     fOutputmycuts->Add(tmpSl);
 348     fOutputmycuts->Add(tmpS27l);
 349     fOutputmycuts->Add(tmpBl);
 350     fOutputmycuts->Add(tmpRl);
 351
 352     fDistr->Add(tmpMS);
 353     fDistr->Add(tmpMB);
 354
 355   }
 356
 357
 358   fNentries=new TH1F("nentriesD0", "nentriesD0->Integral(1,2) = number of AODs *** nentriesD0->Integral(3,4) = number of candidates selected with cuts *** nentriesD0->Integral(5,6) = number of D0 selected with cuts", 6,1.,4.);
 359
 360   return;
 361 }
 362
 363 //________________________________________________________________________
 364 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 365 {
 366   // Execute analysis for current event:
 367   // heavy flavor candidates association to MC truth
 368   //cout<<"I'm in UserExec"<<endl;
 369   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 370
 371   TString bname;
 372   if(fArray==0){ //D0 candidates
 373     // load D0->Kpi candidates
 374     //cout<<"D0 candidates"<<endl;
 375     bname="D0toKpi";
 376   } else { //LikeSign candidates
 377     // load like sign candidates
 378     //cout<<"LS candidates"<<endl;
 379     bname="LikeSign2Prong";
 380   }
 381
 382   TClonesArray *inputArray=0;
 383
 384   if(!aod && AODEvent() && IsStandardAOD()) {
 385     // In case there is an AOD handler writing a standard AOD, use the AOD
 386     // event in memory rather than the input (ESD) event.
 387     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 388     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 389     // have to taken from the AOD event hold by the AliAODExtension
 390     AliAODHandler* aodHandler = (AliAODHandler*)
 391       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 392     if(aodHandler->GetExtensions()) {
 393       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 394       AliAODEvent *aodFromExt = ext->GetAOD();
 395       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 396     }
 397   } else {
 398     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 399   }
 400
 401
 402   if(!inputArray) {
 403     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
 404     return;
 405   }
 406
 407   // AOD primary vertex
 408   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 409   //vtx1->Print();
 410
 411   TClonesArray *mcArray = 0;
 412   AliAODMCHeader *mcHeader = 0;
 413
 414   if(fReadMC) {
 415     // load MC particles
 416     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 417     if(!mcArray) {
 418       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 419       return;
 420     }
 421
 422     // load MC header
 423     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 424     if(!mcHeader) {
 425       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 426       return;
 427     }
 428   }
 429
 430   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
 431
 432   //histogram filled with 1 for every AOD
 433   fNentries->Fill(1);
 434   PostData(3,fNentries);
 435
 436   // loop over candidates
 437   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 438   if(fDebug>1) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 439
 440   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 441     //Int_t nPosPairs=0, nNegPairs=0;
 442     //cout<<"inside the loop"<<endl;
 443     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 444     Bool_t unsetvtx=kFALSE;
 445     if(!d->GetOwnPrimaryVtx()) {
 446       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 447       unsetvtx=kTRUE;
 448     }
 449
 450     //check reco daughter in acceptance
 451     Double_t eta0=d->EtaProng(0);
 452     Double_t eta1=d->EtaProng(1);
 453     Bool_t prongsinacc=kFALSE;
 454     if (TMath::Abs(eta0) < 0.9 && TMath::Abs(eta1) < 0.9) prongsinacc=kTRUE;
 455
 456     //add distr here
 457     UInt_t pdgs[2];
 458
 459     Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 460     pdgs[0]=211;
 461     pdgs[1]=321;
 462     Double_t minvD0 = d->InvMassD0();
 463     pdgs[1]=211;
 464     pdgs[0]=321;
 465     Double_t minvD0bar = d->InvMassD0bar();
 466     //cout<<"inside mass cut"<<endl;
 467
 468     Int_t pdgDgD0toKpi[2]={321,211};
 469     Int_t lab=-9999;
 470     if(fReadMC) lab=d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 471     Double_t pt = d->Pt(); //mother pt
 472
 473     if(lab>=0 && fReadMC){ //signal
 474
 475       //check pdg of the prongs
 476       AliAODTrack *prong0=(AliAODTrack*)d->GetDaughter(0);
 477       AliAODTrack *prong1=(AliAODTrack*)d->GetDaughter(1);
 478       Int_t labprong[2];
 479       labprong[0]=prong0->GetLabel();
 480       labprong[1]=prong1->GetLabel();
 481       AliAODMCParticle *mcprong=0;
 482       Int_t PDGprong[2]={0,0};
 483       for (Int_t iprong=0;iprong<2;iprong++){
 484         if(labprong[iprong]>=0)  mcprong= (AliAODMCParticle*)mcArray->At(labprong[iprong]);
 485         PDGprong[iprong]=mcprong->GetPdgCode();
 486       }
 487
 488       //no mass cut ditributions: ptbis
 489
 490       if(pt>0. && pt<=1.) {
 491         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
 492           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_1"))->Fill(d->PtProng(0));
 493           ((TH1F*)fDistr->FindObject("hptKSnoMcut_1"))->Fill(d->PtProng(1));
 494         }else {
 495           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
 496             ((TH1F*)fDistr->FindObject("hptKSnoMcut_1"))->Fill(d->PtProng(0));
 497             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_1"))->Fill(d->PtProng(1));
 498           }
 499         }
 500       }
 501       if(pt>1. && pt<=2.) {
 502         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
 503           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_2"))->Fill(d->PtProng(0));
 504           ((TH1F*)fDistr->FindObject("hptKSnoMcut_2"))->Fill(d->PtProng(1));
 505         }else {
 506           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
 507             ((TH1F*)fDistr->FindObject("hptKSnoMcut_2"))->Fill(d->PtProng(0));
 508             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_2"))->Fill(d->PtProng(1));
 509           }
 510         }
 511       }
 512       if(pt>2. && pt<=3.) {
 513         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
 514           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_3"))->Fill(d->PtProng(0));
 515           ((TH1F*)fDistr->FindObject("hptKSnoMcut_3"))->Fill(d->PtProng(1));
 516         }else {
 517           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
 518             ((TH1F*)fDistr->FindObject("hptKSnoMcut_3"))->Fill(d->PtProng(0));
 519             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_3"))->Fill(d->PtProng(1));
 520           }
 521         }
 522       }
 523       if(pt>3. && pt<=5.) {
 524
 525         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
 526           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_4"))->Fill(d->PtProng(0));
 527           ((TH1F*)fDistr->FindObject("hptKSnoMcut_4"))->Fill(d->PtProng(1));
 528         }else {
 529           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
 530             ((TH1F*)fDistr->FindObject("hptKSnoMcut_4"))->Fill(d->PtProng(0));
 531             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_4"))->Fill(d->PtProng(1));
 532           }
 533         }
 534       }
 535       if(pt>5.)           {
 536
 537         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
 538           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_5"))->Fill(d->PtProng(0));
 539           ((TH1F*)fDistr->FindObject("hptKSnoMcut_5"))->Fill(d->PtProng(1));
 540         }else {
 541           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
 542             ((TH1F*)fDistr->FindObject("hptKSnoMcut_5"))->Fill(d->PtProng(0));
 543             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_5"))->Fill(d->PtProng(1));
 544           }
 545         }
 546       }
 547
 548       if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421){//D0
 549
 550         //no mass cut ditributions: mass, costhetastar
 551
 552         if(pt>0. && pt<=1.) {
 553           ((TH1F*)fDistr->FindObject("hMassS_1"))->Fill(minvD0);
 554         }
 555         if(pt>1. && pt<=2.) {
 556           ((TH1F*)fDistr->FindObject("hMassS_2"))->Fill(minvD0);
 557
 558         }
 559         if(pt>2. && pt<=3.) {
 560           ((TH1F*)fDistr->FindObject("hMassS_3"))->Fill(minvD0);
 561
 562         }
 563         if(pt>3. && pt<=5.) {
 564           ((TH1F*)fDistr->FindObject("hMassS_4"))->Fill(minvD0);
 565
 566         }
 567         if(pt>5.)           {
 568           ((TH1F*)fDistr->FindObject("hMassS_5"))->Fill(minvD0);
 569
 570         }
 571
 572       }
 573       else { //D0bar
 574
 575         //no mass cut ditributions: mass
 576
 577         if(pt>0. && pt<=1.) {
 578           ((TH1F*)fDistr->FindObject("hMassS_1"))->Fill(minvD0bar);
 579
 580         }
 581         if(pt>1. && pt<=2.) {
 582           ((TH1F*)fDistr->FindObject("hMassS_2"))->Fill(minvD0bar);
 583
 584         }
 585         if(pt>2. && pt<=3.) {
 586           ((TH1F*)fDistr->FindObject("hMassS_3"))->Fill(minvD0bar);
 587
 588         }
 589         if(pt>3. && pt<=5.) {
 590           ((TH1F*)fDistr->FindObject("hMassS_4"))->Fill(minvD0bar);
 591
 592         }
 593         if(pt>5.)           {
 594           ((TH1F*)fDistr->FindObject("hMassS_5"))->Fill(minvD0bar);
 595         }
 596
 597       }
 598
 599       //apply cut on invariant mass on the pair
 600       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
 601
 602         if(fArray==1) cout<<"LS signal: ERROR"<<endl;
 603         for (Int_t iprong=0; iprong<2; iprong++){
 604           AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(iprong);
 605           Int_t labprong=prong->GetLabel();
 606
 607           //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
 608           AliAODMCParticle *mcprong=0;
 609           if(labprong>=0)  mcprong= (AliAODMCParticle*)mcArray->At(labprong);
 610           Int_t pdgprong=mcprong->GetPdgCode();
 611           if(TMath::Abs(pdgprong)==211) {
 612             //cout<<"pi"<<endl;
 613             if(pt>0. && pt<=1.){
 614               ((TH1F*)fDistr->FindObject("hptpiS_1"))->Fill(d->PtProng(iprong));
 615               ((TH1F*)fDistr->FindObject("hd0piS_1"))->Fill(d->Getd0Prong(iprong));
 616             }
 617             if(pt>1. && pt<=2.){
 618               ((TH1F*)fDistr->FindObject("hptpiS_2"))->Fill(d->PtProng(iprong));
 619               ((TH1F*)fDistr->FindObject("hd0piS_2"))->Fill(d->Getd0Prong(iprong));
 620             }
 621
 622             if(pt>2. && pt<=3.){
 623               ((TH1F*)fDistr->FindObject("hptpiS_3"))->Fill(d->PtProng(iprong));
 624               ((TH1F*)fDistr->FindObject("hd0piS_3"))->Fill(d->Getd0Prong(iprong));
 625
 626             }
 627             if(pt>3. && pt<=5.){
 628               ((TH1F*)fDistr->FindObject("hptpiS_4"))->Fill(d->PtProng(iprong));
 629               ((TH1F*)fDistr->FindObject("hd0piS_4"))->Fill(d->Getd0Prong(iprong));
 630
 631             }
 632             if(pt>5.)          {
 633               ((TH1F*)fDistr->FindObject("hptpiS_5"))->Fill(d->PtProng(iprong));
 634               ((TH1F*)fDistr->FindObject("hd0piS_5"))->Fill(d->Getd0Prong(iprong));
 635
 636             }
 637
 638
 639         }
 640
 641         if(TMath::Abs(pdgprong)==321) {
 642             //cout<<"kappa"<<endl;
 643             if(pt>0. && pt<=1.){
 644               ((TH1F*)fDistr->FindObject("hptKS_1"))->Fill(d->PtProng(iprong));
 645               ((TH1F*)fDistr->FindObject("hd0KS_1"))->Fill(d->Getd0Prong(iprong));
 646             }
 647             if(pt>1. && pt<=2.){
 648               ((TH1F*)fDistr->FindObject("hptKS_2"))->Fill(d->PtProng(iprong));
 649               ((TH1F*)fDistr->FindObject("hd0KS_2"))->Fill(d->Getd0Prong(iprong));
 650             }
 651             if(pt>2. && pt<=3.){
 652               ((TH1F*)fDistr->FindObject("hptKS_3"))->Fill(d->PtProng(iprong));
 653               ((TH1F*)fDistr->FindObject("hd0KS_3"))->Fill(d->Getd0Prong(iprong));
 654             }
 655             if(pt>3. && pt<=5.){
 656               ((TH1F*)fDistr->FindObject("hptKS_4"))->Fill(d->PtProng(iprong));
 657               ((TH1F*)fDistr->FindObject("hd0KS_4"))->Fill(d->Getd0Prong(iprong));
 658
 659             }
 660             if(pt>5.)          {
 661               ((TH1F*)fDistr->FindObject("hptKS_5"))->Fill(d->PtProng(iprong));
 662               ((TH1F*)fDistr->FindObject("hd0KS_5"))->Fill(d->Getd0Prong(iprong));
 663
 664             }
 665
 666           }
 667
 668           if(pt>0. && pt<=1.){
 669             ((TH1F*)fDistr->FindObject("hdcaS_1"))->Fill(d->GetDCA());
 670
 671             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 672               ((TH1F*)fDistr->FindObject("hcosthetastarS_1"))->Fill(d->CosThetaStarD0());
 673             else ((TH1F*)fDistr->FindObject("hcosthetastarS_1"))->Fill(d->CosThetaStarD0bar());
 674
 675           }
 676           if(pt>1. && pt<=2.){
 677             ((TH1F*)fDistr->FindObject("hdcaS_2"))->Fill(d->GetDCA());
 678
 679             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 680               ((TH1F*)fDistr->FindObject("hcosthetastarS_2"))->Fill(d->CosThetaStarD0());
 681             else ((TH1F*)fDistr->FindObject("hcosthetastarS_2"))->Fill(d->CosThetaStarD0bar());
 682
 683           }
 684           if(pt>2. && pt<=3.){
 685             ((TH1F*)fDistr->FindObject("hdcaS_3"))->Fill(d->GetDCA());
 686
 687             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 688               ((TH1F*)fDistr->FindObject("hcosthetastarS_3"))->Fill(d->CosThetaStarD0());
 689             else ((TH1F*)fDistr->FindObject("hcosthetastarS_3"))->Fill(d->CosThetaStarD0bar());
 690
 691           }
 692           if(pt>3. && pt<=5.){
 693             ((TH1F*)fDistr->FindObject("hdcaS_4"))->Fill(d->GetDCA());
 694
 695             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 696               ((TH1F*)fDistr->FindObject("hcosthetastarS_4"))->Fill(d->CosThetaStarD0());
 697             else ((TH1F*)fDistr->FindObject("hcosthetastarS_4"))->Fill(d->CosThetaStarD0bar());
 698
 699           }
 700           if(pt>5.)          {
 701             ((TH1F*)fDistr->FindObject("hdcaS_5"))->Fill(d->GetDCA());
 702
 703             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 704               ((TH1F*)fDistr->FindObject("hcosthetastarS_5"))->Fill(d->CosThetaStarD0());
 705             else ((TH1F*)fDistr->FindObject("hcosthetastarS_5"))->Fill(d->CosThetaStarD0bar());
 706
 707           }
 708
 709         }
 710
 711         if(pt>0. && pt<=1.){
 712           ((TH1F*)fDistr->FindObject("hd0d0S_1"))->Fill(d->Prodd0d0());
 713           ((TH1F*)fDistr->FindObject("hcosthetapointS_1"))->Fill(d->CosPointingAngle());
 714           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_1"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 715         }
 716         if(pt>1. && pt<=2.){
 717           ((TH1F*)fDistr->FindObject("hd0d0S_2"))->Fill(d->Prodd0d0());
 718           ((TH1F*)fDistr->FindObject("hcosthetapointS_2"))->Fill(d->CosPointingAngle());
 719           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_2"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 720         }
 721         if(pt>2. && pt<=3.){
 722           ((TH1F*)fDistr->FindObject("hd0d0S_3"))->Fill(d->Prodd0d0());
 723           ((TH1F*)fDistr->FindObject("hcosthetapointS_3"))->Fill(d->CosPointingAngle());
 724           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_3"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 725         }
 726         if(pt>3. && pt<=5.){
 727           ((TH1F*)fDistr->FindObject("hd0d0S_4"))->Fill(d->Prodd0d0());
 728           ((TH1F*)fDistr->FindObject("hcosthetapointS_4"))->Fill(d->CosPointingAngle());
 729           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_4"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 730         }
 731         if(pt>5.)          {
 732           ((TH1F*)fDistr->FindObject("hd0d0S_5"))->Fill(d->Prodd0d0());
 733           ((TH1F*)fDistr->FindObject("hcosthetapointS_5"))->Fill(d->CosPointingAngle());
 734           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_5"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 735         }
 736
 737       } //invmass cut
 738
 739     } else{ //Background or LS
 740       //cout<<"is background"<<endl;
 741       Double_t pt = d->Pt();
 742
 743       //no mass cut distributions: mass, ptbis
 744       if(pt>0. && pt<=1.) {
 745         ((TH1F*)fDistr->FindObject("hMassB_1"))->Fill(minvD0);
 746         ((TH1F*)fDistr->FindObject("hMassB_1"))->Fill(minvD0bar);
 747
 748         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_1"))->Fill(d->PtProng(0));
 749         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_1"))->Fill(d->PtProng(0));
 750         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_1"))->Fill(d->PtProng(1));
 751
 752       }
 753       if(pt>1. && pt<=2.) {
 754         ((TH1F*)fDistr->FindObject("hMassB_2"))->Fill(minvD0);
 755         ((TH1F*)fDistr->FindObject("hMassB_2"))->Fill(minvD0bar);
 756
 757         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_2"))->Fill(d->PtProng(0));
 758         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_2"))->Fill(d->PtProng(0));
 759         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_2"))->Fill(d->PtProng(1));
 760
 761       }
 762       if(pt>2. && pt<=3.) {
 763         ((TH1F*)fDistr->FindObject("hMassB_3"))->Fill(minvD0);
 764         ((TH1F*)fDistr->FindObject("hMassB_3"))->Fill(minvD0bar);
 765
 766         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_3"))->Fill(d->PtProng(0));
 767         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_3"))->Fill(d->PtProng(0));
 768         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_3"))->Fill(d->PtProng(1));
 769
 770       }
 771       if(pt>3. && pt<=5.) {
 772         ((TH1F*)fDistr->FindObject("hMassB_4"))->Fill(minvD0);
 773         ((TH1F*)fDistr->FindObject("hMassB_4"))->Fill(minvD0bar);
 774
 775         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_4"))->Fill(d->PtProng(0));
 776         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_4"))->Fill(d->PtProng(0));
 777         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_4"))->Fill(d->PtProng(1));
 778
 779       }
 780       if(pt>5.)           {
 781         ((TH1F*)fDistr->FindObject("hMassB_5"))->Fill(minvD0);
 782         ((TH1F*)fDistr->FindObject("hMassB_5"))->Fill(minvD0bar);
 783
 784         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_5"))->Fill(d->PtProng(0));
 785         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_5"))->Fill(d->PtProng(0));
 786         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_5"))->Fill(d->PtProng(1));
 787
 788       }
 789
 790       //apply cut on invariant mass on the pair
 791       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
 792
 793
 794         AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(0);
 795         if(!prong) cout<<"No daughter found";
 796         else{
 797           if(prong->Charge()==1) {fTotPosPairs[4]++;} else {fTotNegPairs[4]++;}
 798         }
 799         if(pt>0. && pt<=1.){
 800           //normalise pt distr to half afterwards
 801           ((TH1F*)fDistr->FindObject("hptB_1"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_1"))->Fill(d->PtProng(1));
 802           ((TH1F*)fDistr->FindObject("hd0B_1"))->Fill(d->Getd0Prong(0));
 803           ((TH1F*)fDistr->FindObject("hdcaB_1"))->Fill(d->GetDCA());
 804           ((TH1F*)fDistr->FindObject("hcosthetastarB_1"))->Fill(d->CosThetaStarD0());
 805           ((TH1F*)fDistr->FindObject("hcosthetastarB_1"))->Fill(d->CosThetaStarD0bar());
 806           ((TH1F*)fDistr->FindObject("hd0d0B_1"))->Fill(d->Prodd0d0());
 807           ((TH1F*)fDistr->FindObject("hcosthetapointB_1"))->Fill(d->CosPointingAngle());
 808           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_1"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 809         }
 810         if(pt>1. && pt<=2.){
 811           ((TH1F*)fDistr->FindObject("hptB_2"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_2"))->Fill(d->PtProng(1));
 812           ((TH1F*)fDistr->FindObject("hd0B_2"))->Fill(d->Getd0Prong(0));
 813           ((TH1F*)fDistr->FindObject("hdcaB_2"))->Fill(d->GetDCA());
 814           ((TH1F*)fDistr->FindObject("hcosthetastarB_2"))->Fill(d->CosThetaStarD0());
 815           ((TH1F*)fDistr->FindObject("hcosthetastarB_2"))->Fill(d->CosThetaStarD0bar());
 816           ((TH1F*)fDistr->FindObject("hd0d0B_2"))->Fill(d->Prodd0d0());
 817           ((TH1F*)fDistr->FindObject("hcosthetapointB_2"))->Fill(d->CosPointingAngle());
 818           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_2"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 819         }
 820         if(pt>2. && pt<=3.){
 821           ((TH1F*)fDistr->FindObject("hptB_3"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_3"))->Fill(d->PtProng(1));
 822           ((TH1F*)fDistr->FindObject("hd0B_3"))->Fill(d->Getd0Prong(0));
 823           ((TH1F*)fDistr->FindObject("hdcaB_3"))->Fill(d->GetDCA());
 824           ((TH1F*)fDistr->FindObject("hcosthetastarB_3"))->Fill(d->CosThetaStarD0());
 825           ((TH1F*)fDistr->FindObject("hcosthetastarB_3"))->Fill(d->CosThetaStarD0bar());
 826           ((TH1F*)fDistr->FindObject("hd0d0B_3"))->Fill(d->Prodd0d0());
 827           ((TH1F*)fDistr->FindObject("hcosthetapointB_3"))->Fill(d->CosPointingAngle());
 828           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_3"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 829         }
 830         if(pt>3. && pt<=5.){
 831           ((TH1F*)fDistr->FindObject("hptB_4"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_4"))->Fill(d->PtProng(1));
 832           ((TH1F*)fDistr->FindObject("hd0B_4"))->Fill(d->Getd0Prong(0));
 833           ((TH1F*)fDistr->FindObject("hdcaB_4"))->Fill(d->GetDCA());
 834           ((TH1F*)fDistr->FindObject("hcosthetastarB_4"))->Fill(d->CosThetaStarD0());
 835           ((TH1F*)fDistr->FindObject("hcosthetastarB_4"))->Fill(d->CosThetaStarD0bar());
 836           ((TH1F*)fDistr->FindObject("hd0d0B_4"))->Fill(d->Prodd0d0());
 837           ((TH1F*)fDistr->FindObject("hcosthetapointB_4"))->Fill(d->CosPointingAngle());
 838           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_4"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 839         }
 840         if(pt>5.)         {
 841           ((TH1F*)fDistr->FindObject("hptB_5"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_5"))->Fill(d->PtProng(1));
 842           ((TH1F*)fDistr->FindObject("hd0B_5"))->Fill(d->Getd0Prong(0));
 843           ((TH1F*)fDistr->FindObject("hdcaB_5"))->Fill(d->GetDCA());
 844           ((TH1F*)fDistr->FindObject("hcosthetastarB_5"))->Fill(d->CosThetaStarD0());
 845           ((TH1F*)fDistr->FindObject("hcosthetastarB_5"))->Fill(d->CosThetaStarD0bar());
 846           ((TH1F*)fDistr->FindObject("hd0d0B_5"))->Fill(d->Prodd0d0());
 847           ((TH1F*)fDistr->FindObject("hcosthetapointB_5"))->Fill(d->CosPointingAngle());
 848           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_5"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 849         }
 850
 851
 852       }// end if inv mass cut
 853     }//end if background
 854
 855
 856
 857      //cuts order
 858 //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 859 //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 860 //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 861 //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 862 //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 863 //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 864 //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 865 //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 866 //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 867
 868     Int_t ptbin=0;
 869     Bool_t isInRange=kFALSE;
 870
 871     //cout<<"P_t = "<<pt<<endl;
 872     if (pt>0. && pt<=1.) {
 873       ptbin=1;
 874       isInRange=kTRUE;
 875       /*
 876       //test d0 cut
 877       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.1,-0.0002,0.5);
 878       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.1,-0.00025,0.7);
 879       */
 880       //original
 881       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.5);
 882       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.00025,0.7);
 883
 884 //       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.7);
 885 //       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,1,1,-0.00015,0.5);
 886       //printf("I'm in the bin %d\n",ptbin);
 887
 888     }
 889
 890     if(pt>1. && pt<=3.) {
 891       if(pt>1. && pt<=2.) ptbin=2;
 892       if(pt>2. && pt<=3.) ptbin=3;
 893       isInRange=kTRUE;
 894       /*
 895       //test d0 cut
 896       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.0002,0.6);
 897       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,0.1,-0.00025,0.8);
 898       */
 899       //original
 900       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0002,0.6);
 901       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
 902
 903       //printf("I'm in the bin %d\n",ptbin);
 904     }
 905
 906     if(pt>3. && pt<=5.){
 907       ptbin=4;
 908       isInRange=kTRUE;
 909       /*
 910       //test d0 cut
 911       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.0001,0.8);
 912       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00015,0.8);
 913       */
 914       //original
 915       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0001,0.8);
 916       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.8);
 917
 918       //printf("I'm in the bin %d\n",ptbin);
 919     }
 920     if(pt>5.&& pt<=10.){ //additional upper limit to compare with Correction Framework
 921       ptbin=5;
 922       isInRange=kTRUE;
 923       /*
 924       //test d0 cut
 925       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00005,0.8);
 926       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00015,0.9);
 927       */
 928       //original
 929       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00005,0.8);
 930       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.9);
 931
 932     }//if(pt>5)    if (pt>0. && pt<=1.) {
 933     //printf("I'm in the bin %d\n",ptbin);
 934     //old
 935     //fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed
 936     if(prongsinacc && isInRange){
 937       FillHists(ptbin,d,mcArray,fVHFPPR,fOutputPPR);
 938       FillHists(ptbin,d,mcArray,fVHFmycuts,fOutputmycuts);
 939     }
 940     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 941   } //end for prongs
 942
 943
 944
 945
 946   // Post the data
 947   PostData(1,fOutputPPR);
 948   PostData(2,fOutputmycuts);
 949   PostData(4,fDistr);
 950   return;
 951 }
 952 //____________________________________________________________________________*
 953 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAnalysisVertexingHF *vhf, TList *listout){
 954   //
 955   // function used in UserExec:
 956   //
 957
 958
 959   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 960
 961   Int_t okD0=0,okD0bar=0;
 962   //cout<<"inside FillHist"<<endl;
 963   if(part->SelectD0(vhf->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
 964     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 965     //printf("SELECTED\n");
 966
 967
 968     AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
 969     if(!prong) cout<<"No daughter found";
 970     else{
 971       if(prong->Charge()==1) {fTotPosPairs[ptbin-1]++;} else {fTotNegPairs[ptbin-1]++;}
 972     }
 973
 974     TString fillthis="";
 975     Int_t pdgDgD0toKpi[2]={321,211};
 976     Int_t labD0=-1;
 977     if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 978
 979     //count candidates selected by cuts
 980     fNentries->Fill(2);
 981     //count true D0 selected by cuts
 982     if (fReadMC && labD0>=0) fNentries->Fill(3);
 983     PostData(3,fNentries);
 984
 985     if (okD0==1) {
 986       fillthis="histMass_";
 987       fillthis+=ptbin;
 988       //cout<<"Filling "<<fillthis<<endl;
 989
 990       //printf("Fill mass with D0");
 991       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 992
 993       if(labD0>=0) {
 994         if(fArray==1) cout<<"LS signal ERROR"<<endl;
 995
 996         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 997         Int_t pdgD0 = partD0->GetPdgCode();
 998         //cout<<"pdg = "<<pdgD0<<endl;
 999         if (pdgD0==421){ //D0
1000           //cout<<"Fill S with D0"<<endl;
1001           fillthis="histSgn_";
1002           fillthis+=ptbin;
1003           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1004           if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1005             fillthis="histSgn27_";
1006             fillthis+=ptbin;
1007             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1008           }
1009         } else{ //it was a D0bar
1010           fillthis="histRfl_";
1011           fillthis+=ptbin;
1012           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1013         }
1014       } else {//background
1015         fillthis="histBkg_";
1016         fillthis+=ptbin;
1017         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1018       }
1019
1020     }
1021     if (okD0bar==1) {
1022       fillthis="histMass_";
1023       fillthis+=ptbin;
1024       //printf("Fill mass with D0bar");
1025       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1026
1027       if(labD0>=0) {
1028         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1029         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1030         Int_t pdgD0 = partD0->GetPdgCode();
1031         //cout<<" pdg = "<<pdgD0<<endl;
1032         if (pdgD0==-421){ //D0bar
1033           fillthis="histSgn_";
1034           fillthis+=ptbin;
1035           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1036           if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1037             fillthis="histSgn27_";
1038             fillthis+=ptbin;
1039             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1040           }
1041
1042
1043         } else{
1044           fillthis="histRfl_";
1045           fillthis+=ptbin;
1046           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1047         }
1048       } else {//background or LS
1049         fillthis="histBkg_";
1050         fillthis+=ptbin;
1051         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1052       }
1053     }
1054
1055
1056   } //else cout<<"NOT SELECTED"<<endl;
1057
1058 }
1059 //________________________________________________________________________
1060 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1061 {
1062   // Terminate analysis
1063   //
1064   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1065
1066
1067   fOutputPPR = dynamic_cast<TList*> (GetOutputData(1));
1068   if (!fOutputPPR) {
1069     printf("ERROR: fOutputthight not available\n");
1070     return;
1071   }
1072   fOutputmycuts = dynamic_cast<TList*> (GetOutputData(2));
1073   if (!fOutputmycuts) {
1074     printf("ERROR: fOutputmycuts not available\n");
1075     return;
1076   }
1077   fDistr = dynamic_cast<TList*> (GetOutputData(4));
1078   if (!fDistr) {
1079     printf("ERROR: fDistr not available\n");
1080     return;
1081   }
1082
1083   if(fArray==1){
1084     for(Int_t ipt=0;ipt<4;ipt++){
1085       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1086
1087
1088       if(fLsNormalization>0) {
1089
1090         TString massName="histMass_";
1091         massName+=ipt+1;
1092         ((TH1F*)fOutputPPR->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputPPR->FindObject(massName))->GetEntries());
1093         ((TH1F*)fOutputmycuts->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputmycuts->FindObject(massName))->GetEntries());
1094       }
1095     }
1096
1097     fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1098
1099     if(fLsNormalization>0) {
1100
1101       TString nameDistr="hptB";
1102       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1103       nameDistr="hdcaB";
1104       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1105       nameDistr="hcosthetastarB";
1106       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1107       nameDistr="hd0B";
1108       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1109       nameDistr="hd0d0B";
1110       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1111       nameDistr="hcosthetapointB";
1112       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1113       nameDistr="hcosthpointd0d0B";
1114       ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1115
1116     }
1117   }
1118
1119   TString cvname;
1120
1121   if (fArray==0){
1122     cvname="D0invmass";
1123   } else cvname="LSinvmass";
1124
1125   TCanvas *c1=new TCanvas(cvname,cvname);
1126   c1->cd();
1127   ((TH1F*)fOutputPPR->FindObject("histMass_4"))->Draw();
1128
1129
1130   return;
1131 }
1132