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