]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSED0Mass.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Additional odifications to compile without ZeroMQ
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSED0Mass.cxx
CommitLineData
49061176 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
27de2dfb 16/* $Id$ */
17
49061176 18/////////////////////////////////////////////////////////////
19//
20// AliAnalysisTaskSE for D0 candidates invariant mass histogram
a41f6fad 21// and comparison with the MC truth and cut variables distributions.
49061176 22//
23// Authors: A.Dainese, andrea.dainese@lnl.infn.it
feb73eca 24// Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
25// Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
209003ca 26// Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian
5230fd6a 27////////////////////////////////////////////////////////////
49061176 28
29#include <Riostream.h>
30#include <TClonesArray.h>
527f330b 31#include <TCanvas.h>
49061176 32#include <TNtuple.h>
90c70b48 33#include <TTree.h>
49061176 34#include <TList.h>
35#include <TH1F.h>
a41f6fad 36#include <TH2F.h>
37#include <TDatabasePDG.h>
49061176 38
5b2e5fae 39#include <AliAnalysisDataSlot.h>
40#include <AliAnalysisDataContainer.h>
b557eb43 41#include "AliAnalysisManager.h"
34137226 42#include "AliESDtrack.h"
4e61a020 43#include "AliVertexerTracks.h"
b557eb43 44#include "AliAODHandler.h"
49061176 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"
a96083b9 55#include "AliNormalizationCounter.h"
b557eb43 56
c64cb1f6 57using std::cout;
58using std::endl;
59
49061176 60ClassImp(AliAnalysisTaskSED0Mass)
61
62
63//________________________________________________________________________
64AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
65AliAnalysisTaskSE(),
a220990f 66 fOutputMass(0),
90c70b48 67 fOutputMassPt(0),
13d21bbd 68 fOutputMassY(0),
a220990f 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),
2b35ac47 83 fSys(0),
90c70b48 84 fIsRejectSDDClusters(0),
85 fFillPtHist(kTRUE),
13d21bbd 86 fFillYHist(kFALSE),
90c70b48 87 fFillImpParHist(kFALSE),
0dc0101f 88 fUseSelectionBit(kTRUE),
90c70b48 89 fWriteVariableTree(kFALSE),
90 fVariablesTree(0),
e2aa82b6 91 fCandidateVariables(),
92 fPIDCheck(kFALSE),
93 fDrawDetSignal(kFALSE),
94 fDetSignal(0)
49061176 95{
96 // Default constructor
5da7eaa9 97 for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
90c70b48 98
49061176 99}
100
101//________________________________________________________________________
ea0d8716 102AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
a220990f 103 AliAnalysisTaskSE(name),
90c70b48 104 fOutputMass(0),
105 fOutputMassPt(0),
13d21bbd 106 fOutputMassY(0),
a220990f 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),
2b35ac47 121 fSys(0),
90c70b48 122 fIsRejectSDDClusters(0),
123 fFillPtHist(kTRUE),
13d21bbd 124 fFillYHist(kFALSE),
90c70b48 125 fFillImpParHist(kFALSE),
19f6b9ff 126 fUseSelectionBit(kTRUE),
90c70b48 127 fWriteVariableTree(kFALSE),
128 fVariablesTree(0),
e2aa82b6 129 fCandidateVariables(),
130 fPIDCheck(kFALSE),
131 fDrawDetSignal(kFALSE),
132 fDetSignal(0)
49061176 133{
134 // Default constructor
87020237 135
136 fNPtBins=cuts->GetNPtBins();
4b4a1d25 137
ea0d8716 138 fCuts=cuts;
5da7eaa9 139 for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
ea0d8716 140
141 // Output slot #1 writes into a TList container (mass with cuts)
49061176 142 DefineOutput(1,TList::Class()); //My private output
ea0d8716 143 // Output slot #2 writes into a TList container (distributions)
90c70b48 144 DefineOutput(2,TList::Class()); //My private output
ea0d8716 145 // Output slot #3 writes into a TH1F container (number of events)
a4ae02cd 146 DefineOutput(3,TH1F::Class()); //My private output
700e80e0 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());
90c70b48 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
e2aa82b6 155 // Output slot #8 writes into a TList container (Detector signals)
156 DefineOutput(8, TList::Class()); //My private output
13d21bbd 157 // Output slot #9 stores the mass vs rapidity (y) distributions
158 DefineOutput(9, TList::Class()); //My private output
49061176 159}
160
161//________________________________________________________________________
162AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
163{
a8ce111e 164 if (fOutputMass) {
ea0d8716 165 delete fOutputMass;
166 fOutputMass = 0;
a4ae02cd 167 }
90c70b48 168 if (fOutputMassPt) {
169 delete fOutputMassPt;
170 fOutputMassPt = 0;
171 }
13d21bbd 172 if (fOutputMassY) {
173 delete fOutputMassY;
174 fOutputMassY = 0;
175 }
a41f6fad 176 if (fDistr) {
177 delete fDistr;
178 fDistr = 0;
179 }
ea0d8716 180 if (fCuts) {
181 delete fCuts;
182 fCuts = 0;
49061176 183 }
90c70b48 184 for(Int_t i=0; i<5; i++){
185 if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
186 }
46c96ce5 187 if (fNentries){
188 delete fNentries;
189 fNentries = 0;
190 }
a96083b9 191 if(fCounter){
192 delete fCounter;
193 fCounter=0;
194 }
90c70b48 195 if(fVariablesTree){
196 delete fVariablesTree;
197 fVariablesTree = 0;
198 }
e2aa82b6 199 if (fDetSignal) {
200 delete fDetSignal;
201 fDetSignal = 0;
202 }
ea0d8716 203
49061176 204}
205
206//________________________________________________________________________
207void AliAnalysisTaskSED0Mass::Init()
208{
209 // Initialization
210
211 if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
212
ea0d8716 213
ea0d8716 214 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
700e80e0 215 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
4e61a020 216 copyfCuts->SetName(nameoutput);
ea0d8716 217 // Post the data
700e80e0 218 PostData(4,copyfCuts);
4e61a020 219
49061176 220
49061176 221 return;
222}
223
224//________________________________________________________________________
225void 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
ea0d8716 233 fOutputMass = new TList();
234 fOutputMass->SetOwner();
235 fOutputMass->SetName("listMass");
49061176 236
90c70b48 237 fOutputMassPt = new TList();
238 fOutputMassPt->SetOwner();
239 fOutputMassPt->SetName("listMassPt");
240
13d21bbd 241 fOutputMassY = new TList();
242 fOutputMassY->SetOwner();
243 fOutputMassY->SetName("listMassY");
244
a41f6fad 245 fDistr = new TList();
246 fDistr->SetOwner();
247 fDistr->SetName("distributionslist");
248
e2aa82b6 249 fDetSignal = new TList();
250 fDetSignal->SetOwner();
251 fDetSignal->SetName("detectorsignals");
252
0108fa62 253 TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
90c70b48 254 TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
13d21bbd 255 TString nameMassY="", nameSgnY="", nameBkgY="", nameRflY="";
90c70b48 256 Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
13d21bbd 257 Int_t nbins2dY=60; Float_t binInY=-1.5, binFinY=1.5;
49061176 258
ea0d8716 259 for(Int_t i=0;i<fCuts->GetNPtBins();i++){
b272aebf 260
7646d6da 261 nameMass="histMass_";
ea0d8716 262 nameMass+=i;
9de8c723 263 nameSgn27="histSgn27_";
ea0d8716 264 nameSgn27+=i;
7646d6da 265 nameSgn="histSgn_";
ea0d8716 266 nameSgn+=i;
7646d6da 267 nameBkg="histBkg_";
ea0d8716 268 nameBkg+=i;
a4ae02cd 269 nameRfl="histRfl_";
ea0d8716 270 nameRfl+=i;
0108fa62 271 nameMassNocutsS="hMassS_";
ea0d8716 272 nameMassNocutsS+=i;
0108fa62 273 nameMassNocutsB="hMassB_";
ea0d8716 274 nameMassNocutsB+=i;
7646d6da 275
b272aebf 276 //histograms of cut variable distributions
a8ce111e 277 if(fFillVarHists){
a220990f 278 if(fReadMC){
5230fd6a 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
e8c80c6f 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
a220990f 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);
b272aebf 336
a220990f 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);
b272aebf 348
a220990f 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.);
b272aebf 376
a220990f 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
5230fd6a 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);
e8c80c6f 388 fDistr->Add(hetaphiD0candidateS);
389 fDistr->Add(hetaphiD0candidatesignalregionS);
390 fDistr->Add(hetaphiD0barcandidateS);
391 fDistr->Add(hetaphiD0barcandidatesignalregionS);
a220990f 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 }
b272aebf 417
209003ca 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
209003ca 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());
e8c80c6f 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
a8ce111e 468 // dca
a8ce111e 469 namedistr="hdcaB_";
470 namedistr+=i;
471 TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
b272aebf 472
a8ce111e 473 // impact parameter
a8ce111e 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
a8ce111e 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);
9af24f46 481
a8ce111e 482 //decay lenght
a220990f 483 namedistr="hdeclB_";
9af24f46 484 namedistr+=i;
a220990f 485 TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
9af24f46 486
a220990f 487 namedistr="hnormdeclB_";
9af24f46 488 namedistr+=i;
a220990f 489 TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
a8ce111e 490
a220990f 491 namedistr="hdeclxyB_";
9af24f46 492 namedistr+=i;
a220990f 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.);
9af24f46 497
a220990f 498 namedistr="hdeclxyd0d0B_";
9af24f46 499 namedistr+=i;
a220990f 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);
9af24f46 501
a220990f 502 namedistr="hnormdeclxyd0d0B_";
a8ce111e 503 namedistr+=i;
a220990f 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
a8ce111e 507 namedistr="hcosthetapointB_";
9af24f46 508 namedistr+=i;
a8ce111e 509 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
9af24f46 510
a220990f 511 namedistr="hcosthetapointxyB_";
a8ce111e 512 namedistr+=i;
a220990f 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();
9af24f46 517
a8ce111e 518
209003ca 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);
e8c80c6f 527 fDistr->Add(hetaphiD0candidateB);
528 fDistr->Add(hetaphiD0candidatesignalregionB);
529 fDistr->Add(hetaphiD0barcandidateB);
530 fDistr->Add(hetaphiD0barcandidatesignalregionB);
5230fd6a 531
a8ce111e 532 fDistr->Add(hdcaB);
533
a8ce111e 534 fDistr->Add(hd0B);
535
a8ce111e 536 fDistr->Add(hd0d0B);
537
a8ce111e 538 fDistr->Add(hcosthetapointB);
539
a220990f 540 fDistr->Add(hcosthetapointxyB);
a8ce111e 541
a8ce111e 542 fDistr->Add(hdeclengthB);
543
a8ce111e 544 fDistr->Add(hnormdeclengthB);
545
a8ce111e 546 fDistr->Add(hdeclxyB);
547
a8ce111e 548 fDistr->Add(hnormdeclxyB);
549
a220990f 550 fDistr->Add(hdeclxyd0d0B);
551 fDistr->Add(hnormdeclxyd0d0B);
552
553 fDistr->Add(tmpMB);
554
a8ce111e 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()){
a220990f 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
a8ce111e 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
a8ce111e 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
a8ce111e 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);
9af24f46 611
a8ce111e 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);
9af24f46 615
a8ce111e 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.);
9af24f46 619
a8ce111e 620 fDistr->Add(hd0vB);
621 fDistr->Add(hd0vp0B);
622 fDistr->Add(hd0vp1B);
623 fDistr->Add(hd0vmoresB);
9af24f46 624
a8ce111e 625 fDistr->Add(hd0d0vB);
626 fDistr->Add(hd0d0vmoresB);
9af24f46 627
a8ce111e 628 fDistr->Add(hdeclengthvB);
9af24f46 629
a8ce111e 630 fDistr->Add(hnormdeclengthvB);
631 }
b272aebf 632
a8ce111e 633
634 }
635 //histograms of invariant mass distributions
a4ae02cd 636
0108fa62 637
a220990f 638 //MC signal
639 if(fReadMC){
6756115f 640 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
a4ae02cd 641
a220990f 642 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
643 tmpSt->Sumw2();
644 tmpSl->Sumw2();
a4ae02cd 645
a220990f 646 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
6756115f 647 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
a220990f 648 //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
6756115f 649 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
a220990f 650 //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
651 tmpBt->Sumw2();
652 //tmpBl->Sumw2();
653 tmpRt->Sumw2();
654 //tmpRl->Sumw2();
a4ae02cd 655
a220990f 656 fOutputMass->Add(tmpSt);
657 fOutputMass->Add(tmpRt);
658 fOutputMass->Add(tmpBt);
659
660 }
ea0d8716 661
a8ce111e 662 //mass
6756115f 663 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.5648,2.1648);
a220990f 664 //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
a8ce111e 665 tmpMt->Sumw2();
a220990f 666 //tmpMl->Sumw2();
a8ce111e 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){
a8ce111e 674
a220990f 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
6756115f 719 TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",600,1.5648,2.1648);
a220990f 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_";
a8ce111e 731 namedistr+=i;
a220990f 732 TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
a8ce111e 733
a220990f 734 // costhetastar
735 namedistr="hcosthetastarB_";
a8ce111e 736 namedistr+=i;
a220990f 737 TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
a8ce111e 738
a220990f 739 //pT no mass cut
a8ce111e 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
a220990f 748 fDistr->Add(hptB);
749 fDistr->Add(hcosthetastarB);
750
a8ce111e 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);
0108fa62 767
a8ce111e 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
a8ce111e 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
a8ce111e 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
a220990f 793
a8ce111e 794 fDistr->Add(hcosthpointd0B);
795
a220990f 796
a8ce111e 797 fDistr->Add(hcosthpointd0d0B);
798 }
34137226 799
a8ce111e 800 } //end pp histo
7646d6da 801 }
a4ae02cd 802
a8ce111e 803
804 //for Like sign analysis
805
9af24f46 806 if(fArray==1){
807 namedistr="hpospair";
82487ae7 808 TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
9af24f46 809 namedistr="hnegpair";
82487ae7 810 TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
818c1271 811 fOutputMass->Add(hpospair);
812 fOutputMass->Add(hnegpair);
9af24f46 813 }
449b1302 814
90c70b48 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){
6756115f 826 TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.16484,nbins2dPt,binInPt,binFinPt);
90c70b48 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
6756115f 832 TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
833 TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
90c70b48 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
6756115f 847 TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
90c70b48 848 tmpMtPt->Sumw2();
849
850 fOutputMassPt->Add(tmpMtPt);
851 }
852
853 if(fFillImpParHist) CreateImpactParameterHistos();
13d21bbd 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){
6756115f 869 TH2F* tmpStY = new TH2F(nameSgnY.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.16484,nbins2dY,binInY,binFinY);
13d21bbd 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
6756115f 872 TH2F* tmpRtY = new TH2F(nameRflY.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
873 TH2F* tmpBtY = new TH2F(nameBkgY.Data(), "Background invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
13d21bbd 874 tmpBtY->Sumw2();
875 tmpRtY->Sumw2();
876
877 fOutputMassY->Add(tmpStY);
878 fOutputMassY->Add(tmpRtY);
879 fOutputMassY->Add(tmpBtY);
880 }
6756115f 881 TH2F* tmpMtY = new TH2F(nameMassY.Data(),"D^{0} invariant mass; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
13d21bbd 882 tmpMtY->Sumw2();
883 fOutputMassY->Add(tmpMtY);
884 }
885 }
886
90c70b48 887
5b2e5fae 888 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
4e61a020 889
90c70b48 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);
34137226 891
892 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
4e61a020 893 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
818c1271 894 if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
895 else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
6b3e3c78 896 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
700e80e0 897 fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
898 fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
a8ce111e 899 if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
e2aa82b6 900 if(fFillVarHists || fPIDCheck){
a8ce111e 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 }
700e80e0 910 fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
a3aa1279 911 fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
a220990f 912 if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
2b35ac47 913 if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
90c70b48 914 fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
34137226 915 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
a4ae02cd 916
700e80e0 917 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
38802708 918 fCounter->Init();
a96083b9 919
90c70b48 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];
e2aa82b6 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";
90c70b48 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
e2aa82b6 947 //
948 // Output slot 8 : List for detector response histograms
949 //
e2aa82b6 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
40445ada 963 // Post the data
ea0d8716 964 PostData(1,fOutputMass);
90c70b48 965 PostData(2,fDistr);
40445ada 966 PostData(3,fNentries);
700e80e0 967 PostData(5,fCounter);
90c70b48 968 PostData(6,fOutputMassPt);
969 PostData(7,fVariablesTree);
e2aa82b6 970 PostData(8, fDetSignal);
13d21bbd 971 PostData(9,fOutputMassY);
90c70b48 972
49061176 973 return;
974}
975
976//________________________________________________________________________
977void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
978{
979 // Execute analysis for current event:
980 // heavy flavor candidates association to MC truth
a4ae02cd 981 //cout<<"I'm in UserExec"<<endl;
ea0d8716 982
983
a220990f 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]);
ea0d8716 994
995
49061176 996 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
feb73eca 997
b557eb43 998 TString bname;
feb73eca 999 if(fArray==0){ //D0 candidates
b557eb43 1000 // load D0->Kpi candidates
feb73eca 1001 //cout<<"D0 candidates"<<endl;
b557eb43 1002 bname="D0toKpi";
feb73eca 1003 } else { //LikeSign candidates
feb73eca 1004 // load like sign candidates
b557eb43 1005 //cout<<"LS candidates"<<endl;
1006 bname="LikeSign2Prong";
1007 }
1008
1009 TClonesArray *inputArray=0;
34137226 1010
b557eb43 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());
34137226 1019
b557eb43 1020 if(aodHandler->GetExtensions()) {
1021 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
34137226 1022 AliAODEvent* aodFromExt = ext->GetAOD();
b557eb43 1023 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
feb73eca 1024 }
dc222f77 1025 } else if(aod) {
b557eb43 1026 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
1027 }
feb73eca 1028
b557eb43 1029
dc222f77 1030 if(!inputArray || !aod) {
b557eb43 1031 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
1032 return;
49061176 1033 }
1034
6b3e3c78 1035 // fix for temporary bug in ESDfilter
7c23877d 1036 // the AODs with null vertex pointer didn't pass the PhysSel
5806c290 1037 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
1038
7c23877d 1039
ce39f0ac 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 }
40445ada 1050
ce39f0ac 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 }
49061176 1057 }
1058
b557eb43 1059 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
40445ada 1060
a4ae02cd 1061 //histogram filled with 1 for every AOD
34137226 1062 fNentries->Fill(0);
1879baec 1063 fCounter->StoreEvent(aod,fCuts,fReadMC);
1064 //fCounter->StoreEvent(aod,fReadMC);
a3aa1279 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
d52f7b50 1070 if(!fCuts->IsEventSelected(aod)) {
1071 if(fCuts->GetWhyRejection()==1) // rejected for pileup
700e80e0 1072 fNentries->Fill(13);
a220990f 1073 if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
90c70b48 1074 if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
d52f7b50 1075 return;
1076 }
2b35ac47 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 }
40445ada 1094
1095 // AOD primary vertex
1096 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
b272aebf 1097
40445ada 1098 Bool_t isGoodVtx=kFALSE;
b272aebf 1099
a220990f 1100 //vtx1->Print();
40445ada 1101 TString primTitle = vtx1->GetTitle();
1102 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
1103 isGoodVtx=kTRUE;
1104 fNentries->Fill(3);
1105 }
1106
75638da0 1107 // loop over candidates
1108 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4e61a020 1109 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
75638da0 1110
6b3e3c78 1111 // FILE *f=fopen("4display.txt","a");
1112 // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
a96083b9 1113 Int_t nSelectedloose=0,nSelectedtight=0;
75638da0 1114 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
75638da0 1115 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
818c1271 1116
19f6b9ff 1117 if(fUseSelectionBit && d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
90c70b48 1118 fNentries->Fill(2);
818c1271 1119 continue; //skip the D0 from Dstar
1120 }
a220990f 1121
600faffe 1122 // Bool_t unsetvtx=kFALSE;
1123 // if(!d->GetOwnPrimaryVtx()) {
1124 // d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
1125 // unsetvtx=kTRUE;
1126 // }
ea0d8716 1127
4e61a020 1128
4e61a020 1129 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
a96083b9 1130 nSelectedloose++;
d7688946 1131 nSelectedtight++;
a8ce111e 1132 if(fSys==0){
1133 if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
1134 }
6b3e3c78 1135 Int_t ptbin=fCuts->PtBin(d->Pt());
700e80e0 1136 if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
90c70b48 1137 fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
d7688946 1138 if(fFillVarHists) {
a8ce111e 1139 //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
a220990f 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 //}
d7688946 1150 FillVarHists(aod,d,mcArray,fCuts,fDistr);
1151 }
a8ce111e 1152
e2aa82b6 1153 if (fDrawDetSignal) {
1154 DrawDetSignal(d, fDetSignal);
1155 }
1156
90c70b48 1157 FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
e2aa82b6 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
75638da0 1168 }
ea0d8716 1169
d7688946 1170 fDaughterTracks.Clear();
600faffe 1171 //if(unsetvtx) d->UnsetOwnPrimaryVtx();
75638da0 1172 } //end for prongs
a96083b9 1173 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1174 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
40445ada 1175 // Post the data
ea0d8716 1176 PostData(1,fOutputMass);
90c70b48 1177 PostData(2,fDistr);
40445ada 1178 PostData(3,fNentries);
700e80e0 1179 PostData(5,fCounter);
90c70b48 1180 PostData(6,fOutputMassPt);
1181 PostData(7,fVariablesTree);
e2aa82b6 1182 PostData(8, fDetSignal);
1183
1184 return;
1185}
1186
1187//____________________________________________________________________________
1188void 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;
90c70b48 1221
40445ada 1222 return;
1223}
b272aebf 1224
40445ada 1225//____________________________________________________________________________
4e61a020 1226void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
40445ada 1227 //
1228 // function used in UserExec to fill variable histograms:
1229 //
b272aebf 1230
b272aebf 1231
40445ada 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
a8ce111e 1236 Int_t isSelectedPID=3;
1237 if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
e2aa82b6 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 }
3cc4604b 1245
a8ce111e 1246 if(fCutOnDistr && !fIsSelectedCandidate) return;
d7688946 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;
0108fa62 1261
4e61a020 1262 Double_t invmasscut=0.03;
1263
e8c80c6f 1264 TString fillthispi="",fillthisK="",fillthis="", fillthispt="", fillthisetaphi="";
b272aebf 1265
ea0d8716 1266 Int_t ptbin=cuts->PtBin(part->Pt());
90c70b48 1267 Double_t pt = part->Pt();
6b3e3c78 1268
4e61a020 1269 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1270 dz1[0]=-99; dz2[0]=-99;
9af24f46 1271 Double_t d0[2];
1272 Double_t decl[2]={-99,-99};
1273 Bool_t recalcvtx=kFALSE;
1274
d7688946 1275
1276
9af24f46 1277 if(fCuts->GetIsPrimaryWithoutDaughters()){
1278 recalcvtx=kTRUE;
1279 //recalculate vertex
1280 AliAODVertex *vtxProp=0x0;
d7688946 1281 vtxProp=GetPrimaryVtxSkipped(aod);
9af24f46 1282 if(vtxProp) {
1283 part->SetOwnPrimaryVtx(vtxProp);
1284 //Bool_t unsetvtx=kTRUE;
1285 //Calculate d0 for daughter tracks with Vtx Skipped
d7688946 1286 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1287 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
9af24f46 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
a8ce111e 1298 decl[0]=part->DecayLength2();
1299 decl[1]=part->NormalizedDecayLength2();
9af24f46 1300 part->UnsetOwnPrimaryVtx();
1301
4e61a020 1302 }
1303
d7688946 1304 Double_t cosThetaStarD0 = 99;
1305 Double_t cosThetaStarD0bar = 99;
1306 Double_t cosPointingAngle = 99;
a8ce111e 1307 Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1308 Double_t decayLength2 = -1, decayLengthxy=-1;
d7688946 1309 Double_t ptProng[2]={-99,-99};
a8ce111e 1310 Double_t d0Prong[2]={-99,-99};
e8c80c6f 1311 Double_t etaD = 99.;
1312 Double_t phiD = 99.;
d7688946 1313
0108fa62 1314
a220990f 1315 //disable the PID
1316 if(!fUsePid4Distr) isSelectedPID=0;
1317 if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
9af24f46 1318
a220990f 1319 //check pdg of the prongs
1320 AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1321 AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
6b3e3c78 1322
a220990f 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();
0108fa62 1338 }
a220990f 1339 }
6b3e3c78 1340
a220990f 1341 if(fSys==0){
1342 //no mass cut ditributions: ptbis
40445ada 1343
a220990f 1344 fillthispi="hptpiSnoMcut_";
1345 fillthispi+=ptbin;
a8ce111e 1346
a220990f 1347 fillthisK="hptKSnoMcut_";
1348 fillthisK+=ptbin;
a8ce111e 1349
a220990f 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 }
a8ce111e 1355
a220990f 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());
a8ce111e 1360 }
a220990f 1361 }
a8ce111e 1362
a220990f 1363 //no mass cut ditributions: mass
e8c80c6f 1364
1365 etaD = part->Eta();
1366 phiD = part->Phi();
1367
1368
a220990f 1369 fillthis="hMassS_";
1370 fillthis+=ptbin;
90c70b48 1371 fillthispt="histSgnPt";
40445ada 1372
a220990f 1373 if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1374 || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1375 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
90c70b48 1376 if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
e8c80c6f 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
a220990f 1388 }
1389 else { //D0bar
90c70b48 1390 if(fReadMC || (!fReadMC && isSelectedPID > 1)){
40445ada 1391 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
90c70b48 1392 if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
e8c80c6f 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
90c70b48 1404 }
a220990f 1405 }
4e61a020 1406
a220990f 1407 //apply cut on invariant mass on the pair
1408 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
d7688946 1409
a220990f 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();
a8ce111e 1417
a220990f 1418 ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
a220990f 1419 d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
a8ce111e 1420
a220990f 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();
40445ada 1425
a220990f 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 }
6b3e3c78 1432
a220990f 1433 Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
6b3e3c78 1434
a220990f 1435 if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1436 //cout<<"pi"<<endl;
1437
1438 if(fSys==0){
40445ada 1439 fillthispi="hptpiS_";
1440 fillthispi+=ptbin;
d7688946 1441 ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
a220990f 1442 }
209003ca 1443
a220990f 1444 fillthispi="hd0piS_";
1445 fillthispi+=ptbin;
1446 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1447 if(recalcvtx) {
4e61a020 1448
a220990f 1449 fillthispi="hd0vpiS_";
1450 fillthispi+=ptbin;
1451 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
b272aebf 1452 }
a220990f 1453
1454 }
40445ada 1455
a220990f 1456 if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1457 //cout<<"kappa"<<endl;
1458 if(fSys==0){
40445ada 1459 fillthisK="hptKS_";
1460 fillthisK+=ptbin;
d7688946 1461 ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
a220990f 1462 }
209003ca 1463
209003ca 1464
a220990f 1465 fillthisK="hd0KS_";
1466 fillthisK+=ptbin;
1467 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1468 if (recalcvtx){
1469 fillthisK="hd0vKS_";
40445ada 1470 fillthisK+=ptbin;
a220990f 1471 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
b272aebf 1472 }
a220990f 1473 }
b272aebf 1474
a220990f 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
b272aebf 1481
a220990f 1482 fillthis="hdcaS_";
1483 fillthis+=ptbin;
1484 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
4e61a020 1485
a220990f 1486 fillthis="hcosthetapointS_";
1487 fillthis+=ptbin;
1488 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
4e61a020 1489
a220990f 1490 fillthis="hcosthetapointxyS_";
1491 fillthis+=ptbin;
1492 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
a8ce111e 1493
4e61a020 1494
a220990f 1495 fillthis="hd0d0S_";
1496 fillthis+=ptbin;
1497 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
4e61a020 1498
a220990f 1499 fillthis="hdeclS_";
1500 fillthis+=ptbin;
1501 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
a8ce111e 1502
a220990f 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]);
a8ce111e 1522
a220990f 1523 if(recalcvtx) {
1524 fillthis="hdeclvS_";
a8ce111e 1525 fillthis+=ptbin;
a220990f 1526 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
a2121012 1527
a220990f 1528 fillthis="hnormdeclvS_";
1529 fillthis+=ptbin;
1530 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
a2121012 1531
a220990f 1532 fillthis="hd0d0vS_";
1533 fillthis+=ptbin;
1534 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1535 }
a8ce111e 1536
a220990f 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);
a8ce111e 1544 }
1545
a220990f 1546 fillthis="hcosthpointd0d0S_";
1547 fillthis+=ptbin;
1548 ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1549 }
a8ce111e 1550
5230fd6a 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
4e61a020 1643
a220990f 1644 } else{ //Background or LS
1645 //if(!fReadMC){
1646 //cout<<"is background"<<endl;
e8c80c6f 1647
1648 etaD = part->Eta();
1649 phiD = part->Phi();
1650
a220990f 1651 //no mass cut distributions: mass, ptbis
1652 fillthis="hMassB_";
1653 fillthis+=ptbin;
90c70b48 1654 fillthispt="histBkgPt";
a8ce111e 1655
90c70b48 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);
e8c80c6f 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 }
90c70b48 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);
e8c80c6f 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
90c70b48 1684 }
a220990f 1685 if(fSys==0){
1686 fillthis="hptB1prongnoMcut_";
1687 fillthis+=ptbin;
40445ada 1688
a220990f 1689 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
40445ada 1690
a220990f 1691 fillthis="hptB2prongsnoMcut_";
1692 fillthis+=ptbin;
1693 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
818c1271 1694 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
a220990f 1695 }
209003ca 1696
209003ca 1697
e8c80c6f 1698 //apply cut on invariant mass on the pair
1699 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
a220990f 1700 if(fSys==0){
818c1271 1701 ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
d7688946 1702 cosThetaStarD0 = part->CosThetaStarD0();
1703 cosThetaStarD0bar = part->CosThetaStarD0bar();
a220990f 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);
e8c80c6f 1712
a220990f 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]++;
600faffe 1723 ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
a220990f 1724 } else {
1725 //fTotNegPairs[ptbin]++;
600faffe 1726 ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
449b1302 1727 }
0108fa62 1728 }
a220990f 1729 }
209003ca 1730
1731 //fill pt and phi distrib for prongs with M cut
1732
1733 if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))){
209003ca 1734 for(Int_t it=0; it<2; it++){
5230fd6a 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());
e8c80c6f 1741
209003ca 1742 Int_t nPointsITS = 0;
1743 for (Int_t il=0; il<6; il++){
1744 if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1745 }
5230fd6a 1746 fillthis="hNITSpointsD0vsptB_";
1747 fillthis+=ptbin;
1748 ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
209003ca 1749 fillthis="hNSPDpointsD0B_";
1750 fillthis+=ptbin;
1751 if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
5230fd6a 1752 ((TH1I*)listout->FindObject(fillthis))->Fill(0);
209003ca 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 }
5230fd6a 1763 fillthis="hNclsD0vsptB_";
209003ca 1764 fillthis+=ptbin;
5230fd6a 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);
209003ca 1768 }
1769
1770
1771 }
1772
1773 if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
5230fd6a 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 }
209003ca 1787 }
40445ada 1788
a220990f 1789 fillthis="hd0B_";
1790 fillthis+=ptbin;
1791 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1792 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
0108fa62 1793
a220990f 1794 if(fReadMC){
1795 Int_t pdgMother[2]={0,0};
1796 Double_t factor[2]={1,1};
4e61a020 1797
a220990f 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();
6b3e3c78 1808 }
1809 }
a220990f 1810 }
6b3e3c78 1811
a220990f 1812 if(fSys==0){
a8ce111e 1813
a220990f 1814 fillthis="hd0moresB_";
1815 fillthis+=ptbin;
6b3e3c78 1816
a220990f 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 }
6b3e3c78 1822
a220990f 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;
6b3e3c78 1829
a220990f 1830 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
6b3e3c78 1831
a220990f 1832 if(recalcvtx){
1833 fillthis="hd0vmoresB_";
1834 fillthis+=ptbin;
1835 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
9af24f46 1836 }
a220990f 1837 }
1838 } //loop on prongs
6b3e3c78 1839
a220990f 1840 if(fSys==0){
6b3e3c78 1841 fillthis="hd0d0moresB_";
1842 fillthis+=ptbin;
1843 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1844
a8ce111e 1845 fillthis="hcosthetapointmoresB_";
1846 fillthis+=ptbin;
1847 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1848
9af24f46 1849 if(recalcvtx){
1850 fillthis="hd0d0vmoresB_";
1851 fillthis+=ptbin;
1852 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1853 }
a220990f 1854 }
1855 } //readMC
ce39f0ac 1856
a220990f 1857 if(fSys==0){
1858 //normalise pt distr to half afterwards
1859 fillthis="hptB_";
40445ada 1860 fillthis+=ptbin;
a220990f 1861 ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1862 ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
34137226 1863
40445ada 1864 fillthis="hcosthetastarB_";
1865 fillthis+=ptbin;
d7688946 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);
40445ada 1868
a220990f 1869
1870 fillthis="hd0p0B_";
40445ada 1871 fillthis+=ptbin;
a220990f 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
40445ada 1886
9af24f46 1887 if(recalcvtx){
a220990f 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_";
a2121012 1897 fillthis+=ptbin;
a220990f 1898 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1899 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1900
a2121012 1901 }
1902
a220990f 1903 }
40445ada 1904
a220990f 1905 fillthis="hdcaB_";
1906 fillthis+=ptbin;
1907 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
a2121012 1908
a220990f 1909 fillthis="hd0d0B_";
1910 fillthis+=ptbin;
1911 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
a8ce111e 1912
a220990f 1913 if(recalcvtx){
1914 fillthis="hd0d0vB_";
a8ce111e 1915 fillthis+=ptbin;
a220990f 1916 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1917 }
a8ce111e 1918
a220990f 1919 fillthis="hcosthetapointB_";
1920 fillthis+=ptbin;
1921 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
a2121012 1922
a220990f 1923 fillthis="hcosthetapointxyB_";
1924 fillthis+=ptbin;
1925 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
9af24f46 1926
a220990f 1927 fillthis="hdeclB_";
1928 fillthis+=ptbin;
1929 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
a2121012 1930
a220990f 1931 fillthis="hnormdeclB_";
1932 fillthis+=ptbin;
1933 ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
9af24f46 1934
a220990f 1935 fillthis="hdeclxyB_";
1936 fillthis+=ptbin;
1937 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
a8ce111e 1938
a220990f 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;
e8c80c6f 1949 ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1950
a220990f 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)
d7688946 1966
49061176 1967 return;
1968}
ea0d8716 1969
90c70b48 1970//____________________________________________________________________________
1971void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
49061176 1972 //
40445ada 1973 // function used in UserExec to fill mass histograms:
49061176 1974 //
9de8c723 1975
1976
a8ce111e 1977 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
9de8c723 1978
d7688946 1979 //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
feb73eca 1980
90c70b48 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
ea0d8716 2003 //cout<<"check cuts = "<<endl;
2004 //cuts->PrintAll();
d7688946 2005 if (!fIsSelectedCandidate){
ea0d8716 2006 //cout<<"Not Selected"<<endl;
6b3e3c78 2007 //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
ea0d8716 2008 return;
2009 }
9de8c723 2010
6b3e3c78 2011
4e61a020 2012 if(fDebug>2) cout<<"Candidate selected"<<endl;
a41f6fad 2013
ea0d8716 2014 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2015 //printf("SELECTED\n");
2016 Int_t ptbin=cuts->PtBin(part->Pt());
90c70b48 2017 Double_t pt = part->Pt();
13d21bbd 2018 Double_t y = part->YD0();
90c70b48 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
9de8c723 2025
818c1271 2026 // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2027 // if(!prong) {
2028 // AliDebug(2,"No daughter found");
2029 // return;
2030 // }
2031 // else{
a8ce111e 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 // }
818c1271 2039 // }
700e80e0 2040
9af24f46 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 // }
6b3e3c78 2051
13d21bbd 2052 TString fillthis="", fillthispt="", fillthismasspt="", fillthismassy="";
ea0d8716 2053 Int_t pdgDgD0toKpi[2]={321,211};
2054 Int_t labD0=-1;
90c70b48 2055 Bool_t isPrimary=kTRUE;
ea0d8716 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
e2aa82b6 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];
19678166 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
e2aa82b6 2067 }
2068
ea0d8716 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);
ea0d8716 2073
d7688946 2074 if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
90c70b48 2075
2076 arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
2077
a220990f 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();
2b35ac47 2084 // cout<<"pdg = "<<pdgD0<<endl;
90c70b48 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
a220990f 2091 if (pdgD0==421){ //D0
2b35ac47 2092 // cout<<"Fill S with D0"<<endl;
a220990f 2093 fillthis="histSgn_";
2094 fillthis+=ptbin;
e2aa82b6 2095 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
90c70b48 2096
2097 if(fFillPtHist){
2098 fillthismasspt="histSgnPt";
e2aa82b6 2099 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
90c70b48 2100 }
2101 if(fFillImpParHist){
e2aa82b6 2102 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
90c70b48 2103 else {
e2aa82b6 2104 fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
2105 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
90c70b48 2106 }
2107 }
2108
13d21bbd 2109 if(fFillYHist){
2110 fillthismassy="histSgnY_";
2111 fillthismassy+=ptbin;
2112 ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2113 }
2114
a220990f 2115 if(fSys==0){
2b35ac47 2116 if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
a220990f 2117 fillthis="histSgn27_";
2118 fillthis+=ptbin;
e2aa82b6 2119 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
a220990f 2120 }
2121 }
2122 } else{ //it was a D0bar
2123 fillthis="histRfl_";
a8ce111e 2124 fillthis+=ptbin;
e2aa82b6 2125 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
90c70b48 2126
2127 if(fFillPtHist){
2128 fillthismasspt="histRflPt";
2129 // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
e2aa82b6 2130 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
90c70b48 2131 }
2132
13d21bbd 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
a8ce111e 2140 }
a220990f 2141 } else {//background
2142 fillthis="histBkg_";
a4ae02cd 2143 fillthis+=ptbin;
e2aa82b6 2144 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
90c70b48 2145
2146 if(fFillPtHist){
2147 fillthismasspt="histBkgPt";
2148 // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
e2aa82b6 2149 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
90c70b48 2150 }
e2aa82b6 2151 if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
90c70b48 2152
13d21bbd 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
a4ae02cd 2160 }
a220990f 2161
2162 }else{
2163 fillthis="histMass_";
a4ae02cd 2164 fillthis+=ptbin;
2b35ac47 2165 // cout<<"Filling "<<fillthis<<endl;
a220990f 2166
2b35ac47 2167 // printf("Fill mass with D0");
e2aa82b6 2168 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2169
90c70b48 2170
2171 if(fFillPtHist){
2172 fillthismasspt="histMassPt";
2173 // cout<<"Filling "<<fillthismasspt<<endl;
e2aa82b6 2174 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
90c70b48 2175 }
2176 if(fFillImpParHist) {
2177 // cout << "Filling fHistMassPtImpParTC[0]"<<endl;
e2aa82b6 2178 fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
90c70b48 2179 }
2180
13d21bbd 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
ea0d8716 2188 }
a220990f 2189
ea0d8716 2190 }
d7688946 2191 if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
90c70b48 2192
2193 arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
2194
a220990f 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();
2b35ac47 2200 // cout<<" pdg = "<<pdgD0<<endl;
90c70b48 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
a220990f 2207 if (pdgD0==-421){ //D0bar
2208 fillthis="histSgn_";
2209 fillthis+=ptbin;
e2aa82b6 2210 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
a220990f 2211 // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
2212 // fillthis="histSgn27_";
2213 // fillthis+=ptbin;
2214 // ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
2215 // }
9de8c723 2216
90c70b48 2217 if(fFillPtHist){
2218 fillthismasspt="histSgnPt";
2219 // cout<<" Filling "<< fillthismasspt << endl;
e2aa82b6 2220 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
90c70b48 2221 }
2222 if(fFillImpParHist){
2223 // cout << " Filling impact parameter thnsparse"<<endl;
e2aa82b6 2224 if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
90c70b48 2225 else {
e2aa82b6 2226 fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
2227 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
90c70b48 2228 }
2229 }
13d21bbd 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 }
a4ae02cd 2237
a220990f 2238 } else{
2239 fillthis="histRfl_";
2240 fillthis+=ptbin;
e2aa82b6 2241 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
90c70b48 2242 if(fFillPtHist){
2243 fillthismasspt="histRflPt";
2244 // cout << " Filling "<<fillthismasspt<<endl;
e2aa82b6 2245 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
90c70b48 2246 }
13d21bbd 2247 if(fFillYHist){
2248 fillthismassy="histRflY_";
2249 fillthismassy+=ptbin;
2250 // cout << " Filling "<<fillthismassy<<endl;
2251 ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2252 }
a220990f 2253 }
2254 } else {//background or LS
2255 fillthis="histBkg_";
a4ae02cd 2256 fillthis+=ptbin;
e2aa82b6 2257 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
90c70b48 2258
2259 if(fFillPtHist){
2260 fillthismasspt="histBkgPt";
2261 // cout<<" Filling "<< fillthismasspt << endl;
e2aa82b6 2262 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
90c70b48 2263 }
e2aa82b6 2264 if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
13d21bbd 2265 if(fFillYHist){
2266 fillthismassy="histBkgY_";
2267 fillthismassy+=ptbin;
2268 // cout<<" Filling "<< fillthismassy << endl;
2269 ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2270 }
a4ae02cd 2271 }
a220990f 2272 }else{
2273 fillthis="histMass_";
ea0d8716 2274 fillthis+=ptbin;
2b35ac47 2275 // printf("Fill mass with D0bar");
e2aa82b6 2276
2277 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
2278
a220990f 2279
90c70b48 2280 if(fFillPtHist){
2281 fillthismasspt="histMassPt";
2282 // cout<<" Filling "<< fillthismasspt << endl;
e2aa82b6 2283 ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
90c70b48 2284 }
e2aa82b6 2285 if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
13d21bbd 2286 if(fFillYHist){
2287 fillthismassy="histMassY_";
2288 fillthismassy+=ptbin;
2289 // cout<<" Filling "<< fillthismassy << endl;
2290 ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2291 }
a4ae02cd 2292 }
ea0d8716 2293 }
a4ae02cd 2294
40445ada 2295 return;
49061176 2296}
4e61a020 2297
2298//__________________________________________________________________________
d7688946 2299AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
4e61a020 2300 //Calculate the primary vertex w/o the daughter tracks of the candidate
2301
4e61a020 2302 Int_t skipped[2];
2303 Int_t nTrksToSkip=2;
d7688946 2304 AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
4e61a020 2305 if(!dgTrack){
2306 AliDebug(2,"no daughter found!");
2307 return 0x0;
2308 }
2309 skipped[0]=dgTrack->GetID();
d7688946 2310 dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
4e61a020 2311 if(!dgTrack){
2312 AliDebug(2,"no daughter found!");
2313 return 0x0;
2314 }
2315 skipped[1]=dgTrack->GetID();
2316
6b3e3c78 2317 AliESDVertex *vertexESD=0x0;
2318 AliAODVertex *vertexAOD=0x0;
2319 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
2320
4e61a020 2321 //
2322 vertexer->SetSkipTracks(nTrksToSkip,skipped);
4e61a020 2323 vertexer->SetMinClusters(4);
a2121012 2324 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
4e61a020 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
49061176 2348//________________________________________________________________________
2349void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
2350{
2351 // Terminate analysis
2352 //
2353 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
2354
6321ee46 2355
ea0d8716 2356 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
2357 if (!fOutputMass) {
2358 printf("ERROR: fOutputMass not available\n");
a4ae02cd 2359 return;
2360 }
90c70b48 2361 fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
2362 if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
2363 printf("ERROR: fOutputMass not available\n");
2364 return;
2365 }
2366
a8ce111e 2367 if(fFillVarHists){
2368 fDistr = dynamic_cast<TList*> (GetOutputData(2));
2369 if (!fDistr) {
2370 printf("ERROR: fDistr not available\n");
2371 return;
2372 }
a41f6fad 2373 }
a8ce111e 2374
40445ada 2375 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
5b2e5fae 2376
40445ada 2377 if(!fNentries){
2378 printf("ERROR: fNEntries not available\n");
2379 return;
2380 }
700e80e0 2381 fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
2382 if(!fCuts){
2383 printf("ERROR: fCuts not available\n");
34137226 2384 return;
2385 }
700e80e0 2386 fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
a96083b9 2387 if (!fCounter) {
2388 printf("ERROR: fCounter not available\n");
2389 return;
2390 }
e2aa82b6 2391 if (fDrawDetSignal) {
2392 fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
2393 if (!fDetSignal) {
2394 printf("ERROR: fDetSignal not available\n");
2395 return;
2396 }
2397 }
13d21bbd 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
700e80e0 2406 Int_t nptbins=fCuts->GetNPtBins();
2407 for(Int_t ipt=0;ipt<nptbins;ipt++){
4e61a020 2408
a8ce111e 2409 if(fArray==1 && fFillVarHists){
84f75f2e 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
feb73eca 2411
2412
ea0d8716 2413 if(fLsNormalization>1e-6) {
9de8c723 2414
feb73eca 2415 TString massName="histMass_";
ea0d8716 2416 massName+=ipt;
2417 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
2418
feb73eca 2419 }
40445ada 2420
feb73eca 2421
84f75f2e 2422 fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
700e80e0 2423 //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
40445ada 2424
ea0d8716 2425 if(fLsNormalization>1e-6) {
40445ada 2426
57c5a6f5 2427 TString nameDistr="hdcaB_";
ea0d8716 2428 nameDistr+=ipt;
40445ada 2429 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2430 nameDistr="hd0B_";
ea0d8716 2431 nameDistr+=ipt;
40445ada 2432 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2433 nameDistr="hd0d0B_";
ea0d8716 2434 nameDistr+=ipt;
40445ada 2435 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2436 nameDistr="hcosthetapointB_";
ea0d8716 2437 nameDistr+=ipt;
40445ada 2438 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
a220990f 2439 if(fSys==0){
57c5a6f5 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());
a220990f 2446 nameDistr="hcosthpointd0d0B_";
2447 nameDistr+=ipt;
2448 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
2449 }
40445ada 2450 }
feb73eca 2451 }
2452 }
4e61a020 2453 TString cvname,cstname;
527f330b 2454
2455 if (fArray==0){
2456 cvname="D0invmass";
4e61a020 2457 cstname="cstat0";
2458 } else {
2459 cvname="LSinvmass";
2460 cstname="cstat1";
2461 }
527f330b 2462
34137226 2463 TCanvas *cMass=new TCanvas(cvname,cvname);
2464 cMass->cd();
ea0d8716 2465 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
527f330b 2466
4e61a020 2467 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
34137226 2468 cStat->cd();
2469 cStat->SetGridy();
2470 fNentries->Draw("htext0");
527f330b 2471
4e61a020 2472 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
2473 // ccheck->cd();
2474
49061176 2475 return;
2476}
2477
90c70b48 2478
2479//________________________________________________________________________
2480void 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//_________________________________________________________________________________________________
2514Float_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//_________________________________________________________________________________________________
2557Int_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}