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 <TNtuple.h>
  29 #include <TList.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TDatabasePDG.h>
  33
  34 #include "AliAODEvent.h"
  35 #include "AliAODVertex.h"
  36 #include "AliAODTrack.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODRecoDecayHF2Prong.h"
  40 #include "AliAODRecoCascadeHF.h"
  41 #include "AliAnalysisVertexingHF.h"
  42 #include "AliAnalysisTaskSE.h"
  43 #include "AliAnalysisTaskSED0Mass.h"
  44
  45 ClassImp(AliAnalysisTaskSED0Mass)
  46
  47
  48 //________________________________________________________________________
  49 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  50 AliAnalysisTaskSE(),
  51 fOutputPPR(0),
  52 fOutputloose(0),
  53 fDistr(0),
  54 fNentries(0),
  55 fVHFPPR(0),
  56 fVHFloose(0),
  57 fArray(0),
  58 fLsNormalization(1.)
  59
  60 {
  61   // Default constructor
  62   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
  67 AliAnalysisTaskSE(name),
  68 fOutputPPR(0),
  69 fOutputloose(0),
  70 fDistr(0),
  71 fNentries(0),
  72 fVHFPPR(0),
  73 fVHFloose(0),
  74 fArray(0),
  75 fLsNormalization(1.)
  76 {
  77   // Default constructor
  78   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  79
  80   // Output slot #1 writes into a TList container
  81   DefineOutput(1,TList::Class());  //My private output
  82   // Output slot #2 writes into a TList container
  83   DefineOutput(2,TList::Class());  //My private output
  84   // Output slot #3 writes into a TH1F container
  85   DefineOutput(3,TH1F::Class());  //My private output
  86   // Output slot #4 writes into a TList container
  87   DefineOutput(4,TList::Class());  //My private output
  88 }
  89
  90 //________________________________________________________________________
  91 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
  92 {
  93    if (fOutputPPR) {
  94     delete fOutputPPR;
  95     fOutputPPR = 0;
  96   }
  97   if (fVHFPPR) {
  98     delete fVHFPPR;
  99     fVHFPPR = 0;
 100   }
 101   if (fOutputloose) {
 102     delete fOutputloose;
 103     fOutputloose = 0;
 104   }
 105
 106   if (fDistr) {
 107     delete fDistr;
 108     fDistr = 0;
 109   }
 110
 111   if (fVHFloose) {
 112     delete fVHFloose;
 113     fVHFloose = 0;
 114   }
 115   if (fNentries){
 116     delete fNentries;
 117     fNentries = 0;
 118   }
 119
 120
 121 }
 122
 123 //________________________________________________________________________
 124 void AliAnalysisTaskSED0Mass::Init()
 125 {
 126   // Initialization
 127
 128   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 129
 130   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
 131
 132   // 2 sets of dedidcated cuts -- defined in UserExec
 133   fVHFPPR = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 134   fVHFloose = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 135
 136   return;
 137 }
 138
 139 //________________________________________________________________________
 140 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 141 {
 142
 143   // Create the output container
 144   //
 145   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 146
 147   // Several histograms are more conveniently managed in a TList
 148   fOutputPPR = new TList();
 149   fOutputPPR->SetOwner();
 150   fOutputPPR->SetName("listPPR");
 151
 152   fOutputloose = new TList();
 153   fOutputloose->SetOwner();
 154   fOutputloose->SetName("listloose");
 155
 156   fDistr = new TList();
 157   fDistr->SetOwner();
 158   fDistr->SetName("distributionslist");
 159
 160   const Int_t nhist=5;
 161
 162   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ", namedistr=" ";
 163
 164   for(Int_t i=0;i<nhist;i++){
 165     nameMass="histMass_";
 166     nameMass+=i+1;
 167     nameSgn27="histSgn27_";
 168     nameSgn27+=i+1;
 169     nameSgn="histSgn_";
 170     nameSgn+=i+1;
 171     nameBkg="histBkg_";
 172     nameBkg+=i+1;
 173     nameRfl="histRfl_";
 174     nameRfl+=i+1;
 175
 176     //histograms of invariant mass distributions
 177
 178     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
 179     TH1F *tmpMl=(TH1F*)tmpMt->Clone();
 180     tmpMt->Sumw2();
 181     tmpMl->Sumw2();
 182
 183     //to compare with AliAnalysisTaskCharmFraction
 184     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);
 185     TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
 186     tmpS27t->Sumw2();
 187     tmpS27l->Sumw2();
 188
 189     TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 190     TH1F *tmpSl=(TH1F*)tmpSt->Clone();
 191     tmpSt->Sumw2();
 192     tmpSl->Sumw2();
 193
 194     TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 195     TH1F *tmpBl=(TH1F*)tmpBt->Clone();
 196     tmpBt->Sumw2();
 197     tmpBl->Sumw2();
 198
 199     //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 200     TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
 201     TH1F *tmpRl=(TH1F*)tmpRt->Clone();
 202     tmpRt->Sumw2();
 203     tmpRl->Sumw2();
 204   //  printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
 205
 206     fOutputPPR->Add(tmpMt);
 207     fOutputPPR->Add(tmpSt);
 208     fOutputPPR->Add(tmpS27t);
 209     fOutputPPR->Add(tmpBt);
 210     fOutputPPR->Add(tmpRt);
 211
 212     fOutputloose->Add(tmpMl);
 213     fOutputloose->Add(tmpSl);
 214     fOutputloose->Add(tmpS27l);
 215     fOutputloose->Add(tmpBl);
 216     fOutputloose->Add(tmpRl);
 217
 218   }
 219
 220   //histograms of cut variable distributions
 221   //  pT
 222   namedistr="hptpiS";
 223   TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 224
 225   namedistr="hptKS";
 226   TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 227
 228   namedistr="hptB";
 229   TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 230
 231   //  dca
 232   namedistr="hdcaS";
 233   TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 234   namedistr="hdcaB";
 235   TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 236
 237   //  costhetastar
 238   namedistr="hcosthetastarS";
 239   TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 240   namedistr="hcosthetastarB";
 241   TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 242
 243   // impact parameter
 244   namedistr="hd0piS";
 245   TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
 246
 247   namedistr="hd0KS";
 248   TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
 249   namedistr="hd0B";
 250   TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution;d0 [cm]",200,-0.1,0.1);
 251
 252   namedistr="hd0d0S";
 253   TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 254   namedistr="hd0d0B";
 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   TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 260   namedistr="hcosthetapointB";
 261   TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 262
 263   namedistr="hcosthpointd0d0S";
 264   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);
 265   namedistr="hcosthpointd0d0B";
 266   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);
 267
 268   fDistr->Add(hptpiS);
 269   fDistr->Add(hptKS);
 270   fDistr->Add(hptB);
 271
 272   fDistr->Add(hdcaS);
 273   fDistr->Add(hdcaB);
 274
 275   fDistr->Add(hd0piS);
 276   fDistr->Add(hd0KS);
 277   fDistr->Add(hd0B);
 278
 279   fDistr->Add(hd0d0S);
 280   fDistr->Add(hd0d0B);
 281
 282   fDistr->Add(hcosthetastarS);
 283   fDistr->Add(hcosthetastarB);
 284
 285   fDistr->Add(hcosthetapointS);
 286   fDistr->Add(hcosthetapointB);
 287
 288   fDistr->Add(hcosthpointd0d0S);
 289   fDistr->Add(hcosthpointd0d0B);
 290
 291   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.);
 292
 293   return;
 294 }
 295
 296 //________________________________________________________________________
 297 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 298 {
 299   // Execute analysis for current event:
 300   // heavy flavor candidates association to MC truth
 301   //cout<<"I'm in UserExec"<<endl;
 302   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 303
 304   TClonesArray *inputArray=0;
 305
 306   if(fArray==0){ //D0 candidates
 307     // load D0->Kpi candidates
 308     //cout<<"D0 candidates"<<endl;
 309     inputArray =
 310       (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 311     if(!inputArray) {
 312       printf("AliAnalysisTaskSECompareHFpt::UserExec: D0toKpi branch not found!\n");
 313       return;
 314     }
 315
 316   } else { //LikeSign candidates
 317
 318     //cout<<"LS candidates"<<endl;
 319     // load like sign candidates
 320     inputArray =
 321       (TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
 322     if(!inputArray) {
 323       printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: LikeSign2Prong branch not found!\n");
 324       return;
 325     }
 326
 327   }
 328
 329   // AOD primary vertex
 330   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 331   //vtx1->Print();
 332
 333   // load MC particles
 334   TClonesArray *mcArray =
 335     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 336   if(!mcArray) {
 337     printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 338     return;
 339   }
 340
 341   // load MC header
 342   AliAODMCHeader *mcHeader =
 343     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 344   if(!mcHeader) {
 345     printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 346     return;
 347   }
 348
 349
 350   //histogram filled with 1 for every AOD
 351   fNentries->Fill(1);
 352   PostData(3,fNentries);
 353
 354   // loop over candidates
 355   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 356   if(fDebug>1) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 357
 358   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 359     //Int_t nPosPairs=0, nNegPairs=0;
 360     //cout<<"inside the loop"<<endl;
 361     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 362     Bool_t unsetvtx=kFALSE;
 363     if(!d->GetOwnPrimaryVtx()) {
 364       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 365       unsetvtx=kTRUE;
 366     }
 367
 368
 369
 370     //check reco daughter in acceptance
 371     Double_t eta0=d->EtaProng(0);
 372     Double_t eta1=d->EtaProng(1);
 373     Bool_t prongsinacc=kFALSE;
 374     if (TMath::Abs(eta0) < 0.9 && TMath::Abs(eta1) < 0.9) prongsinacc=kTRUE;
 375
 376     //add distr here
 377     UInt_t pdgs[2];
 378
 379     Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 380     pdgs[0]=211;
 381     pdgs[1]=321;
 382     Double_t minvD0 = d->InvMassD0();
 383     pdgs[1]=211;
 384     pdgs[0]=321;
 385     Double_t minvD0bar = d->InvMassD0bar();
 386     //apply cut on invariant mass on the pair
 387     if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
 388       //cout<<"inside mass cut"<<endl;
 389       Int_t pdgDgD0toKpi[2]={321,211};
 390       Int_t lab=d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 391
 392       if(lab>=0){ //signal
 393         if(fArray==1) cout<<"LS signal: ERROR"<<endl;
 394         for (Int_t iprong=0; iprong<2; iprong++){
 395           AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(iprong);
 396           Int_t labprong=prong->GetLabel();
 397
 398           //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
 399           AliAODMCParticle *mcprong=0;
 400           if(labprong>=0)  mcprong= (AliAODMCParticle*)mcArray->At(labprong);
 401           Int_t pdgprong=mcprong->GetPdgCode();
 402           if(TMath::Abs(pdgprong)==211) {
 403             //cout<<"pi"<<endl;
 404             ((TH1F*)fDistr->FindObject("hptpiS"))->Fill(d->PtProng(iprong));
 405             ((TH1F*)fDistr->FindObject("hd0piS"))->Fill(d->Getd0Prong(iprong));
 406           }
 407
 408           if(TMath::Abs(pdgprong)==321) {
 409             //cout<<"kappa"<<endl;
 410             ((TH1F*)fDistr->FindObject("hptKS"))->Fill(d->PtProng(iprong));
 411             ((TH1F*)fDistr->FindObject("hd0KS"))->Fill(d->Getd0Prong(iprong));
 412           }
 413           ((TH1F*)fDistr->FindObject("hdcaS"))->Fill(d->GetDCA());
 414
 415         }
 416
 417         if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
 418           ((TH1F*)fDistr->FindObject("hcosthetastarS"))->Fill(d->CosThetaStarD0());
 419         else ((TH1F*)fDistr->FindObject("hcosthetastarS"))->Fill(d->CosThetaStarD0bar());
 420
 421         ((TH1F*)fDistr->FindObject("hd0d0S"))->Fill(d->Prodd0d0());
 422
 423         ((TH1F*)fDistr->FindObject("hcosthetapointS"))->Fill(d->CosPointingAngle());
 424         ((TH1F*)fDistr->FindObject("hcosthpointd0d0S"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 425
 426         //cout<<"ok point"<<endl;
 427
 428
 429       } else{ //Background or LS
 430         //cout<<"is background"<<endl;
 431         AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(0);
 432         if(!prong) cout<<"No daughter found";
 433         else{
 434           if(prong->Charge()==1) {fTotPosPairs[4]++;} else {fTotNegPairs[4]++;}
 435         }
 436         ((TH1F*)fDistr->FindObject("hptB"))->Fill(d->PtProng(0));
 437         //cout<<"ptok"<<endl;
 438         ((TH1F*)fDistr->FindObject("hd0B"))->Fill(d->Getd0Prong(0));
 439         //cout<<"d0ok"<<endl;
 440         ((TH1F*)fDistr->FindObject("hdcaB"))->Fill(d->GetDCA());
 441         //cout<<"dcaok"<<endl;
 442         ((TH1F*)fDistr->FindObject("hcosthetastarB"))->Fill(d->CosThetaStarD0());
 443         ((TH1F*)fDistr->FindObject("hcosthetastarB"))->Fill(d->CosThetaStarD0bar());
 444         ((TH1F*)fDistr->FindObject("hd0d0B"))->Fill(d->Prodd0d0());
 445         //cout<<"d0d0ok"<<endl;
 446         ((TH1F*)fDistr->FindObject("hcosthetapointB"))->Fill(d->CosPointingAngle());
 447         ((TH1F*)fDistr->FindObject("hcosthpointd0d0B"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
 448
 449         //cout<<"pointok"<<endl;
 450
 451       }
 452
 453     } //inv mass cut
 454
 455      //cuts order
 456 //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 457 //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 458 //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 459 //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 460 //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 461 //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 462 //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 463 //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 464 //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 465
 466     Double_t pt = d->Pt();
 467
 468     Int_t ptbin=0;
 469
 470     //cout<<"P_t = "<<pt<<endl;
 471     if (pt>0. && pt<=1.) {
 472       ptbin=1;
 473       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.5);
 474       fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.00025,0.7);
 475 //       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.7);
 476 //       fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,1,1,-0.00015,0.5);
 477       //printf("I'm in the bin %d\n",ptbin);
 478
 479     }
 480
 481     if(pt>1. && pt<=3.) {
 482       if(pt>1. && pt<=2.) ptbin=2;
 483       if(pt>2. && pt<=3.) ptbin=3;
 484       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0002,0.8);
 485       fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
 486       //printf("I'm in the bin %d\n",ptbin);
 487     }
 488
 489     if(pt>3. && pt<=5.){
 490         ptbin=4;
 491         fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0001,0.8);
 492         fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.8);
 493         //printf("I'm in the bin %d\n",ptbin);
 494     }
 495     if(pt>5.){
 496       ptbin=5;
 497       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00005,0.8);
 498       fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.9);
 499     }//if(pt>5)
 500
 501     //printf("I'm in the bin %d\n",ptbin);
 502     //old
 503     //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
 504     if(prongsinacc){
 505       FillHists(ptbin,d,mcArray,fVHFPPR,fOutputPPR);
 506       FillHists(ptbin,d,mcArray,fVHFloose,fOutputloose);
 507     }
 508     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 509   }
 510
 511
 512
 513
 514   // Post the data
 515   PostData(1,fOutputPPR);
 516   PostData(2,fOutputloose);
 517   PostData(4,fDistr);
 518   return;
 519 }
 520 //____________________________________________________________________________*
 521 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAnalysisVertexingHF *vhf, TList *listout){
 522   //
 523   // function used in UserExec:
 524   //
 525
 526
 527   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 528
 529   Int_t okD0=0,okD0bar=0;
 530   //cout<<"inside FillHist"<<endl;
 531   if(part->SelectD0(vhf->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
 532     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 533     //printf("SELECTED\n");
 534
 535
 536     AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
 537     if(!prong) cout<<"No daughter found";
 538     else{
 539       if(prong->Charge()==1) {fTotPosPairs[ptbin-1]++;} else {fTotNegPairs[ptbin-1]++;}
 540     }
 541
 542     TString fillthis="";
 543     Int_t pdgDgD0toKpi[2]={321,211};
 544     Int_t labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
 545
 546     //count candidates selected by cuts
 547     fNentries->Fill(2);
 548     //count true D0 selected by cuts
 549     if (labD0>=0) fNentries->Fill(3);
 550     PostData(3,fNentries);
 551
 552     if (okD0==1) {
 553       fillthis="histMass_";
 554       fillthis+=ptbin;
 555       //cout<<"Filling "<<fillthis<<endl;
 556
 557       //printf("Fill mass with D0");
 558       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 559
 560       if(labD0>=0) {
 561         if(fArray==1) cout<<"LS signal ERROR"<<endl;
 562
 563         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 564         Int_t pdgD0 = partD0->GetPdgCode();
 565         //cout<<"pdg = "<<pdgD0<<endl;
 566         if (pdgD0==421){ //D0
 567           //cout<<"Fill S with D0"<<endl;
 568           fillthis="histSgn_";
 569           fillthis+=ptbin;
 570           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 571           if(TMath::Abs(invmassD0 - mPDG) < 0.027){
 572             fillthis="histSgn27_";
 573             fillthis+=ptbin;
 574             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 575           }
 576         } else{ //it was a D0bar
 577           fillthis="histRfl_";
 578           fillthis+=ptbin;
 579           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 580         }
 581       } else {//background
 582         fillthis="histBkg_";
 583         fillthis+=ptbin;
 584         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 585       }
 586
 587     }
 588     if (okD0bar==1) {
 589       fillthis="histMass_";
 590       fillthis+=ptbin;
 591       //printf("Fill mass with D0bar");
 592       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
 593
 594       if(labD0>=0) {
 595         if(fArray==1) cout<<"LS signal ERROR"<<endl;
 596         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 597         Int_t pdgD0 = partD0->GetPdgCode();
 598         //cout<<" pdg = "<<pdgD0<<endl;
 599         if (pdgD0==-421){ //D0bar
 600           fillthis="histSgn_";
 601           fillthis+=ptbin;
 602           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 603           if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
 604             fillthis="histSgn27_";
 605             fillthis+=ptbin;
 606             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 607           }
 608
 609
 610         } else{
 611           fillthis="histRfl_";
 612           fillthis+=ptbin;
 613           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 614         }
 615       } else {//background or LS
 616         fillthis="histBkg_";
 617         fillthis+=ptbin;
 618         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 619       }
 620     }
 621
 622
 623   } //else cout<<"NOT SELECTED"<<endl;
 624
 625 }
 626 //________________________________________________________________________
 627 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
 628 {
 629   // Terminate analysis
 630   //
 631   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
 632
 633   fOutputPPR = dynamic_cast<TList*> (GetOutputData(1));
 634   if (!fOutputPPR) {
 635     printf("ERROR: fOutputthight not available\n");
 636     return;
 637   }
 638   fOutputloose = dynamic_cast<TList*> (GetOutputData(2));
 639   if (!fOutputloose) {
 640     printf("ERROR: fOutputloose not available\n");
 641     return;
 642   }
 643   fDistr = dynamic_cast<TList*> (GetOutputData(4));
 644   if (!fDistr) {
 645     printf("ERROR: fDistr not available\n");
 646     return;
 647   }
 648
 649   if(fArray==1){
 650     for(Int_t ipt=0;ipt<4;ipt++){
 651       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
 652
 653
 654       if(fLsNormalization>0) {
 655
 656         TString massName="histMass_";
 657         massName+=ipt+1;
 658         ((TH1F*)fOutputPPR->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputPPR->FindObject(massName))->GetEntries());
 659         ((TH1F*)fOutputloose->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputloose->FindObject(massName))->GetEntries());
 660       }
 661     }
 662
 663     fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
 664
 665     if(fLsNormalization>0) {
 666
 667       TString nameDistr="hptB";
 668       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 669       nameDistr="hdcaB";
 670       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 671       nameDistr="hcosthetastarB";
 672       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 673       nameDistr="hd0B";
 674       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 675       nameDistr="hd0d0B";
 676       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 677       nameDistr="hcosthetapointB";
 678       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 679       nameDistr="hcosthpointd0d0B";
 680       ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
 681
 682     }
 683   }
 684
 685   return;
 686 }
 687