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