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