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 <AliAnalysisDataSlot.h>
  36 #include <AliAnalysisDataContainer.h>
  37 #include "AliAnalysisManager.h"
  38 #include "AliESDtrack.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAODRecoDecayHF2Prong.h"
  46 #include "AliAODRecoCascadeHF.h"
  47 #include "AliAnalysisVertexingHF.h"
  48 #include "AliAnalysisTaskSE.h"
  49 #include "AliAnalysisTaskSED0Mass.h"
  50
  51
  52 ClassImp(AliAnalysisTaskSED0Mass)
  53
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  57 AliAnalysisTaskSE(),
  58 fOutputMass(0),
  59 fDistr(0),
  60 fNentries(0),
  61 fChecks(0),
  62 fCutList(0),
  63 fCuts(0),
  64 fArray(0),
  65 fReadMC(0),
  66 fCutOnDistr(0),
  67 fLsNormalization(1.)
  68
  69
  70 {
  71   // Default constructor
  72   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  73 }
  74
  75 //________________________________________________________________________
  76 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
  77 AliAnalysisTaskSE(name),
  78 fOutputMass(0),
  79 fDistr(0),
  80 fNentries(0),
  81 fChecks(0),
  82 fCutList(0),
  83 fCuts(0),
  84 fArray(0),
  85 fReadMC(0),
  86 fCutOnDistr(0),
  87 fLsNormalization(1.)
  88
  89
  90 {
  91   // Default constructor
  92   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
  93
  94   //fCuts=new AliRDHFCutsD0toKpi(*cuts);
  95   fCuts=cuts;
  96
  97   // Output slot #1 writes into a TList container (mass with cuts)
  98   DefineOutput(1,TList::Class());  //My private output
  99   // Output slot #2 writes into a TList container (distributions)
 100   DefineOutput(2,TList::Class());  //My private output
 101   // Output slot #3 writes into a TH1F container (number of events)
 102   DefineOutput(3,TH1F::Class());  //My private output
 103   // Output slot #4 writes into a TList container (quality check)
 104   DefineOutput(4,TList::Class());  //My private output
 105   // Output slot #5 writes into a TList container (cuts)
 106   DefineOutput(5,AliRDHFCutsD0toKpi::Class());  //My private output
 107   //DefineOutput(5,TList::Class());  //My private output
 108
 109 }
 110
 111 //________________________________________________________________________
 112 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
 113 {
 114    if (fOutputMass) {
 115     delete fOutputMass;
 116     fOutputMass = 0;
 117   }
 118   if (fDistr) {
 119     delete fDistr;
 120     fDistr = 0;
 121   }
 122   if (fChecks) {
 123     delete fChecks;
 124     fChecks = 0;
 125   }
 126   if (fCuts) {
 127     delete fCuts;
 128     fCuts = 0;
 129   }
 130   if (fNentries){
 131     delete fNentries;
 132     fNentries = 0;
 133   }
 134
 135   if(fCutList){
 136     delete fCutList;
 137     fCutList=0;
 138   }
 139
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskSED0Mass::Init()
 144 {
 145   // Initialization
 146
 147   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 148
 149
 150   fCutList=new TList();
 151   fCutList->SetOwner();
 152   fCutList->SetName("listofCutsObj");
 153   fCutList->Add(fCuts);
 154
 155   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
 156   // Post the data
 157   PostData(5,copyfCuts);
 158   //PostData(5,fCutList);
 159
 160   return;
 161 }
 162
 163 //________________________________________________________________________
 164 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 165 {
 166
 167   // Create the output container
 168   //
 169   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 170
 171   // Several histograms are more conveniently managed in a TList
 172   fOutputMass = new TList();
 173   fOutputMass->SetOwner();
 174   fOutputMass->SetName("listMass");
 175
 176   fDistr = new TList();
 177   fDistr->SetOwner();
 178   fDistr->SetName("distributionslist");
 179
 180   fChecks = new TList();
 181   fChecks->SetOwner();
 182   fChecks->SetName("checkHistograms");
 183
 184   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
 185
 186   for(Int_t i=0;i<fCuts->GetNPtBins();i++){
 187
 188     nameMass="histMass_";
 189     nameMass+=i;
 190     nameSgn27="histSgn27_";
 191     nameSgn27+=i;
 192     nameSgn="histSgn_";
 193     nameSgn+=i;
 194     nameBkg="histBkg_";
 195     nameBkg+=i;
 196     nameRfl="histRfl_";
 197     nameRfl+=i;
 198     nameMassNocutsS="hMassS_";
 199     nameMassNocutsS+=i;
 200     nameMassNocutsB="hMassB_";
 201     nameMassNocutsB+=i;
 202
 203     //histograms of cut variable distributions
 204
 205     //  pT
 206     namedistr="hptpiS_";
 207     namedistr+=i;
 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;
 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;
 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;
 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;
 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;
 229     TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 230
 231     namedistr="hptB2prongsnoMcut_";
 232     namedistr+=i;
 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;
 239     TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 240     namedistr="hdcaB_";
 241     namedistr+=i;
 242     TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 243
 244     //  costhetastar
 245     namedistr="hcosthetastarS_";
 246     namedistr+=i;
 247     TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 248     namedistr="hcosthetastarB_";
 249     namedistr+=i;
 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;
 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;
 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;
 262     TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution;d0 [cm]",200,-0.1,0.1);
 263
 264     namedistr="hd0d0S_";
 265     namedistr+=i;
 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;
 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;
 274     TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 275     namedistr="hcosthetapointB_";
 276     namedistr+=i;
 277     TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 278
 279     namedistr="hcosthpointd0d0S_";
 280     namedistr+=i;
 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;
 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     fOutputMass->Add(tmpMt);
 355     fOutputMass->Add(tmpSt);
 356     fOutputMass->Add(tmpS27t);
 357     fOutputMass->Add(tmpBt);
 358     fOutputMass->Add(tmpRt);
 359
 360     fDistr->Add(tmpMS);
 361     fDistr->Add(tmpMB);
 362
 363
 364   }
 365
 366   //histograms for vertex checking
 367   TString checkname="hptGoodTr";
 368
 369   TH1F* hptGoodTr=new TH1F(checkname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Number",200,0.,8.);
 370   hptGoodTr->SetTitleOffset(1.3,"Y");
 371   checkname="hdistrGoodTr";
 372
 373   TH1F* hdistrGoodTr=new TH1F(checkname.Data(),"Distribution of number of good tracks per event;no.good-tracks/ev;Entries",31,0,31);
 374   hdistrGoodTr->SetTitleOffset(1.3,"Y");
 375
 376   //conta gli eventi con vertice buoni e almeno due tracce utilizzabili
 377   fChecks->Add(hptGoodTr);
 378   fChecks->Add(hdistrGoodTr);
 379
 380   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 381   cout<<nameoutput<<endl;
 382   fNentries=new TH1F(nameoutput, "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.);
 383
 384   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 385   fNentries->GetXaxis()->SetBinLabel(2,"nCandidatesSelected");
 386   fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
 387   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtx");
 388   fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
 389   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 390
 391
 392   // Post the data
 393   PostData(1,fOutputMass);
 394   PostData(2,fDistr);
 395   PostData(3,fNentries);
 396   PostData(4,fChecks);
 397
 398   return;
 399 }
 400
 401 //________________________________________________________________________
 402 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 403 {
 404   // Execute analysis for current event:
 405   // heavy flavor candidates association to MC truth
 406   //cout<<"I'm in UserExec"<<endl;
 407
 408
 409     //cuts order
 410     //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 411     //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 412     //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 413     //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 414     //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 415     //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 416     //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 417     //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 418     //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 419
 420
 421   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 422
 423   TString bname;
 424   if(fArray==0){ //D0 candidates
 425     // load D0->Kpi candidates
 426     //cout<<"D0 candidates"<<endl;
 427     bname="D0toKpi";
 428   } else { //LikeSign candidates
 429     // load like sign candidates
 430     //cout<<"LS candidates"<<endl;
 431     bname="LikeSign2Prong";
 432   }
 433
 434   TClonesArray *inputArray=0;
 435
 436   if(!aod && AODEvent() && IsStandardAOD()) {
 437     // In case there is an AOD handler writing a standard AOD, use the AOD
 438     // event in memory rather than the input (ESD) event.
 439     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 440     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 441     // have to taken from the AOD event hold by the AliAODExtension
 442     AliAODHandler* aodHandler = (AliAODHandler*)
 443       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 444
 445     if(aodHandler->GetExtensions()) {
 446       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 447       AliAODEvent* aodFromExt = ext->GetAOD();
 448       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 449     }
 450   } else {
 451     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 452   }
 453
 454
 455   if(!inputArray) {
 456     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
 457     return;
 458   }
 459
 460
 461   TClonesArray *mcArray = 0;
 462   AliAODMCHeader *mcHeader = 0;
 463
 464   if(fReadMC) {
 465     // load MC particles
 466     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 467     if(!mcArray) {
 468       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 469       return;
 470     }
 471
 472     // load MC header
 473     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 474     if(!mcHeader) {
 475       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 476       return;
 477     }
 478   }
 479
 480   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
 481
 482   //histogram filled with 1 for every AOD
 483   fNentries->Fill(0);
 484
 485
 486
 487   // AOD primary vertex
 488   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 489
 490   Bool_t isGoodVtx=kFALSE;
 491
 492    //vtx1->Print();
 493   TString primTitle = vtx1->GetTitle();
 494   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 495     isGoodVtx=kTRUE;
 496     fNentries->Fill(3);
 497   }
 498
 499   //cout<<"Start checks"<<endl;
 500   Int_t ntracks=0,isGoodTrack=0;
 501
 502   if(aod) ntracks=aod->GetNTracks();
 503
 504   //cout<<"ntracks = "<<ntracks<<endl;
 505
 506   //loop on tracks in the event
 507   for (Int_t k=0;k<ntracks;k++){
 508     AliAODTrack* track=aod->GetTrack(k);
 509     //cout<<"in loop"<<endl;
 510     //check clusters of the tracks
 511     Int_t nclsTot=0,nclsSPD=0;
 512
 513     for(Int_t l=0;l<6;l++) {
 514       if(TESTBIT(track->GetITSClusterMap(),l)) {
 515         nclsTot++; if(l<2) nclsSPD++;
 516       }
 517     }
 518
 519     if (track->Pt()>0.3 &&
 520         track->GetStatus()&AliESDtrack::kTPCrefit &&
 521         track->GetStatus()&AliESDtrack::kITSrefit &&
 522         nclsTot>3 &&
 523         nclsSPD>0) {//fill hist good tracks
 524       //cout<<"in if"<<endl;
 525       ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
 526       isGoodTrack++;
 527     }
 528     //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
 529     ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
 530   }
 531   //number of events with good vertex and at least 2 good tracks
 532   if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
 533
 534
 535   // loop over candidates
 536   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 537   if(fDebug>1) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 538
 539   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 540     //Int_t nPosPairs=0, nNegPairs=0;
 541     //cout<<"inside the loop"<<endl;
 542     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 543     Bool_t unsetvtx=kFALSE;
 544     if(!d->GetOwnPrimaryVtx()) {
 545       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 546       unsetvtx=kTRUE;
 547     }
 548
 549     //check reco daughter in acceptance
 550     Double_t eta0=d->EtaProng(0);
 551     Double_t eta1=d->EtaProng(1);
 552
 553     if (TMath::Abs(eta0) < 0.9 && TMath::Abs(eta1) < 0.9) {
 554       FillVarHists(d,mcArray,fCuts,fDistr);
 555       FillMassHists(d,mcArray,fCuts,fOutputMass);
 556     }
 557
 558     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 559   } //end for prongs
 560
 561
 562   // Post the data
 563   PostData(1,fOutputMass);
 564   PostData(2,fDistr);
 565   PostData(3,fNentries);
 566   PostData(4,fChecks);
 567
 568   cout<<"Other PostData"<<endl;
 569   return;
 570 }
 571
 572 //____________________________________________________________________________
 573 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
 574   //
 575   // function used in UserExec to fill variable histograms:
 576   //
 577
 578   //add distr here
 579   UInt_t pdgs[2];
 580
 581   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 582   pdgs[0]=211;
 583   pdgs[1]=321;
 584   Double_t minvD0 = part->InvMassD0();
 585   pdgs[1]=211;
 586   pdgs[0]=321;
 587   Double_t minvD0bar = part->InvMassD0bar();
 588   //cout<<"inside mass cut"<<endl;
 589
 590   Int_t pdgDgD0toKpi[2]={321,211};
 591   Int_t lab=-9999;
 592   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)
 593   //Double_t pt = d->Pt(); //mother pt
 594   Bool_t isSelected=kFALSE;
 595
 596   if(fCutOnDistr){
 597     isSelected = cuts->IsSelected(part,AliRDHFCuts::kAll);
 598     if (!isSelected){
 599       //cout<<"Not Selected"<<endl;
 600       return;
 601     }
 602   }
 603
 604   TString fillthispi="",fillthisK="",fillthis="";
 605
 606   Int_t ptbin=cuts->PtBin(part->Pt());
 607   if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
 608     //printf("\nif no cuts or cuts passed\n");
 609     if(lab>=0 && fReadMC){ //signal
 610
 611       //check pdg of the prongs
 612       AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
 613       AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
 614       Int_t labprong[2];
 615       labprong[0]=prong0->GetLabel();
 616       labprong[1]=prong1->GetLabel();
 617       AliAODMCParticle *mcprong=0;
 618       Int_t pdgProng[2]={0,0};
 619       for (Int_t iprong=0;iprong<2;iprong++){
 620         if(labprong[iprong]>=0)  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
 621         pdgProng[iprong]=mcprong->GetPdgCode();
 622       }
 623
 624       //no mass cut ditributions: ptbis
 625
 626       fillthispi="hptpiSnoMcut_";
 627       fillthispi+=ptbin;
 628
 629       fillthisK="hptKSnoMcut_";
 630       fillthisK+=ptbin;
 631
 632       if (TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321){
 633         ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
 634         ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
 635       }else {
 636         if (TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211){
 637           ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
 638           ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
 639         }
 640       }
 641
 642       //no mass cut ditributions: mass
 643       fillthis="hMassS_";
 644       fillthis+=ptbin;
 645
 646       if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
 647         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
 648       }
 649       else { //D0bar
 650         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
 651       }
 652
 653       //apply cut on invariant mass on the pair
 654       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
 655
 656         if(fArray==1) cout<<"LS signal: ERROR"<<endl;
 657         for (Int_t iprong=0; iprong<2; iprong++){
 658           AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
 659           labprong[iprong]=prong->GetLabel();
 660
 661           //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
 662           if(labprong[iprong]>=0)  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
 663           Int_t pdgprong=mcprong->GetPdgCode();
 664
 665           if(TMath::Abs(pdgprong)==211) {
 666             //cout<<"pi"<<endl;
 667
 668             fillthispi="hptpiS_";
 669             fillthispi+=ptbin;
 670             ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
 671
 672             fillthisK="hd0KS_";
 673             fillthisK+=ptbin;
 674             ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
 675           }
 676
 677           if(TMath::Abs(pdgprong)==321) {
 678             //cout<<"kappa"<<endl;
 679
 680             fillthisK="hptKS_";
 681             fillthisK+=ptbin;
 682             ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
 683
 684             fillthisK="hd0KS_";
 685             fillthisK+=ptbin;
 686             ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
 687           }
 688
 689           fillthis="hdcaS_";
 690           fillthis+=ptbin;
 691           ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
 692
 693           fillthis="hcosthetapointS_";
 694           fillthis+=ptbin;
 695           ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
 696
 697           fillthis="hcosthpointd0d0S_";
 698           fillthis+=ptbin;
 699           ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
 700
 701           fillthis="hcosthetastarS_";
 702           fillthis+=ptbin;
 703           if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
 704           else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
 705
 706           fillthis="hd0d0S_";
 707           fillthis+=ptbin;
 708           ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
 709
 710         }
 711       }
 712     } else{ //Background or LS
 713       //cout<<"is background"<<endl;
 714
 715       //no mass cut distributions: mass, ptbis
 716       fillthis="hMassB_";
 717       fillthis+=ptbin;
 718
 719       if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
 720       if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
 721
 722       fillthis="hptB1prongnoMcut_";
 723       fillthis+=ptbin;
 724
 725       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
 726
 727       fillthis="hptB2prongsnoMcut_";
 728       fillthis+=ptbin;
 729       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
 730       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
 731
 732       //apply cut on invariant mass on the pair
 733       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
 734
 735
 736         AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
 737         if(!prong) cout<<"No daughter found";
 738         else{
 739           if(prong->Charge()==1) {fTotPosPairs[4]++;} else {fTotNegPairs[4]++;}
 740         }
 741
 742         //normalise pt distr to half afterwards
 743         fillthis="hptB_";
 744         fillthis+=ptbin;
 745
 746         ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));((TH1F*)listout->FindObject("hptB_1"))->Fill(part->PtProng(1));
 747
 748         fillthis="hd0B_";
 749         fillthis+=ptbin;
 750         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
 751
 752         fillthis="hdcaB_";
 753         fillthis+=ptbin;
 754         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
 755
 756         fillthis="hcosthetastarB_";
 757         fillthis+=ptbin;
 758         if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
 759             if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
 760
 761         fillthis="hd0d0B_";
 762         fillthis+=ptbin;
 763         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
 764
 765         fillthis="hcosthetapointB_";
 766         fillthis+=ptbin;
 767         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
 768
 769         fillthis="hcosthpointd0d0B_";
 770         fillthis+=ptbin;
 771         ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
 772
 773       }
 774     }
 775   }
 776   return;
 777 }
 778 //____________________________________________________________________________
 779
 780 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
 781   //
 782   // function used in UserExec to fill mass histograms:
 783   //
 784
 785
 786   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 787
 788   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kAll); //selected
 789   //cout<<"is selected = "<<isSelected<<endl;
 790
 791   //cout<<"check cuts = "<<endl;
 792   //cuts->PrintAll();
 793   if (!isSelected){
 794     //cout<<"Not Selected"<<endl;
 795     return;
 796   }
 797
 798   cout<<"Candidate selected"<<endl;
 799
 800   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
 801   //printf("SELECTED\n");
 802   Int_t ptbin=cuts->PtBin(part->Pt());
 803
 804   AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
 805   if(!prong) cout<<"No daughter found";
 806   else{
 807     if(prong->Charge()==1) {fTotPosPairs[ptbin-1]++;} else {fTotNegPairs[ptbin-1]++;}
 808   }
 809
 810   TString fillthis="";
 811   Int_t pdgDgD0toKpi[2]={321,211};
 812   Int_t labD0=-1;
 813   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)
 814
 815   //count candidates selected by cuts
 816   fNentries->Fill(1);
 817   //count true D0 selected by cuts
 818   if (fReadMC && labD0>=0) fNentries->Fill(2);
 819   //PostData(3,fNentries);
 820
 821   if (isSelected==1 || isSelected==3) { //D0
 822     fillthis="histMass_";
 823     fillthis+=ptbin;
 824     //cout<<"Filling "<<fillthis<<endl;
 825
 826     //printf("Fill mass with D0");
 827     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 828
 829     if(labD0>=0) {
 830       if(fArray==1) cout<<"LS signal ERROR"<<endl;
 831
 832       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 833       Int_t pdgD0 = partD0->GetPdgCode();
 834       //cout<<"pdg = "<<pdgD0<<endl;
 835       if (pdgD0==421){ //D0
 836         //cout<<"Fill S with D0"<<endl;
 837         fillthis="histSgn_";
 838         fillthis+=ptbin;
 839         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 840         if(TMath::Abs(invmassD0 - mPDG) < 0.027){
 841           fillthis="histSgn27_";
 842           fillthis+=ptbin;
 843           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 844         }
 845       } else{ //it was a D0bar
 846         fillthis="histRfl_";
 847         fillthis+=ptbin;
 848         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 849       }
 850     } else {//background
 851       fillthis="histBkg_";
 852       fillthis+=ptbin;
 853       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
 854     }
 855
 856   }
 857   if (isSelected>1) { //D0bar
 858     fillthis="histMass_";
 859     fillthis+=ptbin;
 860     //printf("Fill mass with D0bar");
 861     ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
 862
 863     if(labD0>=0) {
 864       if(fArray==1) cout<<"LS signal ERROR"<<endl;
 865       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
 866       Int_t pdgD0 = partD0->GetPdgCode();
 867       //cout<<" pdg = "<<pdgD0<<endl;
 868       if (pdgD0==-421){ //D0bar
 869         fillthis="histSgn_";
 870         fillthis+=ptbin;
 871         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 872         if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
 873           fillthis="histSgn27_";
 874           fillthis+=ptbin;
 875           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 876         }
 877
 878
 879       } else{
 880         fillthis="histRfl_";
 881         fillthis+=ptbin;
 882         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 883       }
 884     } else {//background or LS
 885       fillthis="histBkg_";
 886       fillthis+=ptbin;
 887       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
 888     }
 889   }
 890
 891   return;
 892 }
 893 //________________________________________________________________________
 894 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
 895 {
 896   // Terminate analysis
 897   //
 898   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
 899
 900
 901   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
 902   if (!fOutputMass) {
 903     printf("ERROR: fOutputMass not available\n");
 904     return;
 905   }
 906   fDistr = dynamic_cast<TList*> (GetOutputData(2));
 907   if (!fDistr) {
 908     printf("ERROR: fDistr not available\n");
 909     return;
 910   }
 911
 912 //   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
 913 //   if(!fNentries){
 914 //     printf("ERROR: fNEntries not available\n");
 915 //     return;
 916 //   }
 917
 918   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
 919
 920   if(!fNentries){
 921     printf("ERROR: fNEntries not available\n");
 922     return;
 923   }
 924
 925   fChecks = dynamic_cast<TList*> (GetOutputData(4));
 926   if (!fChecks) {
 927     printf("ERROR: fChecks not available\n");
 928     return;
 929   }
 930
 931   if(fArray==1){
 932     for(Int_t ipt=0;ipt<5;ipt++){
 933       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
 934
 935
 936       if(fLsNormalization>1e-6) {
 937
 938         TString massName="histMass_";
 939         massName+=ipt;
 940         ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
 941
 942       }
 943
 944
 945       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
 946
 947       if(fLsNormalization>1e-6) {
 948
 949         TString nameDistr="hptB_";
 950         nameDistr+=ipt;
 951         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 952         nameDistr="hdcaB_";
 953         nameDistr+=ipt;
 954         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 955         nameDistr="hcosthetastarB_";
 956         nameDistr+=ipt;
 957         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 958         nameDistr="hd0B_";
 959         nameDistr+=ipt;
 960         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 961         nameDistr="hd0d0B_";
 962         nameDistr+=ipt;
 963         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 964         nameDistr="hcosthetapointB_";
 965         nameDistr+=ipt;
 966         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
 967         nameDistr="hcosthpointd0d0B_";
 968         nameDistr+=ipt;
 969         ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
 970
 971       }
 972     }
 973   }
 974   TString cvname;
 975
 976   if (fArray==0){
 977     cvname="D0invmass";
 978   } else cvname="LSinvmass";
 979
 980   TCanvas *cMass=new TCanvas(cvname,cvname);
 981   cMass->cd();
 982   ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
 983
 984   TCanvas* cStat=new TCanvas("cstat","Stat");
 985   cStat->cd();
 986   cStat->SetGridy();
 987   //((TH1F*)fDistr->FindObject("nEntriesD0"))->Draw("htext0");
 988   fNentries->Draw("htext0");
 989
 990   return;
 991 }
 992
 993