PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx

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