PWGHF/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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
  21 // and comparison with the MC truth and cut variables distributions.
  22 //
  23 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
  24 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
  25 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
  26 // Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian)
  27 /////////////////////////////////////////////////////////////
  28
  29 #include <Riostream.h>
  30 #include <TClonesArray.h>
  31 #include <TCanvas.h>
  32 #include <TNtuple.h>
  33 #include <TTree.h>
  34 #include <TList.h>
  35 #include <TH1F.h>
  36 #include <TH2F.h>
  37 #include <TDatabasePDG.h>
  38
  39 #include <AliAnalysisDataSlot.h>
  40 #include <AliAnalysisDataContainer.h>
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDtrack.h"
  43 #include "AliVertexerTracks.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODEvent.h"
  46 #include "AliAODVertex.h"
  47 #include "AliAODTrack.h"
  48 #include "AliAODMCHeader.h"
  49 #include "AliAODMCParticle.h"
  50 #include "AliAODRecoDecayHF2Prong.h"
  51 #include "AliAODRecoCascadeHF.h"
  52 #include "AliAnalysisVertexingHF.h"
  53 #include "AliAnalysisTaskSE.h"
  54 #include "AliAnalysisTaskSED0Mass.h"
  55 #include "AliNormalizationCounter.h"
  56
  57 using std::cout;
  58 using std::endl;
  59
  60 ClassImp(AliAnalysisTaskSED0Mass)
  61
  62
  63 //________________________________________________________________________
  64 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
  65 AliAnalysisTaskSE(),
  66   fOutputMass(0),
  67   fOutputMassPt(0),
  68   fDistr(0),
  69   fNentries(0),
  70   fCuts(0),
  71   fArray(0),
  72   fReadMC(0),
  73   fCutOnDistr(0),
  74   fUsePid4Distr(0),
  75   fCounter(0),
  76   fNPtBins(1),
  77   fLsNormalization(1.),
  78   fFillOnlyD0D0bar(0),
  79   fDaughterTracks(),
  80   fIsSelectedCandidate(0),
  81   fFillVarHists(kTRUE),
  82   fSys(0),
  83   fIsRejectSDDClusters(0),
  84   fFillPtHist(kTRUE),
  85   fFillImpParHist(kFALSE),
  86   fUseSelectionBit(kTRUE),
  87   fWriteVariableTree(kFALSE),
  88   fVariablesTree(0),
  89   fCandidateVariables(),
  90   fPIDCheck(kFALSE),
  91   fDrawDetSignal(kFALSE),
  92   fDetSignal(0)
  93 {
  94   // Default constructor
  95   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
  96
  97 }
  98
  99 //________________________________________________________________________
 100 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
 101   AliAnalysisTaskSE(name),
 102   fOutputMass(0),
 103   fOutputMassPt(0),
 104   fDistr(0),
 105   fNentries(0),
 106   fCuts(0),
 107   fArray(0),
 108   fReadMC(0),
 109   fCutOnDistr(0),
 110   fUsePid4Distr(0),
 111   fCounter(0),
 112   fNPtBins(1),
 113   fLsNormalization(1.),
 114   fFillOnlyD0D0bar(0),
 115   fDaughterTracks(),
 116   fIsSelectedCandidate(0),
 117   fFillVarHists(kTRUE),
 118   fSys(0),
 119   fIsRejectSDDClusters(0),
 120   fFillPtHist(kTRUE),
 121   fFillImpParHist(kFALSE),
 122   fUseSelectionBit(kTRUE),
 123   fWriteVariableTree(kFALSE),
 124   fVariablesTree(0),
 125   fCandidateVariables(),
 126   fPIDCheck(kFALSE),
 127   fDrawDetSignal(kFALSE),
 128   fDetSignal(0)
 129 {
 130   // Default constructor
 131
 132   fNPtBins=cuts->GetNPtBins();
 133
 134   fCuts=cuts;
 135   for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
 136
 137   // Output slot #1 writes into a TList container (mass with cuts)
 138   DefineOutput(1,TList::Class());  //My private output
 139   // Output slot #2 writes into a TList container (distributions)
 140   DefineOutput(2,TList::Class());  //My private output
 141   // Output slot #3 writes into a TH1F container (number of events)
 142   DefineOutput(3,TH1F::Class());  //My private output
 143   // Output slot #4 writes into a TList container (cuts)
 144   DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
 145   // Output slot #5 writes Normalization Counter
 146   DefineOutput(5,AliNormalizationCounter::Class());
 147   // Output slot #6 stores the mass vs pt and impact parameter distributions
 148   DefineOutput(6,TList::Class());  //My private output
 149   // Output slot #7 keeps a tree of the candidate variables after track selection
 150   DefineOutput(7,TTree::Class());  //My private outpu
 151   // Output slot #8 writes into a TList container (Detector signals)
 152   DefineOutput(8, TList::Class()); //My private output
 153 }
 154
 155 //________________________________________________________________________
 156 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
 157 {
 158   if (fOutputMass) {
 159     delete fOutputMass;
 160     fOutputMass = 0;
 161   }
 162   if (fOutputMassPt) {
 163     delete fOutputMassPt;
 164     fOutputMassPt = 0;
 165   }
 166   if (fDistr) {
 167     delete fDistr;
 168     fDistr = 0;
 169   }
 170   if (fCuts) {
 171     delete fCuts;
 172     fCuts = 0;
 173   }
 174   for(Int_t i=0; i<5; i++){
 175     if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
 176   }
 177   if (fNentries){
 178     delete fNentries;
 179     fNentries = 0;
 180   }
 181   if(fCounter){
 182     delete fCounter;
 183     fCounter=0;
 184   }
 185   if(fVariablesTree){
 186     delete fVariablesTree;
 187     fVariablesTree = 0;
 188   }
 189   if (fDetSignal) {
 190     delete fDetSignal;
 191     fDetSignal = 0;
 192   }
 193
 194 }
 195
 196 //________________________________________________________________________
 197 void AliAnalysisTaskSED0Mass::Init()
 198 {
 199   // Initialization
 200
 201   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
 202
 203
 204   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
 205   const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 206   copyfCuts->SetName(nameoutput);
 207   // Post the data
 208   PostData(4,copyfCuts);
 209
 210
 211   return;
 212 }
 213
 214 //________________________________________________________________________
 215 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
 216 {
 217
 218   // Create the output container
 219   //
 220   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
 221
 222   // Several histograms are more conveniently managed in a TList
 223   fOutputMass = new TList();
 224   fOutputMass->SetOwner();
 225   fOutputMass->SetName("listMass");
 226
 227   fOutputMassPt = new TList();
 228   fOutputMassPt->SetOwner();
 229   fOutputMassPt->SetName("listMassPt");
 230
 231   fDistr = new TList();
 232   fDistr->SetOwner();
 233   fDistr->SetName("distributionslist");
 234
 235   fDetSignal = new TList();
 236   fDetSignal->SetOwner();
 237   fDetSignal->SetName("detectorsignals");
 238
 239   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
 240   TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
 241   Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
 242
 243   for(Int_t i=0;i<fCuts->GetNPtBins();i++){
 244
 245     nameMass="histMass_";
 246     nameMass+=i;
 247     nameSgn27="histSgn27_";
 248     nameSgn27+=i;
 249     nameSgn="histSgn_";
 250     nameSgn+=i;
 251     nameBkg="histBkg_";
 252     nameBkg+=i;
 253     nameRfl="histRfl_";
 254     nameRfl+=i;
 255     nameMassNocutsS="hMassS_";
 256     nameMassNocutsS+=i;
 257     nameMassNocutsB="hMassB_";
 258     nameMassNocutsB+=i;
 259
 260     //histograms of cut variable distributions
 261     if(fFillVarHists){
 262       if(fReadMC){
 263         //  dca
 264         namedistr="hdcaS_";
 265         namedistr+=i;
 266         TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 267         // impact parameter
 268         namedistr="hd0piS_";
 269         namedistr+=i;
 270         TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
 271
 272         namedistr="hd0KS_";
 273         namedistr+=i;
 274         TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
 275         namedistr="hd0d0S_";
 276         namedistr+=i;
 277         TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 278
 279         //decay lenght
 280         namedistr="hdeclS_";
 281         namedistr+=i;
 282         TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
 283
 284         namedistr="hnormdeclS_";
 285         namedistr+=i;
 286         TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
 287
 288         namedistr="hdeclxyS_";
 289         namedistr+=i;
 290         TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 291         namedistr="hnormdeclxyS_";
 292         namedistr+=i;
 293         TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 294
 295         namedistr="hdeclxyd0d0S_";
 296         namedistr+=i;
 297         TH2F* hdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation decay Length XY - d_{0}#times d_{0};Decay Length XY [cm];d_{0}#times d_{0}[cm^{2}]",200,0,0.15,200,-0.001,0.001);
 298
 299         namedistr="hnormdeclxyd0d0S_";
 300         namedistr+=i;
 301         TH2F* hnormdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
 302
 303         //  costhetapoint
 304         namedistr="hcosthetapointS_";
 305         namedistr+=i;
 306         TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 307
 308         namedistr="hcosthetapointxyS_";
 309         namedistr+=i;
 310         TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
 311
 312         TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 313         tmpMS->Sumw2();
 314
 315
 316         fDistr->Add(hdcaS);
 317
 318         fDistr->Add(hd0piS);
 319         fDistr->Add(hd0KS);
 320
 321         fDistr->Add(hd0d0S);
 322
 323         fDistr->Add(hcosthetapointS);
 324
 325         fDistr->Add(hcosthetapointxyS);
 326
 327         fDistr->Add(hdeclengthS);
 328
 329         fDistr->Add(hnormdeclengthS);
 330
 331         fDistr->Add(hdeclxyS);
 332
 333         fDistr->Add(hnormdeclxyS);
 334
 335         fDistr->Add(hdeclxyd0d0S);
 336         fDistr->Add(hnormdeclxyd0d0S);
 337
 338         fDistr->Add(tmpMS);
 339       }
 340
 341       //  dca
 342       namedistr="hdcaB_";
 343       namedistr+=i;
 344       TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
 345
 346       // impact parameter
 347       namedistr="hd0B_";
 348       namedistr+=i;
 349       TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 350
 351       namedistr="hd0d0B_";
 352       namedistr+=i;
 353       TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 354
 355       //decay lenght
 356       namedistr="hdeclB_";
 357       namedistr+=i;
 358       TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
 359
 360       namedistr="hnormdeclB_";
 361       namedistr+=i;
 362       TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
 363
 364       namedistr="hdeclxyB_";
 365       namedistr+=i;
 366       TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
 367       namedistr="hnormdeclxyB_";
 368       namedistr+=i;
 369       TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
 370
 371       namedistr="hdeclxyd0d0B_";
 372       namedistr+=i;
 373       TH2F* hdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation decay Length XY - d_{0}#times d_{0};Decay Length XY [cm];d_{0}#times d_{0}[cm^{2}]",200,0,0.15,200,-0.001,0.001);
 374
 375       namedistr="hnormdeclxyd0d0B_";
 376       namedistr+=i;
 377       TH2F* hnormdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
 378
 379       //  costhetapoint
 380       namedistr="hcosthetapointB_";
 381       namedistr+=i;
 382       TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 383
 384       namedistr="hcosthetapointxyB_";
 385       namedistr+=i;
 386       TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
 387
 388       TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
 389       tmpMB->Sumw2();
 390
 391
 392
 393       fDistr->Add(hdcaB);
 394
 395       fDistr->Add(hd0B);
 396
 397       fDistr->Add(hd0d0B);
 398
 399       fDistr->Add(hcosthetapointB);
 400
 401       fDistr->Add(hcosthetapointxyB);
 402
 403       fDistr->Add(hdeclengthB);
 404
 405       fDistr->Add(hnormdeclengthB);
 406
 407       fDistr->Add(hdeclxyB);
 408
 409       fDistr->Add(hnormdeclxyB);
 410
 411       fDistr->Add(hdeclxyd0d0B);
 412       fDistr->Add(hnormdeclxyd0d0B);
 413
 414       fDistr->Add(tmpMB);
 415
 416       //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
 417
 418       if(fCuts->GetIsPrimaryWithoutDaughters()){
 419         if(fReadMC){
 420           namedistr="hd0vpiS_";
 421           namedistr+=i;
 422           TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
 423
 424           namedistr="hd0vKS_";
 425           namedistr+=i;
 426           TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
 427
 428           namedistr="hd0d0vS_";
 429           namedistr+=i;
 430           TH1F *hd0d0vS = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 431
 432           namedistr="hdeclvS_";
 433           namedistr+=i;
 434           TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 435
 436           namedistr="hnormdeclvS_";
 437           namedistr+=i;
 438           TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 439
 440           fDistr->Add(hd0vpiS);
 441           fDistr->Add(hd0vKS);
 442           fDistr->Add(hd0d0vS);
 443           fDistr->Add(hdeclengthvS);
 444           fDistr->Add(hnormdeclengthvS);
 445
 446         }
 447
 448         namedistr="hd0vmoresB_";
 449         namedistr+=i;
 450         TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 451
 452         namedistr="hd0d0vmoresB_";
 453         namedistr+=i;
 454         TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 455
 456
 457         namedistr="hd0vB_";
 458         namedistr+=i;
 459         TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 460
 461         namedistr="hd0vp0B_";
 462         namedistr+=i;
 463         TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 464
 465         namedistr="hd0vp1B_";
 466         namedistr+=i;
 467         TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
 468
 469         namedistr="hd0d0vB_";
 470         namedistr+=i;
 471         TH1F *hd0d0vB = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
 472
 473         namedistr="hdeclvB_";
 474         namedistr+=i;
 475         TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
 476
 477         namedistr="hnormdeclvB_";
 478         namedistr+=i;
 479         TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
 480
 481         fDistr->Add(hd0vB);
 482         fDistr->Add(hd0vp0B);
 483         fDistr->Add(hd0vp1B);
 484         fDistr->Add(hd0vmoresB);
 485
 486         fDistr->Add(hd0d0vB);
 487         fDistr->Add(hd0d0vmoresB);
 488
 489         fDistr->Add(hdeclengthvB);
 490
 491         fDistr->Add(hnormdeclengthvB);
 492       }
 493
 494
 495     }
 496     //histograms of invariant mass distributions
 497
 498
 499     //MC signal
 500     if(fReadMC){
 501       TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 502
 503       TH1F *tmpSl=(TH1F*)tmpSt->Clone();
 504       tmpSt->Sumw2();
 505       tmpSl->Sumw2();
 506
 507       //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 508       TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 509       //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
 510       TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
 511       //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
 512       tmpBt->Sumw2();
 513       //tmpBl->Sumw2();
 514       tmpRt->Sumw2();
 515       //tmpRl->Sumw2();
 516
 517       fOutputMass->Add(tmpSt);
 518       fOutputMass->Add(tmpRt);
 519       fOutputMass->Add(tmpBt);
 520
 521     }
 522
 523     //mass
 524     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.5648,2.1648);
 525     //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
 526     tmpMt->Sumw2();
 527     //tmpMl->Sumw2();
 528     //distribution w/o cuts large range
 529     // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
 530
 531     fOutputMass->Add(tmpMt);
 532
 533     if(fSys==0){ //histograms filled only in pp to save time in PbPb
 534       if(fFillVarHists){
 535
 536         if(fReadMC){
 537           //  pT
 538           namedistr="hptpiS_";
 539           namedistr+=i;
 540           TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 541
 542           namedistr="hptKS_";
 543           namedistr+=i;
 544           TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 545
 546           //  costhetastar
 547           namedistr="hcosthetastarS_";
 548           namedistr+=i;
 549           TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 550
 551           //pT no mass cut
 552
 553           namedistr="hptpiSnoMcut_";
 554           namedistr+=i;
 555           TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
 556
 557           namedistr="hptKSnoMcut_";
 558           namedistr+=i;
 559           TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
 560
 561           fDistr->Add(hptpiS);
 562           fDistr->Add(hptKS);
 563           fDistr->Add(hcosthetastarS);
 564
 565           fDistr->Add(hptpiSnoMcut);
 566           fDistr->Add(hptKSnoMcut);
 567
 568           //  costhetapoint vs d0 or d0d0
 569           namedistr="hcosthpointd0S_";
 570           namedistr+=i;
 571           TH2F *hcosthpointd0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 572           namedistr="hcosthpointd0d0S_";
 573           namedistr+=i;
 574           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);
 575
 576           fDistr->Add(hcosthpointd0S);
 577           fDistr->Add(hcosthpointd0d0S);
 578
 579           //to compare with AliAnalysisTaskCharmFraction
 580           TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.5648,2.1648);
 581           TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
 582           tmpS27t->Sumw2();
 583           tmpS27l->Sumw2();
 584
 585           fOutputMass->Add(tmpS27t);
 586           fOutputMass->Add(tmpS27l);
 587
 588         }
 589
 590         //  pT
 591         namedistr="hptB_";
 592         namedistr+=i;
 593         TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 594
 595         //  costhetastar
 596         namedistr="hcosthetastarB_";
 597         namedistr+=i;
 598         TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
 599
 600         //pT no mass cut
 601         namedistr="hptB1prongnoMcut_";
 602         namedistr+=i;
 603         TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 604
 605         namedistr="hptB2prongsnoMcut_";
 606         namedistr+=i;
 607         TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
 608
 609         fDistr->Add(hptB);
 610         fDistr->Add(hcosthetastarB);
 611
 612         fDistr->Add(hptB1pnoMcut);
 613         fDistr->Add(hptB2pnoMcut);
 614
 615         //impact parameter of negative/positive track
 616         namedistr="hd0p0B_";
 617         namedistr+=i;
 618         TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
 619
 620         namedistr="hd0p1B_";
 621         namedistr+=i;
 622         TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
 623
 624         //impact parameter corrected for strangeness
 625         namedistr="hd0moresB_";
 626         namedistr+=i;
 627         TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
 628
 629         namedistr="hd0d0moresB_";
 630         namedistr+=i;
 631         TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
 632
 633
 634         namedistr="hcosthetapointmoresB_";
 635         namedistr+=i;
 636         TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
 637
 638         //  costhetapoint vs d0 or d0d0
 639         namedistr="hcosthpointd0B_";
 640         namedistr+=i;
 641         TH2F *hcosthpointd0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
 642
 643         namedistr="hcosthpointd0d0B_";
 644         namedistr+=i;
 645         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);
 646
 647         fDistr->Add(hd0p0B);
 648         fDistr->Add(hd0p1B);
 649
 650         fDistr->Add(hd0moresB);
 651         fDistr->Add(hd0d0moresB);
 652         fDistr->Add(hcosthetapointmoresB);
 653
 654
 655         fDistr->Add(hcosthpointd0B);
 656
 657
 658         fDistr->Add(hcosthpointd0d0B);
 659       }
 660
 661     } //end pp histo
 662   }
 663
 664
 665   //for Like sign analysis
 666
 667   if(fArray==1){
 668     namedistr="hpospair";
 669     TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 670     namedistr="hnegpair";
 671     TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
 672     fOutputMass->Add(hpospair);
 673     fOutputMass->Add(hnegpair);
 674   }
 675
 676
 677   // 2D Pt distributions and impact parameter histograms
 678   if(fFillPtHist) {
 679
 680     nameMassPt="histMassPt";
 681     nameSgnPt="histSgnPt";
 682     nameBkgPt="histBkgPt";
 683     nameRflPt="histRflPt";
 684
 685     //MC signal
 686     if(fReadMC){
 687       TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.16484,nbins2dPt,binInPt,binFinPt);
 688       TH2F *tmpSlPt=(TH2F*)tmpStPt->Clone();
 689       tmpStPt->Sumw2();
 690       tmpSlPt->Sumw2();
 691
 692       //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
 693       TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 694       TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 695       tmpBtPt->Sumw2();
 696       tmpRtPt->Sumw2();
 697
 698       fOutputMassPt->Add(tmpStPt);
 699       fOutputMassPt->Add(tmpRtPt);
 700       fOutputMassPt->Add(tmpBtPt);
 701
 702       //       cout<<endl<<endl<<"***************************************"<<endl;
 703       //       cout << " created and added to the list "<< nameSgnPt.Data() <<" "<< tmpStPt <<
 704       //        ", "<<nameRflPt.Data() <<" " << tmpRtPt<<", "<<nameBkgPt.Data()<< " " << tmpBtPt <<endl;
 705       //       cout<<"***************************************"<<endl<<endl;
 706     }
 707
 708     TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",200,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
 709     tmpMtPt->Sumw2();
 710
 711     fOutputMassPt->Add(tmpMtPt);
 712   }
 713
 714   if(fFillImpParHist) CreateImpactParameterHistos();
 715
 716   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
 717
 718   fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 18,-0.5,17.5);
 719
 720   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 721   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
 722   if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
 723   else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
 724   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
 725   fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
 726   fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
 727   if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
 728   if(fFillVarHists || fPIDCheck){
 729     fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
 730     fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
 731     fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
 732     fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
 733   }
 734   if(fReadMC && fSys==0){
 735     fNentries->GetXaxis()->SetBinLabel(12,"K");
 736     fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
 737   }
 738   fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
 739   fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
 740   if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
 741   if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
 742   fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
 743   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
 744
 745   fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
 746   fCounter->Init();
 747
 748   //
 749   // Output slot 7 : tree of the candidate variables after track selection
 750   //
 751   nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
 752   fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
 753   Int_t nVar = 15;
 754   fCandidateVariables = new Double_t [nVar];
 755   TString * fCandidateVariableNames = new TString[nVar];
 756   fCandidateVariableNames[0] = "massD0";
 757   fCandidateVariableNames[1] = "massD0bar";
 758   fCandidateVariableNames[2] = "pt";
 759   fCandidateVariableNames[3] = "dca";
 760   fCandidateVariableNames[4] = "costhsD0";
 761   fCandidateVariableNames[5] = "costhsD0bar";
 762   fCandidateVariableNames[6] = "ptk";
 763   fCandidateVariableNames[7] = "ptpi";
 764   fCandidateVariableNames[8] = "d0k";
 765   fCandidateVariableNames[9] = "d0pi";
 766   fCandidateVariableNames[10] = "d0xd0";
 767   fCandidateVariableNames[11] = "costhp";
 768   fCandidateVariableNames[12] = "costhpxy";
 769   fCandidateVariableNames[13] = "lxy";
 770   fCandidateVariableNames[14] = "specialcuts";
 771   for(Int_t ivar=0; ivar<nVar; ivar++){
 772     fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
 773   }
 774
 775   //
 776   // Output slot 8 : List for detector response histograms
 777   //
 778   if (fDrawDetSignal) {
 779     TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
 780     TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
 781
 782     TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
 783     TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
 784
 785     fDetSignal->Add(TOFSigBefPID);
 786     fDetSignal->Add(TOFSigAftPID);
 787     fDetSignal->Add(TPCSigBefPID);
 788     fDetSignal->Add(TPCSigAftPID);
 789   }
 790
 791   // Post the data
 792   PostData(1,fOutputMass);
 793   PostData(2,fDistr);
 794   PostData(3,fNentries);
 795   PostData(5,fCounter);
 796   PostData(6,fOutputMassPt);
 797   PostData(7,fVariablesTree);
 798   PostData(8, fDetSignal);
 799
 800   return;
 801 }
 802
 803 //________________________________________________________________________
 804 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
 805 {
 806   // Execute analysis for current event:
 807   // heavy flavor candidates association to MC truth
 808   //cout<<"I'm in UserExec"<<endl;
 809
 810
 811   //cuts order
 812   //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
 813   //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
 814   //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
 815   //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
 816   //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
 817   //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
 818   //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
 819   //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
 820   //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
 821
 822
 823   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 824
 825   TString bname;
 826   if(fArray==0){ //D0 candidates
 827     // load D0->Kpi candidates
 828     //cout<<"D0 candidates"<<endl;
 829     bname="D0toKpi";
 830   } else { //LikeSign candidates
 831     // load like sign candidates
 832     //cout<<"LS candidates"<<endl;
 833     bname="LikeSign2Prong";
 834   }
 835
 836   TClonesArray *inputArray=0;
 837
 838   if(!aod && AODEvent() && IsStandardAOD()) {
 839     // In case there is an AOD handler writing a standard AOD, use the AOD
 840     // event in memory rather than the input (ESD) event.
 841     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 842     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 843     // have to taken from the AOD event hold by the AliAODExtension
 844     AliAODHandler* aodHandler = (AliAODHandler*)
 845       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 846
 847     if(aodHandler->GetExtensions()) {
 848       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 849       AliAODEvent* aodFromExt = ext->GetAOD();
 850       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
 851     }
 852   } else if(aod) {
 853     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
 854   }
 855
 856
 857   if(!inputArray || !aod) {
 858     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
 859     return;
 860   }
 861
 862   // fix for temporary bug in ESDfilter
 863   // the AODs with null vertex pointer didn't pass the PhysSel
 864   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 865
 866
 867   TClonesArray *mcArray = 0;
 868   AliAODMCHeader *mcHeader = 0;
 869
 870   if(fReadMC) {
 871     // load MC particles
 872     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 873     if(!mcArray) {
 874       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
 875       return;
 876     }
 877
 878     // load MC header
 879     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 880     if(!mcHeader) {
 881       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
 882       return;
 883     }
 884   }
 885
 886   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
 887
 888   //histogram filled with 1 for every AOD
 889   fNentries->Fill(0);
 890   fCounter->StoreEvent(aod,fCuts,fReadMC);
 891   //fCounter->StoreEvent(aod,fReadMC);
 892
 893   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 894   TString trigclass=aod->GetFiredTriggerClasses();
 895   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
 896
 897   if(!fCuts->IsEventSelected(aod)) {
 898     if(fCuts->GetWhyRejection()==1) // rejected for pileup
 899       fNentries->Fill(13);
 900     if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
 901     if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
 902     return;
 903   }
 904
 905   // Check the Nb of SDD clusters
 906   if (fIsRejectSDDClusters) {
 907     Bool_t skipEvent = kFALSE;
 908     Int_t ntracks = 0;
 909     if (aod) ntracks = aod->GetNTracks();
 910     for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
 911       //    ... get the track
 912       AliAODTrack * track = aod->GetTrack(itrack);
 913       if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
 914         skipEvent=kTRUE;
 915         fNentries->Fill(16);
 916         break;
 917       }
 918     }
 919     if (skipEvent) return;
 920   }
 921
 922   // AOD primary vertex
 923   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 924
 925   Bool_t isGoodVtx=kFALSE;
 926
 927   //vtx1->Print();
 928   TString primTitle = vtx1->GetTitle();
 929   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 930     isGoodVtx=kTRUE;
 931     fNentries->Fill(3);
 932   }
 933
 934   // loop over candidates
 935   Int_t nInD0toKpi = inputArray->GetEntriesFast();
 936   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
 937
 938   // FILE *f=fopen("4display.txt","a");
 939   // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
 940   Int_t nSelectedloose=0,nSelectedtight=0;
 941   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
 942     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
 943
 944     if(fUseSelectionBit && d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
 945         fNentries->Fill(2);
 946         continue; //skip the D0 from Dstar
 947       }
 948
 949     // Bool_t unsetvtx=kFALSE;
 950     // if(!d->GetOwnPrimaryVtx()) {
 951     //   d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 952     //   unsetvtx=kTRUE;
 953     // }
 954
 955
 956     if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
 957       nSelectedloose++;
 958       nSelectedtight++;
 959       if(fSys==0){
 960         if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
 961       }
 962       Int_t ptbin=fCuts->PtBin(d->Pt());
 963       if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
 964       fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
 965       if(fFillVarHists) {
 966         //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
 967         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
 968         fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
 969         //check daughters
 970         if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
 971           AliDebug(1,"at least one daughter not found!");
 972           fNentries->Fill(5);
 973           fDaughterTracks.Clear();
 974           continue;
 975         }
 976         //}
 977         FillVarHists(aod,d,mcArray,fCuts,fDistr);
 978       }
 979
 980       if (fDrawDetSignal) {
 981         DrawDetSignal(d, fDetSignal);
 982       }
 983
 984       FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
 985       if (fPIDCheck) {
 986         Int_t isSelectedPIDfill = 3;
 987         if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both
 988
 989         if (isSelectedPIDfill == 0)fNentries->Fill(7);
 990         if (isSelectedPIDfill == 1)fNentries->Fill(8);
 991         if (isSelectedPIDfill == 2)fNentries->Fill(9);
 992         if (isSelectedPIDfill == 3)fNentries->Fill(10);
 993       }
 994
 995     }
 996
 997     fDaughterTracks.Clear();
 998     //if(unsetvtx) d->UnsetOwnPrimaryVtx();
 999   } //end for prongs
1000   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1001   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1002   // Post the data
1003   PostData(1,fOutputMass);
1004   PostData(2,fDistr);
1005   PostData(3,fNentries);
1006   PostData(5,fCounter);
1007   PostData(6,fOutputMassPt);
1008   PostData(7,fVariablesTree);
1009   PostData(8, fDetSignal);
1010
1011   return;
1012 }
1013
1014 //____________________________________________________________________________
1015 void AliAnalysisTaskSED0Mass::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
1016 {
1017   //
1018   // Function called in UserExec for drawing detector signal histograms:
1019   //
1020   fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
1021   fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
1022
1023   AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1024   AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1025
1026
1027   //For filling detector histograms
1028   Int_t isSelectedPIDfill = 3;
1029   if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1030
1031
1032   //fill "before PID" histos with every daughter
1033   ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1034   ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1035   ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1036   ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1037
1038   if (isSelectedPIDfill != 0)  { //fill "After PID" with everything that's not rejected
1039     ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1040     ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1041     ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1042     ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1043
1044   }
1045
1046   delete esdtrack1;
1047   delete esdtrack2;
1048
1049   return;
1050 }
1051
1052 //____________________________________________________________________________
1053 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
1054   //
1055   // function used in UserExec to fill variable histograms:
1056   //
1057
1058
1059   Int_t pdgDgD0toKpi[2]={321,211};
1060   Int_t lab=-9999;
1061   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)
1062   //Double_t pt = d->Pt(); //mother pt
1063   Int_t isSelectedPID=3;
1064   if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1065   if (!fPIDCheck) {  //if fPIDCheck, this is already filled elsewhere
1066     if (isSelectedPID==0)fNentries->Fill(7);
1067     if (isSelectedPID==1)fNentries->Fill(8);
1068     if (isSelectedPID==2)fNentries->Fill(9);
1069     if (isSelectedPID==3)fNentries->Fill(10);
1070     //fNentries->Fill(8+isSelectedPID);
1071   }
1072
1073   if(fCutOnDistr && !fIsSelectedCandidate) return;
1074   //printf("\nif no cuts or cuts passed\n");
1075
1076
1077   //add distr here
1078   UInt_t pdgs[2];
1079
1080   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1081   pdgs[0]=211;
1082   pdgs[1]=321;
1083   Double_t minvD0 = part->InvMassD0();
1084   pdgs[1]=211;
1085   pdgs[0]=321;
1086   Double_t minvD0bar = part->InvMassD0bar();
1087   //cout<<"inside mass cut"<<endl;
1088
1089   Double_t invmasscut=0.03;
1090
1091   TString fillthispi="",fillthisK="",fillthis="", fillthispt="";
1092
1093   Int_t ptbin=cuts->PtBin(part->Pt());
1094   Double_t pt = part->Pt();
1095
1096   Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1097   dz1[0]=-99; dz2[0]=-99;
1098   Double_t d0[2];
1099   Double_t decl[2]={-99,-99};
1100   Bool_t recalcvtx=kFALSE;
1101
1102
1103
1104   if(fCuts->GetIsPrimaryWithoutDaughters()){
1105     recalcvtx=kTRUE;
1106     //recalculate vertex
1107     AliAODVertex *vtxProp=0x0;
1108     vtxProp=GetPrimaryVtxSkipped(aod);
1109     if(vtxProp) {
1110       part->SetOwnPrimaryVtx(vtxProp);
1111       //Bool_t unsetvtx=kTRUE;
1112       //Calculate d0 for daughter tracks with Vtx Skipped
1113       AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1114       AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1115       esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1116       esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1117       delete vtxProp; vtxProp=NULL;
1118       delete esdtrack1;
1119       delete esdtrack2;
1120     }
1121
1122     d0[0]=dz1[0];
1123     d0[1]=dz2[0];
1124
1125     decl[0]=part->DecayLength2();
1126     decl[1]=part->NormalizedDecayLength2();
1127     part->UnsetOwnPrimaryVtx();
1128
1129   }
1130
1131   Double_t cosThetaStarD0 = 99;
1132   Double_t cosThetaStarD0bar = 99;
1133   Double_t cosPointingAngle = 99;
1134   Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1135   Double_t decayLength2 = -1, decayLengthxy=-1;
1136   Double_t ptProng[2]={-99,-99};
1137   Double_t d0Prong[2]={-99,-99};
1138
1139
1140   //disable the PID
1141   if(!fUsePid4Distr) isSelectedPID=0;
1142   if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
1143
1144     //check pdg of the prongs
1145     AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1146     AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1147
1148     if(!prong0 || !prong1) {
1149       return;
1150     }
1151     Int_t labprong[2];
1152     if(fReadMC){
1153       labprong[0]=prong0->GetLabel();
1154       labprong[1]=prong1->GetLabel();
1155     }
1156     AliAODMCParticle *mcprong=0;
1157     Int_t pdgProng[2]={0,0};
1158
1159     for (Int_t iprong=0;iprong<2;iprong++){
1160       if(fReadMC && labprong[iprong]>=0) {
1161         mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1162         pdgProng[iprong]=mcprong->GetPdgCode();
1163       }
1164     }
1165
1166     if(fSys==0){
1167       //no mass cut ditributions: ptbis
1168
1169       fillthispi="hptpiSnoMcut_";
1170       fillthispi+=ptbin;
1171
1172       fillthisK="hptKSnoMcut_";
1173       fillthisK+=ptbin;
1174
1175       if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
1176           || (isSelectedPID==1 || isSelectedPID==3)){
1177         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
1178         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
1179       }
1180
1181       if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
1182           || (isSelectedPID==2 || isSelectedPID==3)){
1183         ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
1184         ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
1185       }
1186     }
1187
1188     //no mass cut ditributions: mass
1189     fillthis="hMassS_";
1190     fillthis+=ptbin;
1191     fillthispt="histSgnPt";
1192
1193     if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1194         || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1195       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1196       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1197     }
1198     else { //D0bar
1199       if(fReadMC || (!fReadMC && isSelectedPID > 1)){
1200         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1201         if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1202       }
1203     }
1204
1205     //apply cut on invariant mass on the pair
1206     if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1207
1208       cosThetaStarD0 = part->CosThetaStarD0();
1209       cosThetaStarD0bar = part->CosThetaStarD0bar();
1210       cosPointingAngle = part->CosPointingAngle();
1211       normalizedDecayLength2 = part->NormalizedDecayLength2();
1212       decayLength2 = part->DecayLength2();
1213       decayLengthxy = part->DecayLengthXY();
1214       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1215
1216       ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
1217       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1218
1219       if(fArray==1) cout<<"LS signal: ERROR"<<endl;
1220       for (Int_t iprong=0; iprong<2; iprong++){
1221         AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1222         if (fReadMC) labprong[iprong]=prong->GetLabel();
1223
1224         //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
1225         Int_t pdgprong=0;
1226         if(fReadMC && labprong[iprong]>=0) {
1227           mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1228           pdgprong=mcprong->GetPdgCode();
1229         }
1230
1231         Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
1232
1233         if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1234           //cout<<"pi"<<endl;
1235
1236           if(fSys==0){
1237             fillthispi="hptpiS_";
1238             fillthispi+=ptbin;
1239             ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
1240           }
1241           fillthispi="hd0piS_";
1242           fillthispi+=ptbin;
1243           ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1244           if(recalcvtx) {
1245
1246             fillthispi="hd0vpiS_";
1247             fillthispi+=ptbin;
1248             ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
1249           }
1250
1251         }
1252
1253         if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1254           //cout<<"kappa"<<endl;
1255           if(fSys==0){
1256             fillthisK="hptKS_";
1257             fillthisK+=ptbin;
1258             ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
1259           }
1260           fillthisK="hd0KS_";
1261           fillthisK+=ptbin;
1262           ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1263           if (recalcvtx){
1264             fillthisK="hd0vKS_";
1265             fillthisK+=ptbin;
1266             ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
1267           }
1268         }
1269
1270         if(fSys==0){
1271           fillthis="hcosthpointd0S_";
1272           fillthis+=ptbin;
1273           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1274         }
1275       } //end loop on prongs
1276
1277       fillthis="hdcaS_";
1278       fillthis+=ptbin;
1279       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1280
1281       fillthis="hcosthetapointS_";
1282       fillthis+=ptbin;
1283       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1284
1285       fillthis="hcosthetapointxyS_";
1286       fillthis+=ptbin;
1287       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1288
1289
1290       fillthis="hd0d0S_";
1291       fillthis+=ptbin;
1292       ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1293
1294       fillthis="hdeclS_";
1295       fillthis+=ptbin;
1296       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1297
1298       fillthis="hnormdeclS_";
1299       fillthis+=ptbin;
1300       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1301
1302       fillthis="hdeclxyS_";
1303       fillthis+=ptbin;
1304       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1305
1306       fillthis="hnormdeclxyS_";
1307       fillthis+=ptbin;
1308       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1309
1310       fillthis="hdeclxyd0d0S_";
1311       fillthis+=ptbin;
1312       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1313
1314       fillthis="hnormdeclxyd0d0S_";
1315       fillthis+=ptbin;
1316       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1317
1318       if(recalcvtx) {
1319         fillthis="hdeclvS_";
1320         fillthis+=ptbin;
1321         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1322
1323         fillthis="hnormdeclvS_";
1324         fillthis+=ptbin;
1325         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1326
1327         fillthis="hd0d0vS_";
1328         fillthis+=ptbin;
1329         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1330       }
1331
1332       if(fSys==0){
1333         fillthis="hcosthetastarS_";
1334         fillthis+=ptbin;
1335         if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1336         else {
1337           if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1338           if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1339         }
1340
1341         fillthis="hcosthpointd0d0S_";
1342         fillthis+=ptbin;
1343         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1344       }
1345
1346     } //end mass cut
1347
1348   } else{ //Background or LS
1349     //if(!fReadMC){
1350     //cout<<"is background"<<endl;
1351
1352     //no mass cut distributions: mass, ptbis
1353     fillthis="hMassB_";
1354     fillthis+=ptbin;
1355     fillthispt="histBkgPt";
1356
1357     if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) {
1358       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1359       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1360     }
1361     if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1362       ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1363       if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1364     }
1365     if(fSys==0){
1366       fillthis="hptB1prongnoMcut_";
1367       fillthis+=ptbin;
1368
1369       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1370
1371       fillthis="hptB2prongsnoMcut_";
1372       fillthis+=ptbin;
1373       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1374       ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1375     }
1376     //apply cut on invariant mass on the pair
1377     if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1378       if(fSys==0){
1379         ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
1380         cosThetaStarD0 = part->CosThetaStarD0();
1381         cosThetaStarD0bar = part->CosThetaStarD0bar();
1382       }
1383
1384       cosPointingAngle = part->CosPointingAngle();
1385       normalizedDecayLength2 = part->NormalizedDecayLength2();
1386       decayLength2 = part->DecayLength2();
1387       decayLengthxy = part->DecayLengthXY();
1388       normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1389       d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1390
1391
1392       AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1393       if(!prongg) {
1394         if(fDebug>2) cout<<"No daughter found";
1395         return;
1396       }
1397       else{
1398         if(fArray==1){
1399           if(prongg->Charge()==1) {
1400             //fTotPosPairs[ptbin]++;
1401             ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
1402           } else {
1403             //fTotNegPairs[ptbin]++;
1404             ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
1405           }
1406         }
1407       }
1408
1409       fillthis="hd0B_";
1410       fillthis+=ptbin;
1411       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1412       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1413
1414       if(fReadMC){
1415         Int_t pdgMother[2]={0,0};
1416         Double_t factor[2]={1,1};
1417
1418         for(Int_t iprong=0;iprong<2;iprong++){
1419           AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1420           lab=prong->GetLabel();
1421           if(lab>=0){
1422             AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1423             if(mcprong){
1424               Int_t labmom=mcprong->GetMother();
1425               if(labmom>=0){
1426                 AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1427                 if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1428               }
1429             }
1430           }
1431
1432           if(fSys==0){
1433
1434             fillthis="hd0moresB_";
1435             fillthis+=ptbin;
1436
1437             if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1438               if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1439               else factor[iprong]=1./.6;
1440               fNentries->Fill(11);
1441             }
1442
1443             if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1444               factor[iprong]=1./0.25;
1445               fNentries->Fill(12);
1446             }
1447             fillthis="hd0moresB_";
1448             fillthis+=ptbin;
1449
1450             ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
1451
1452             if(recalcvtx){
1453               fillthis="hd0vmoresB_";
1454               fillthis+=ptbin;
1455               ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1456             }
1457           }
1458         } //loop on prongs
1459
1460         if(fSys==0){
1461           fillthis="hd0d0moresB_";
1462           fillthis+=ptbin;
1463           ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1464
1465           fillthis="hcosthetapointmoresB_";
1466           fillthis+=ptbin;
1467           ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1468
1469           if(recalcvtx){
1470             fillthis="hd0d0vmoresB_";
1471             fillthis+=ptbin;
1472             ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1473           }
1474         }
1475       } //readMC
1476
1477       if(fSys==0){
1478         //normalise pt distr to half afterwards
1479         fillthis="hptB_";
1480         fillthis+=ptbin;
1481         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1482         ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
1483
1484         fillthis="hcosthetastarB_";
1485         fillthis+=ptbin;
1486         if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1487         if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1488
1489
1490         fillthis="hd0p0B_";
1491         fillthis+=ptbin;
1492         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1493         fillthis="hd0p1B_";
1494         fillthis+=ptbin;
1495         ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1496
1497         fillthis="hcosthpointd0d0B_";
1498         fillthis+=ptbin;
1499         ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1500
1501         fillthis="hcosthpointd0B_";
1502         fillthis+=ptbin;
1503         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1504         ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1505
1506
1507         if(recalcvtx){
1508
1509           fillthis="hd0vp0B_";
1510           fillthis+=ptbin;
1511           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1512           fillthis="hd0vp1B_";
1513           fillthis+=ptbin;
1514           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1515
1516           fillthis="hd0vB_";
1517           fillthis+=ptbin;
1518           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1519           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1520
1521         }
1522
1523       }
1524
1525       fillthis="hdcaB_";
1526       fillthis+=ptbin;
1527       ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1528
1529       fillthis="hd0d0B_";
1530       fillthis+=ptbin;
1531       ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
1532
1533       if(recalcvtx){
1534         fillthis="hd0d0vB_";
1535         fillthis+=ptbin;
1536         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1537       }
1538
1539       fillthis="hcosthetapointB_";
1540       fillthis+=ptbin;
1541       ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1542
1543       fillthis="hcosthetapointxyB_";
1544       fillthis+=ptbin;
1545       ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1546
1547       fillthis="hdeclB_";
1548       fillthis+=ptbin;
1549       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1550
1551       fillthis="hnormdeclB_";
1552       fillthis+=ptbin;
1553       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1554
1555       fillthis="hdeclxyB_";
1556       fillthis+=ptbin;
1557       ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1558
1559       fillthis="hnormdeclxyB_";
1560       fillthis+=ptbin;
1561       ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1562
1563       fillthis="hdeclxyd0d0B_";
1564       fillthis+=ptbin;
1565       ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1566
1567       fillthis="hnormdeclxyd0d0B_";
1568       fillthis+=ptbin;
1569       ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1570
1571       if(recalcvtx) {
1572
1573         fillthis="hdeclvB_";
1574         fillthis+=ptbin;
1575         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1576
1577         fillthis="hnormdeclvB_";
1578         fillthis+=ptbin;
1579         ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1580
1581
1582       }
1583     }//mass cut
1584   }//else (background)
1585
1586   return;
1587 }
1588
1589 //____________________________________________________________________________
1590 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
1591   //
1592   // function used in UserExec to fill mass histograms:
1593   //
1594
1595
1596   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1597
1598   //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
1599
1600   // Fill candidate variable Tree (track selection, no candidate selection)
1601   if( fWriteVariableTree && !part->HasBadDaughters()
1602       && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
1603     fCandidateVariables[0] = part->InvMassD0();
1604     fCandidateVariables[1] = part->InvMassD0bar();
1605     fCandidateVariables[2] = part->Pt();
1606     fCandidateVariables[3] = part->GetDCA();
1607     Double_t ctsD0=0. ,ctsD0bar=0.; part->CosThetaStarD0(ctsD0,ctsD0bar);
1608     fCandidateVariables[4] = ctsD0;
1609     fCandidateVariables[5] = ctsD0bar;
1610     fCandidateVariables[6] = part->Pt2Prong(0);
1611     fCandidateVariables[7] = part->Pt2Prong(1);
1612     fCandidateVariables[8] = part->Getd0Prong(0);
1613     fCandidateVariables[9] = part->Getd0Prong(1);
1614     fCandidateVariables[10] = part->Prodd0d0();
1615     fCandidateVariables[11] = part->CosPointingAngle();
1616     fCandidateVariables[12] = part->CosPointingAngleXY();
1617     fCandidateVariables[13] = part->NormalizedDecayLengthXY();
1618     fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
1619     fVariablesTree->Fill();
1620   }
1621
1622   //cout<<"check cuts = "<<endl;
1623   //cuts->PrintAll();
1624   if (!fIsSelectedCandidate){
1625     //cout<<"Not Selected"<<endl;
1626     //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
1627     return;
1628   }
1629
1630
1631   if(fDebug>2)  cout<<"Candidate selected"<<endl;
1632
1633   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1634   //printf("SELECTED\n");
1635   Int_t ptbin=cuts->PtBin(part->Pt());
1636   Double_t pt = part->Pt();
1637
1638   Double_t impparXY=part->ImpParXY()*10000.;
1639   Double_t trueImpParXY=0.;
1640   Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
1641   Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
1642
1643
1644   // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1645   // if(!prong) {
1646   //   AliDebug(2,"No daughter found");
1647   //   return;
1648   // }
1649   // else{
1650     // if(prong->Charge()==1) {
1651     //   ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
1652     //   //fTotPosPairs[ptbin]++;
1653     // } else {
1654     //   ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
1655     //   //fTotNegPairs[ptbin]++;
1656     // }
1657   //  }
1658
1659   // for(Int_t it=0;it<2;it++){
1660
1661   //    //request on spd points to be addes
1662   //   if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
1663   //     FILE *f=fopen("4display.txt","a");
1664   //     fprintf(f,"pt: %f \n Rapidity: %f \t Period Number: %x \t Run Number: %d \t BunchCrossNumb: %d \t OrbitNumber: %d\n",part->Pt(),part->Y(421),aod->GetPeriodNumber(),aod->GetRunNumber(),aod->GetBunchCrossNumber(),aod->GetOrbitNumber());
1665   //     fclose(f);
1666   //     //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
1667   //   }
1668   // }
1669
1670   TString fillthis="", fillthispt="", fillthismasspt="";
1671   Int_t pdgDgD0toKpi[2]={321,211};
1672   Int_t labD0=-1;
1673   Bool_t isPrimary=kTRUE;
1674   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)
1675
1676   //Define weights for Bayesian (if appropriate)
1677
1678   Double_t weigD0=1.;
1679   Double_t weigD0bar=1.;
1680   if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
1681     weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
1682     weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
1683     if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
1684     if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
1685   }
1686
1687   //count candidates selected by cuts
1688   fNentries->Fill(1);
1689   //count true D0 selected by cuts
1690   if (fReadMC && labD0>=0) fNentries->Fill(2);
1691
1692   if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
1693
1694     arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
1695
1696     if(fReadMC){
1697       if(labD0>=0) {
1698         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1699
1700         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1701         Int_t pdgD0 = partD0->GetPdgCode();
1702         //      cout<<"pdg = "<<pdgD0<<endl;
1703
1704         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
1705         if(!isPrimary)
1706           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
1707         arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
1708
1709         if (pdgD0==421){ //D0
1710           //      cout<<"Fill S with D0"<<endl;
1711           fillthis="histSgn_";
1712           fillthis+=ptbin;
1713           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
1714
1715           if(fFillPtHist){
1716             fillthismasspt="histSgnPt";
1717             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
1718           }
1719           if(fFillImpParHist){
1720             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
1721             else {
1722               fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
1723               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
1724             }
1725           }
1726
1727           if(fSys==0){
1728             if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
1729               fillthis="histSgn27_";
1730               fillthis+=ptbin;
1731               ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
1732             }
1733           }
1734         } else{ //it was a D0bar
1735           fillthis="histRfl_";
1736           fillthis+=ptbin;
1737           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
1738
1739           if(fFillPtHist){
1740             fillthismasspt="histRflPt";
1741             //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
1742             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
1743           }
1744
1745         }
1746       } else {//background
1747         fillthis="histBkg_";
1748         fillthis+=ptbin;
1749         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
1750
1751         if(fFillPtHist){
1752           fillthismasspt="histBkgPt";
1753           //      cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
1754           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
1755         }
1756         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
1757
1758       }
1759
1760     }else{
1761       fillthis="histMass_";
1762       fillthis+=ptbin;
1763       //      cout<<"Filling "<<fillthis<<endl;
1764
1765       //      printf("Fill mass with D0");
1766        ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
1767
1768
1769       if(fFillPtHist){
1770         fillthismasspt="histMassPt";
1771         //      cout<<"Filling "<<fillthismasspt<<endl;
1772         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
1773       }
1774       if(fFillImpParHist) {
1775         //      cout << "Filling fHistMassPtImpParTC[0]"<<endl;
1776         fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
1777       }
1778
1779     }
1780
1781   }
1782   if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
1783
1784     arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
1785
1786     if(fReadMC){
1787       if(labD0>=0) {
1788         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1789         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1790         Int_t pdgD0 = partD0->GetPdgCode();
1791         //      cout<<" pdg = "<<pdgD0<<endl;
1792
1793         if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
1794         if(!isPrimary)
1795           trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
1796         arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
1797
1798         if (pdgD0==-421){ //D0bar
1799           fillthis="histSgn_";
1800           fillthis+=ptbin;
1801           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
1802           // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1803           //   fillthis="histSgn27_";
1804           //   fillthis+=ptbin;
1805           //   ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1806           // }
1807
1808           if(fFillPtHist){
1809             fillthismasspt="histSgnPt";
1810             //      cout<<" Filling "<< fillthismasspt << endl;
1811             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
1812           }
1813           if(fFillImpParHist){
1814             //      cout << " Filling impact parameter thnsparse"<<endl;
1815             if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
1816             else {
1817               fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
1818               fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
1819             }
1820           }
1821
1822         } else{
1823           fillthis="histRfl_";
1824           fillthis+=ptbin;
1825           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
1826           if(fFillPtHist){
1827             fillthismasspt="histRflPt";
1828             //      cout << " Filling "<<fillthismasspt<<endl;
1829             ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
1830           }
1831         }
1832       } else {//background or LS
1833         fillthis="histBkg_";
1834         fillthis+=ptbin;
1835         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
1836
1837         if(fFillPtHist){
1838           fillthismasspt="histBkgPt";
1839           //      cout<<" Filling "<< fillthismasspt << endl;
1840           ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
1841         }
1842         if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
1843
1844       }
1845     }else{
1846       fillthis="histMass_";
1847       fillthis+=ptbin;
1848       //      printf("Fill mass with D0bar");
1849
1850         ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
1851
1852
1853       if(fFillPtHist){
1854         fillthismasspt="histMassPt";
1855         //      cout<<" Filling "<< fillthismasspt << endl;
1856         ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
1857       }
1858       if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
1859
1860     }
1861   }
1862
1863   return;
1864 }
1865
1866 //__________________________________________________________________________
1867 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
1868   //Calculate the primary vertex w/o the daughter tracks of the candidate
1869
1870   Int_t skipped[2];
1871   Int_t nTrksToSkip=2;
1872   AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1873   if(!dgTrack){
1874     AliDebug(2,"no daughter found!");
1875     return 0x0;
1876   }
1877   skipped[0]=dgTrack->GetID();
1878   dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1879   if(!dgTrack){
1880     AliDebug(2,"no daughter found!");
1881     return 0x0;
1882   }
1883   skipped[1]=dgTrack->GetID();
1884
1885   AliESDVertex *vertexESD=0x0;
1886   AliAODVertex *vertexAOD=0x0;
1887   AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1888
1889   //
1890   vertexer->SetSkipTracks(nTrksToSkip,skipped);
1891   vertexer->SetMinClusters(4);
1892   vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
1893   if(!vertexESD) return vertexAOD;
1894   if(vertexESD->GetNContributors()<=0) {
1895     AliDebug(2,"vertexing failed");
1896     delete vertexESD; vertexESD=NULL;
1897     return vertexAOD;
1898   }
1899
1900   delete vertexer; vertexer=NULL;
1901
1902
1903   // convert to AliAODVertex
1904   Double_t pos[3],cov[6],chi2perNDF;
1905   vertexESD->GetXYZ(pos); // position
1906   vertexESD->GetCovMatrix(cov); //covariance matrix
1907   chi2perNDF = vertexESD->GetChi2toNDF();
1908   delete vertexESD; vertexESD=NULL;
1909
1910   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1911   return vertexAOD;
1912
1913 }
1914
1915
1916 //________________________________________________________________________
1917 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1918 {
1919   // Terminate analysis
1920   //
1921   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1922
1923
1924   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1925   if (!fOutputMass) {
1926     printf("ERROR: fOutputMass not available\n");
1927     return;
1928   }
1929   fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
1930   if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
1931     printf("ERROR: fOutputMass not available\n");
1932     return;
1933   }
1934
1935   if(fFillVarHists){
1936     fDistr = dynamic_cast<TList*> (GetOutputData(2));
1937     if (!fDistr) {
1938       printf("ERROR: fDistr not available\n");
1939       return;
1940     }
1941   }
1942
1943   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
1944
1945   if(!fNentries){
1946     printf("ERROR: fNEntries not available\n");
1947     return;
1948   }
1949   fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
1950   if(!fCuts){
1951     printf("ERROR: fCuts not available\n");
1952     return;
1953   }
1954   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
1955   if (!fCounter) {
1956     printf("ERROR: fCounter not available\n");
1957     return;
1958   }
1959   if (fDrawDetSignal) {
1960     fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
1961     if (!fDetSignal) {
1962       printf("ERROR: fDetSignal not available\n");
1963       return;
1964     }
1965   }
1966   Int_t nptbins=fCuts->GetNPtBins();
1967   for(Int_t ipt=0;ipt<nptbins;ipt++){
1968
1969     if(fArray==1 && fFillVarHists){
1970       fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)); //after cuts
1971
1972
1973       if(fLsNormalization>1e-6) {
1974
1975         TString massName="histMass_";
1976         massName+=ipt;
1977         ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1978
1979       }
1980
1981
1982       fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
1983       //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1984
1985       if(fLsNormalization>1e-6) {
1986
1987         TString nameDistr="hptB_";
1988         nameDistr+=ipt;
1989         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1990         nameDistr="hdcaB_";
1991         nameDistr+=ipt;
1992         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1993         nameDistr="hcosthetastarB_";
1994         nameDistr+=ipt;
1995         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1996         nameDistr="hd0B_";
1997         nameDistr+=ipt;
1998         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1999         nameDistr="hd0d0B_";
2000         nameDistr+=ipt;
2001         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2002         nameDistr="hcosthetapointB_";
2003         nameDistr+=ipt;
2004         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2005         if(fSys==0){
2006           nameDistr="hcosthpointd0d0B_";
2007           nameDistr+=ipt;
2008           ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
2009         }
2010       }
2011     }
2012   }
2013   TString cvname,cstname;
2014
2015   if (fArray==0){
2016     cvname="D0invmass";
2017     cstname="cstat0";
2018   } else {
2019     cvname="LSinvmass";
2020     cstname="cstat1";
2021   }
2022
2023   TCanvas *cMass=new TCanvas(cvname,cvname);
2024   cMass->cd();
2025   ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
2026
2027   TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
2028   cStat->cd();
2029   cStat->SetGridy();
2030   fNentries->Draw("htext0");
2031
2032   // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
2033   // ccheck->cd();
2034
2035   return;
2036 }
2037
2038
2039 //________________________________________________________________________
2040 void AliAnalysisTaskSED0Mass::CreateImpactParameterHistos(){
2041   // Histos for impact paramter study
2042
2043   Int_t nmassbins=200;
2044   Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
2045   Int_t fNImpParBins=400;
2046   Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
2047   Int_t nbins[3]={nmassbins,200,fNImpParBins};
2048   Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
2049   Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
2050
2051
2052   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
2053                                         "Mass vs. pt vs.imppar - All",
2054                                         3,nbins,xmin,xmax);
2055   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
2056                                         "Mass vs. pt vs.imppar - promptD",
2057                                         3,nbins,xmin,xmax);
2058   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
2059                                         "Mass vs. pt vs.imppar - DfromB",
2060                                         3,nbins,xmin,xmax);
2061   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
2062                                         "Mass vs. pt vs.true imppar -DfromB",
2063                                         3,nbins,xmin,xmax);
2064   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
2065                                         "Mass vs. pt vs.imppar - backgr.",
2066                                         3,nbins,xmin,xmax);
2067
2068   for(Int_t i=0; i<5;i++){
2069     fOutputMassPt->Add(fHistMassPtImpParTC[i]);
2070   }
2071 }
2072
2073 //_________________________________________________________________________________________________
2074 Float_t AliAnalysisTaskSED0Mass::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
2075   // true impact parameter calculation
2076
2077   printf(" AliAnalysisTaskSED0MassV1::GetTrueImpactParameter() \n");
2078
2079   Double_t vtxTrue[3];
2080   mcHeader->GetVertex(vtxTrue);
2081   Double_t origD[3];
2082   partD0->XvYvZv(origD);
2083   Short_t charge=partD0->Charge();
2084   Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
2085   for(Int_t iDau=0; iDau<2; iDau++){
2086     pXdauTrue[iDau]=0.;
2087     pYdauTrue[iDau]=0.;
2088     pZdauTrue[iDau]=0.;
2089   }
2090
2091   //  Int_t nDau=partD0->GetNDaughters();
2092   Int_t labelFirstDau = partD0->GetDaughter(0);
2093
2094   for(Int_t iDau=0; iDau<2; iDau++){
2095     Int_t ind = labelFirstDau+iDau;
2096     AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
2097     if(!part) continue;
2098     Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
2099     if(!part){
2100       AliError("Daughter particle not found in MC array");
2101       return 99999.;
2102     }
2103     if(pdgCode==211 || pdgCode==321){
2104       pXdauTrue[iDau]=part->Px();
2105       pYdauTrue[iDau]=part->Py();
2106       pZdauTrue[iDau]=part->Pz();
2107     }
2108   }
2109
2110   Double_t d0dummy[2]={0.,0.};
2111   AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
2112   return aodDzeroMC.ImpParXY();
2113
2114 }
2115
2116 //_________________________________________________________________________________________________
2117 Int_t AliAnalysisTaskSED0Mass::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2118   //
2119   // checking whether the mother of the particles come from a charm or a bottom quark
2120   //
2121   printf(" AliAnalysisTaskSED0MassV1::CheckOrigin() \n");
2122
2123   Int_t pdgGranma = 0;
2124   Int_t mother = 0;
2125   mother = mcPartCandidate->GetMother();
2126   Int_t istep = 0;
2127   Int_t abspdgGranma =0;
2128   Bool_t isFromB=kFALSE;
2129   Bool_t isQuarkFound=kFALSE;
2130   while (mother >0 ){
2131     istep++;
2132     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2133     if (mcGranma){
2134       pdgGranma = mcGranma->GetPdgCode();
2135       abspdgGranma = TMath::Abs(pdgGranma);
2136       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2137         isFromB=kTRUE;
2138       }
2139       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
2140       mother = mcGranma->GetMother();
2141     }else{
2142       AliError("Failed casting the mother particle!");
2143       break;
2144     }
2145   }
2146
2147   if(isFromB) return 5;
2148   else return 4;
2149 }