PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx

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