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