PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx

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