]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSED0Mass.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for centrality selection usage (Davide)
[u/mrichter/AliRoot.git] / PWG3 / 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)
49061176 26/////////////////////////////////////////////////////////////
27
28#include <Riostream.h>
29#include <TClonesArray.h>
527f330b 30#include <TCanvas.h>
49061176 31#include <TNtuple.h>
32#include <TList.h>
33#include <TH1F.h>
a41f6fad 34#include <TH2F.h>
35#include <TDatabasePDG.h>
49061176 36
5b2e5fae 37#include <AliAnalysisDataSlot.h>
38#include <AliAnalysisDataContainer.h>
b557eb43 39#include "AliAnalysisManager.h"
34137226 40#include "AliESDtrack.h"
4e61a020 41#include "AliVertexerTracks.h"
b557eb43 42#include "AliAODHandler.h"
49061176 43#include "AliAODEvent.h"
44#include "AliAODVertex.h"
45#include "AliAODTrack.h"
46#include "AliAODMCHeader.h"
47#include "AliAODMCParticle.h"
48#include "AliAODRecoDecayHF2Prong.h"
49#include "AliAODRecoCascadeHF.h"
50#include "AliAnalysisVertexingHF.h"
51#include "AliAnalysisTaskSE.h"
52#include "AliAnalysisTaskSED0Mass.h"
a96083b9 53#include "AliNormalizationCounter.h"
b557eb43 54
49061176 55ClassImp(AliAnalysisTaskSED0Mass)
56
57
58//________________________________________________________________________
59AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
60AliAnalysisTaskSE(),
a220990f 61 fOutputMass(0),
62 fDistr(0),
63 fNentries(0),
64 fCuts(0),
65 fArray(0),
66 fReadMC(0),
67 fCutOnDistr(0),
68 fUsePid4Distr(0),
69 fCounter(0),
70 fNPtBins(1),
71 fLsNormalization(1.),
72 fFillOnlyD0D0bar(0),
73 fDaughterTracks(),
74 fIsSelectedCandidate(0),
75 fFillVarHists(kTRUE),
76 fSys(0)
49061176 77{
78 // Default constructor
79}
80
81//________________________________________________________________________
ea0d8716 82AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
a220990f 83 AliAnalysisTaskSE(name),
84 fOutputMass(0),
85 fDistr(0),
86 fNentries(0),
87 fCuts(0),
88 fArray(0),
89 fReadMC(0),
90 fCutOnDistr(0),
91 fUsePid4Distr(0),
92 fCounter(0),
93 fNPtBins(1),
94 fLsNormalization(1.),
95 fFillOnlyD0D0bar(0),
96 fDaughterTracks(),
97 fIsSelectedCandidate(0),
98 fFillVarHists(kTRUE),
99 fSys(0)
49061176 100{
101 // Default constructor
87020237 102
103 fNPtBins=cuts->GetNPtBins();
4b4a1d25 104
ea0d8716 105 fCuts=cuts;
106
107 // Output slot #1 writes into a TList container (mass with cuts)
49061176 108 DefineOutput(1,TList::Class()); //My private output
ea0d8716 109 // Output slot #2 writes into a TList container (distributions)
a8ce111e 110 if(fFillVarHists) DefineOutput(2,TList::Class()); //My private output
ea0d8716 111 // Output slot #3 writes into a TH1F container (number of events)
a4ae02cd 112 DefineOutput(3,TH1F::Class()); //My private output
700e80e0 113 // Output slot #4 writes into a TList container (cuts)
114 DefineOutput(4,AliRDHFCutsD0toKpi::Class()); //My private output
115 // Output slot #5 writes Normalization Counter
116 DefineOutput(5,AliNormalizationCounter::Class());
49061176 117}
118
119//________________________________________________________________________
120AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
121{
a8ce111e 122 if (fOutputMass) {
ea0d8716 123 delete fOutputMass;
124 fOutputMass = 0;
a4ae02cd 125 }
a41f6fad 126 if (fDistr) {
127 delete fDistr;
128 fDistr = 0;
129 }
ea0d8716 130 if (fCuts) {
131 delete fCuts;
132 fCuts = 0;
49061176 133 }
46c96ce5 134 if (fNentries){
135 delete fNentries;
136 fNentries = 0;
137 }
a96083b9 138 if(fCounter){
139 delete fCounter;
140 fCounter=0;
141 }
ea0d8716 142
49061176 143}
144
145//________________________________________________________________________
146void AliAnalysisTaskSED0Mass::Init()
147{
148 // Initialization
149
150 if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
151
ea0d8716 152
ea0d8716 153 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
700e80e0 154 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
4e61a020 155 copyfCuts->SetName(nameoutput);
ea0d8716 156 // Post the data
700e80e0 157 PostData(4,copyfCuts);
4e61a020 158
49061176 159
49061176 160 return;
161}
162
163//________________________________________________________________________
164void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
165{
166
167 // Create the output container
168 //
169 if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
170
171 // Several histograms are more conveniently managed in a TList
ea0d8716 172 fOutputMass = new TList();
173 fOutputMass->SetOwner();
174 fOutputMass->SetName("listMass");
49061176 175
a41f6fad 176 fDistr = new TList();
177 fDistr->SetOwner();
178 fDistr->SetName("distributionslist");
179
0108fa62 180 TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
49061176 181
ea0d8716 182 for(Int_t i=0;i<fCuts->GetNPtBins();i++){
b272aebf 183
7646d6da 184 nameMass="histMass_";
ea0d8716 185 nameMass+=i;
9de8c723 186 nameSgn27="histSgn27_";
ea0d8716 187 nameSgn27+=i;
7646d6da 188 nameSgn="histSgn_";
ea0d8716 189 nameSgn+=i;
7646d6da 190 nameBkg="histBkg_";
ea0d8716 191 nameBkg+=i;
a4ae02cd 192 nameRfl="histRfl_";
ea0d8716 193 nameRfl+=i;
0108fa62 194 nameMassNocutsS="hMassS_";
ea0d8716 195 nameMassNocutsS+=i;
0108fa62 196 nameMassNocutsB="hMassB_";
ea0d8716 197 nameMassNocutsB+=i;
7646d6da 198
b272aebf 199 //histograms of cut variable distributions
a8ce111e 200 if(fFillVarHists){
a220990f 201 if(fReadMC){
202 // dca
203 namedistr="hdcaS_";
204 namedistr+=i;
205 TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
206 // impact parameter
207 namedistr="hd0piS_";
208 namedistr+=i;
209 TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
b272aebf 210
a220990f 211 namedistr="hd0KS_";
212 namedistr+=i;
213 TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
214 namedistr="hd0d0S_";
215 namedistr+=i;
216 TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
217
218 //decay lenght
219 namedistr="hdeclS_";
220 namedistr+=i;
221 TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
b272aebf 222
a220990f 223 namedistr="hnormdeclS_";
224 namedistr+=i;
225 TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
226
227 namedistr="hdeclxyS_";
228 namedistr+=i;
229 TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
230 namedistr="hnormdeclxyS_";
231 namedistr+=i;
232 TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
233
234 namedistr="hdeclxyd0d0S_";
235 namedistr+=i;
236 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);
237
238 namedistr="hnormdeclxyd0d0S_";
239 namedistr+=i;
240 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);
241
242 // costhetapoint
243 namedistr="hcosthetapointS_";
244 namedistr+=i;
245 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
246
247 namedistr="hcosthetapointxyS_";
248 namedistr+=i;
249 TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
b272aebf 250
a220990f 251 TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
252 tmpMS->Sumw2();
253
254
255 fDistr->Add(hdcaS);
256
257 fDistr->Add(hd0piS);
258 fDistr->Add(hd0KS);
259
260 fDistr->Add(hd0d0S);
261
262 fDistr->Add(hcosthetapointS);
263
264 fDistr->Add(hcosthetapointxyS);
265
266 fDistr->Add(hdeclengthS);
267
268 fDistr->Add(hnormdeclengthS);
269
270 fDistr->Add(hdeclxyS);
271
272 fDistr->Add(hnormdeclxyS);
273
274 fDistr->Add(hdeclxyd0d0S);
275 fDistr->Add(hnormdeclxyd0d0S);
276
277 fDistr->Add(tmpMS);
278 }
b272aebf 279
a8ce111e 280 // dca
a8ce111e 281 namedistr="hdcaB_";
282 namedistr+=i;
283 TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
b272aebf 284
a8ce111e 285 // impact parameter
a8ce111e 286 namedistr="hd0B_";
287 namedistr+=i;
288 TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
289
a8ce111e 290 namedistr="hd0d0B_";
291 namedistr+=i;
292 TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
9af24f46 293
a8ce111e 294 //decay lenght
a220990f 295 namedistr="hdeclB_";
9af24f46 296 namedistr+=i;
a220990f 297 TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
9af24f46 298
a220990f 299 namedistr="hnormdeclB_";
9af24f46 300 namedistr+=i;
a220990f 301 TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
a8ce111e 302
a220990f 303 namedistr="hdeclxyB_";
9af24f46 304 namedistr+=i;
a220990f 305 TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
306 namedistr="hnormdeclxyB_";
307 namedistr+=i;
308 TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
9af24f46 309
a220990f 310 namedistr="hdeclxyd0d0B_";
9af24f46 311 namedistr+=i;
a220990f 312 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 313
a220990f 314 namedistr="hnormdeclxyd0d0B_";
a8ce111e 315 namedistr+=i;
a220990f 316 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);
317
318 // costhetapoint
a8ce111e 319 namedistr="hcosthetapointB_";
9af24f46 320 namedistr+=i;
a8ce111e 321 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
9af24f46 322
a220990f 323 namedistr="hcosthetapointxyB_";
a8ce111e 324 namedistr+=i;
a220990f 325 TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
326
327 TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
328 tmpMB->Sumw2();
9af24f46 329
a8ce111e 330
a8ce111e 331
a8ce111e 332 fDistr->Add(hdcaB);
333
a8ce111e 334 fDistr->Add(hd0B);
335
a8ce111e 336 fDistr->Add(hd0d0B);
337
a8ce111e 338 fDistr->Add(hcosthetapointB);
339
a220990f 340 fDistr->Add(hcosthetapointxyB);
a8ce111e 341
a8ce111e 342 fDistr->Add(hdeclengthB);
343
a8ce111e 344 fDistr->Add(hnormdeclengthB);
345
a8ce111e 346 fDistr->Add(hdeclxyB);
347
a8ce111e 348 fDistr->Add(hnormdeclxyB);
349
a220990f 350 fDistr->Add(hdeclxyd0d0B);
351 fDistr->Add(hnormdeclxyd0d0B);
352
353 fDistr->Add(tmpMB);
354
a8ce111e 355 //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
356
357 if(fCuts->GetIsPrimaryWithoutDaughters()){
a220990f 358 if(fReadMC){
359 namedistr="hd0vpiS_";
360 namedistr+=i;
361 TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
362
363 namedistr="hd0vKS_";
364 namedistr+=i;
365 TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
366
367 namedistr="hd0d0vS_";
368 namedistr+=i;
369 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);
370
371 namedistr="hdeclvS_";
372 namedistr+=i;
373 TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
374
375 namedistr="hnormdeclvS_";
376 namedistr+=i;
377 TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
378
379 fDistr->Add(hd0vpiS);
380 fDistr->Add(hd0vKS);
381 fDistr->Add(hd0d0vS);
382 fDistr->Add(hdeclengthvS);
383 fDistr->Add(hnormdeclengthvS);
384
385 }
386
a8ce111e 387 namedistr="hd0vmoresB_";
388 namedistr+=i;
389 TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
390
391 namedistr="hd0d0vmoresB_";
392 namedistr+=i;
393 TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
394
395
a8ce111e 396 namedistr="hd0vB_";
397 namedistr+=i;
398 TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
399
400 namedistr="hd0vp0B_";
401 namedistr+=i;
402 TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
403
404 namedistr="hd0vp1B_";
405 namedistr+=i;
406 TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
407
a8ce111e 408 namedistr="hd0d0vB_";
409 namedistr+=i;
410 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 411
a8ce111e 412 namedistr="hdeclvB_";
413 namedistr+=i;
414 TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
9af24f46 415
a8ce111e 416 namedistr="hnormdeclvB_";
417 namedistr+=i;
418 TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
9af24f46 419
a8ce111e 420 fDistr->Add(hd0vB);
421 fDistr->Add(hd0vp0B);
422 fDistr->Add(hd0vp1B);
423 fDistr->Add(hd0vmoresB);
9af24f46 424
a8ce111e 425 fDistr->Add(hd0d0vB);
426 fDistr->Add(hd0d0vmoresB);
9af24f46 427
a8ce111e 428 fDistr->Add(hdeclengthvB);
9af24f46 429
a8ce111e 430 fDistr->Add(hnormdeclengthvB);
431 }
b272aebf 432
a8ce111e 433
434 }
435 //histograms of invariant mass distributions
a4ae02cd 436
0108fa62 437
a220990f 438 //MC signal
439 if(fReadMC){
440 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
a4ae02cd 441
a220990f 442 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
443 tmpSt->Sumw2();
444 tmpSl->Sumw2();
a4ae02cd 445
a220990f 446 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
447 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
448 //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
449 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.5648,2.1648);
450 //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
451 tmpBt->Sumw2();
452 //tmpBl->Sumw2();
453 tmpRt->Sumw2();
454 //tmpRl->Sumw2();
a4ae02cd 455
a220990f 456 fOutputMass->Add(tmpSt);
457 fOutputMass->Add(tmpRt);
458 fOutputMass->Add(tmpBt);
459
460 }
ea0d8716 461
a8ce111e 462 //mass
463 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.5648,2.1648);
a220990f 464 //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
a8ce111e 465 tmpMt->Sumw2();
a220990f 466 //tmpMl->Sumw2();
a8ce111e 467 //distribution w/o cuts large range
468 // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
469
470 fOutputMass->Add(tmpMt);
471
472 if(fSys==0){ //histograms filled only in pp to save time in PbPb
473 if(fFillVarHists){
a8ce111e 474
a220990f 475 if(fReadMC){
476 // pT
477 namedistr="hptpiS_";
478 namedistr+=i;
479 TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
480
481 namedistr="hptKS_";
482 namedistr+=i;
483 TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
484
485 // costhetastar
486 namedistr="hcosthetastarS_";
487 namedistr+=i;
488 TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
489
490 //pT no mass cut
491
492 namedistr="hptpiSnoMcut_";
493 namedistr+=i;
494 TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
495
496 namedistr="hptKSnoMcut_";
497 namedistr+=i;
498 TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
499
500 fDistr->Add(hptpiS);
501 fDistr->Add(hptKS);
502 fDistr->Add(hcosthetastarS);
503
504 fDistr->Add(hptpiSnoMcut);
505 fDistr->Add(hptKSnoMcut);
506
507 // costhetapoint vs d0 or d0d0
508 namedistr="hcosthpointd0S_";
509 namedistr+=i;
510 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);
511 namedistr="hcosthpointd0d0S_";
512 namedistr+=i;
513 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);
514
515 fDistr->Add(hcosthpointd0S);
516 fDistr->Add(hcosthpointd0d0S);
517
518 //to compare with AliAnalysisTaskCharmFraction
519 TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.5648,2.1648);
520 TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
521 tmpS27t->Sumw2();
522 tmpS27l->Sumw2();
523
524 fOutputMass->Add(tmpS27t);
525 fOutputMass->Add(tmpS27l);
526
527 }
528
529 // pT
530 namedistr="hptB_";
a8ce111e 531 namedistr+=i;
a220990f 532 TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
a8ce111e 533
a220990f 534 // costhetastar
535 namedistr="hcosthetastarB_";
a8ce111e 536 namedistr+=i;
a220990f 537 TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
a8ce111e 538
a220990f 539 //pT no mass cut
a8ce111e 540 namedistr="hptB1prongnoMcut_";
541 namedistr+=i;
542 TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
543
544 namedistr="hptB2prongsnoMcut_";
545 namedistr+=i;
546 TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
547
a220990f 548 fDistr->Add(hptB);
549 fDistr->Add(hcosthetastarB);
550
a8ce111e 551 fDistr->Add(hptB1pnoMcut);
552 fDistr->Add(hptB2pnoMcut);
553
554 //impact parameter of negative/positive track
555 namedistr="hd0p0B_";
556 namedistr+=i;
557 TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
558
559 namedistr="hd0p1B_";
560 namedistr+=i;
561 TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
562
563 //impact parameter corrected for strangeness
564 namedistr="hd0moresB_";
565 namedistr+=i;
566 TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
0108fa62 567
a8ce111e 568 namedistr="hd0d0moresB_";
569 namedistr+=i;
570 TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
571
572
573 namedistr="hcosthetapointmoresB_";
574 namedistr+=i;
575 TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
576
577 // costhetapoint vs d0 or d0d0
a8ce111e 578 namedistr="hcosthpointd0B_";
579 namedistr+=i;
580 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);
581
a8ce111e 582 namedistr="hcosthpointd0d0B_";
583 namedistr+=i;
584 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);
585
586 fDistr->Add(hd0p0B);
587 fDistr->Add(hd0p1B);
588
589 fDistr->Add(hd0moresB);
590 fDistr->Add(hd0d0moresB);
591 fDistr->Add(hcosthetapointmoresB);
592
a220990f 593
a8ce111e 594 fDistr->Add(hcosthpointd0B);
595
a220990f 596
a8ce111e 597 fDistr->Add(hcosthpointd0d0B);
598 }
34137226 599
a8ce111e 600 } //end pp histo
7646d6da 601 }
a4ae02cd 602
a8ce111e 603
604 //for Like sign analysis
605
9af24f46 606 if(fArray==1){
607 namedistr="hpospair";
82487ae7 608 TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
9af24f46 609 namedistr="hnegpair";
82487ae7 610 TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
818c1271 611 fOutputMass->Add(hpospair);
612 fOutputMass->Add(hnegpair);
9af24f46 613 }
449b1302 614
5b2e5fae 615 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
4e61a020 616
a220990f 617 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", 16,-0.5,15.5);
34137226 618
619 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
4e61a020 620 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
818c1271 621 if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
622 else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
6b3e3c78 623 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
700e80e0 624 fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
625 fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
a8ce111e 626 if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
627 if(fFillVarHists){
628 fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
629 fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
630 fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
631 fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
632 }
633 if(fReadMC && fSys==0){
634 fNentries->GetXaxis()->SetBinLabel(12,"K");
635 fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
636 }
700e80e0 637 fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
a3aa1279 638 fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
a220990f 639 if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
34137226 640 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
a4ae02cd 641
700e80e0 642 fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
a96083b9 643
40445ada 644 // Post the data
ea0d8716 645 PostData(1,fOutputMass);
a8ce111e 646 if(fFillVarHists) PostData(2,fDistr);
40445ada 647 PostData(3,fNentries);
700e80e0 648 PostData(5,fCounter);
49061176 649 return;
650}
651
652//________________________________________________________________________
653void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
654{
655 // Execute analysis for current event:
656 // heavy flavor candidates association to MC truth
a4ae02cd 657 //cout<<"I'm in UserExec"<<endl;
ea0d8716 658
659
a220990f 660 //cuts order
661 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
662 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
663 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
664 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
665 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
666 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
667 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
668 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
669 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
ea0d8716 670
671
49061176 672 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
feb73eca 673
b557eb43 674 TString bname;
feb73eca 675 if(fArray==0){ //D0 candidates
b557eb43 676 // load D0->Kpi candidates
feb73eca 677 //cout<<"D0 candidates"<<endl;
b557eb43 678 bname="D0toKpi";
feb73eca 679 } else { //LikeSign candidates
feb73eca 680 // load like sign candidates
b557eb43 681 //cout<<"LS candidates"<<endl;
682 bname="LikeSign2Prong";
683 }
684
685 TClonesArray *inputArray=0;
34137226 686
b557eb43 687 if(!aod && AODEvent() && IsStandardAOD()) {
688 // In case there is an AOD handler writing a standard AOD, use the AOD
689 // event in memory rather than the input (ESD) event.
690 aod = dynamic_cast<AliAODEvent*> (AODEvent());
691 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
692 // have to taken from the AOD event hold by the AliAODExtension
693 AliAODHandler* aodHandler = (AliAODHandler*)
694 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
34137226 695
b557eb43 696 if(aodHandler->GetExtensions()) {
697 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
34137226 698 AliAODEvent* aodFromExt = ext->GetAOD();
b557eb43 699 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
feb73eca 700 }
dc222f77 701 } else if(aod) {
b557eb43 702 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
703 }
feb73eca 704
b557eb43 705
dc222f77 706 if(!inputArray || !aod) {
b557eb43 707 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
708 return;
49061176 709 }
710
6b3e3c78 711 // fix for temporary bug in ESDfilter
7c23877d 712 // the AODs with null vertex pointer didn't pass the PhysSel
5806c290 713 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
714
7c23877d 715
ce39f0ac 716 TClonesArray *mcArray = 0;
717 AliAODMCHeader *mcHeader = 0;
718
719 if(fReadMC) {
720 // load MC particles
721 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
722 if(!mcArray) {
723 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
724 return;
725 }
40445ada 726
ce39f0ac 727 // load MC header
728 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
729 if(!mcHeader) {
730 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
731 return;
732 }
49061176 733 }
734
b557eb43 735 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
40445ada 736
a4ae02cd 737 //histogram filled with 1 for every AOD
34137226 738 fNentries->Fill(0);
a3aa1279 739 fCounter->StoreEvent(aod,fReadMC);
740
741 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
742 TString trigclass=aod->GetFiredTriggerClasses();
743 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
744
d52f7b50 745 if(!fCuts->IsEventSelected(aod)) {
746 if(fCuts->GetWhyRejection()==1) // rejected for pileup
700e80e0 747 fNentries->Fill(13);
a220990f 748 if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
d52f7b50 749 return;
750 }
40445ada 751
752 // AOD primary vertex
753 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
b272aebf 754
40445ada 755 Bool_t isGoodVtx=kFALSE;
b272aebf 756
a220990f 757 //vtx1->Print();
40445ada 758 TString primTitle = vtx1->GetTitle();
759 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
760 isGoodVtx=kTRUE;
761 fNentries->Fill(3);
762 }
763
75638da0 764 // loop over candidates
765 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4e61a020 766 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
75638da0 767
6b3e3c78 768 // FILE *f=fopen("4display.txt","a");
769 // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
a96083b9 770 Int_t nSelectedloose=0,nSelectedtight=0;
75638da0 771 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
75638da0 772 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
818c1271 773
774 if(d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
775 fNentries->Fill(2);
776 continue; //skip the D0 from Dstar
777 }
a220990f 778
d7688946 779 Bool_t unsetvtx=kFALSE;
780 if(!d->GetOwnPrimaryVtx()) {
781 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
782 unsetvtx=kTRUE;
783 }
ea0d8716 784
4e61a020 785
4e61a020 786 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
a96083b9 787 nSelectedloose++;
d7688946 788 nSelectedtight++;
a8ce111e 789 if(fSys==0){
790 if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
791 }
6b3e3c78 792 Int_t ptbin=fCuts->PtBin(d->Pt());
700e80e0 793 if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
d7688946 794 fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod); //selected
795 if(fFillVarHists) {
a8ce111e 796 //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
a220990f 797 fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
798 fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
799 //check daughters
800 if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
801 AliDebug(1,"at least one daughter not found!");
802 fNentries->Fill(5);
803 fDaughterTracks.Clear();
804 continue;
805 }
806 //}
d7688946 807 FillVarHists(aod,d,mcArray,fCuts,fDistr);
808 }
a8ce111e 809
d7688946 810 FillMassHists(d,mcArray,fCuts,fOutputMass);
75638da0 811 }
ea0d8716 812
d7688946 813 fDaughterTracks.Clear();
814 if(unsetvtx) d->UnsetOwnPrimaryVtx();
75638da0 815 } //end for prongs
a96083b9 816 fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
817 fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
40445ada 818 // Post the data
ea0d8716 819 PostData(1,fOutputMass);
a8ce111e 820 if(fFillVarHists) PostData(2,fDistr);
40445ada 821 PostData(3,fNentries);
700e80e0 822 PostData(5,fCounter);
40445ada 823 return;
824}
b272aebf 825
40445ada 826//____________________________________________________________________________
4e61a020 827void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
40445ada 828 //
829 // function used in UserExec to fill variable histograms:
830 //
b272aebf 831
b272aebf 832
40445ada 833 Int_t pdgDgD0toKpi[2]={321,211};
834 Int_t lab=-9999;
835 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)
836 //Double_t pt = d->Pt(); //mother pt
a8ce111e 837 Int_t isSelectedPID=3;
838 if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
700e80e0 839 if (isSelectedPID==0)fNentries->Fill(7);
840 if (isSelectedPID==1)fNentries->Fill(8);
841 if (isSelectedPID==2)fNentries->Fill(9);
842 if (isSelectedPID==3)fNentries->Fill(10);
a220990f 843 //fNentries->Fill(8+isSelectedPID);
3cc4604b 844
a8ce111e 845 if(fCutOnDistr && !fIsSelectedCandidate) return;
d7688946 846 //printf("\nif no cuts or cuts passed\n");
847
848
849 //add distr here
850 UInt_t pdgs[2];
851
852 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
853 pdgs[0]=211;
854 pdgs[1]=321;
855 Double_t minvD0 = part->InvMassD0();
856 pdgs[1]=211;
857 pdgs[0]=321;
858 Double_t minvD0bar = part->InvMassD0bar();
859 //cout<<"inside mass cut"<<endl;
0108fa62 860
4e61a020 861 Double_t invmasscut=0.03;
862
40445ada 863 TString fillthispi="",fillthisK="",fillthis="";
b272aebf 864
ea0d8716 865 Int_t ptbin=cuts->PtBin(part->Pt());
6b3e3c78 866
4e61a020 867 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
868 dz1[0]=-99; dz2[0]=-99;
9af24f46 869 Double_t d0[2];
870 Double_t decl[2]={-99,-99};
871 Bool_t recalcvtx=kFALSE;
872
d7688946 873
874
9af24f46 875 if(fCuts->GetIsPrimaryWithoutDaughters()){
876 recalcvtx=kTRUE;
877 //recalculate vertex
878 AliAODVertex *vtxProp=0x0;
d7688946 879 vtxProp=GetPrimaryVtxSkipped(aod);
9af24f46 880 if(vtxProp) {
881 part->SetOwnPrimaryVtx(vtxProp);
882 //Bool_t unsetvtx=kTRUE;
883 //Calculate d0 for daughter tracks with Vtx Skipped
d7688946 884 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
885 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
9af24f46 886 esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
887 esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
888 delete vtxProp; vtxProp=NULL;
889 delete esdtrack1;
890 delete esdtrack2;
891 }
892
893 d0[0]=dz1[0];
894 d0[1]=dz2[0];
895
a8ce111e 896 decl[0]=part->DecayLength2();
897 decl[1]=part->NormalizedDecayLength2();
9af24f46 898 part->UnsetOwnPrimaryVtx();
899
4e61a020 900 }
901
d7688946 902 Double_t cosThetaStarD0 = 99;
903 Double_t cosThetaStarD0bar = 99;
904 Double_t cosPointingAngle = 99;
a8ce111e 905 Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
906 Double_t decayLength2 = -1, decayLengthxy=-1;
d7688946 907 Double_t ptProng[2]={-99,-99};
a8ce111e 908 Double_t d0Prong[2]={-99,-99};
d7688946 909
0108fa62 910
a220990f 911 //disable the PID
912 if(!fUsePid4Distr) isSelectedPID=0;
913 if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
9af24f46 914
a220990f 915 //check pdg of the prongs
916 AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
917 AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
6b3e3c78 918
a220990f 919 if(!prong0 || !prong1) {
920 return;
921 }
922 Int_t labprong[2];
923 if(fReadMC){
924 labprong[0]=prong0->GetLabel();
925 labprong[1]=prong1->GetLabel();
926 }
927 AliAODMCParticle *mcprong=0;
928 Int_t pdgProng[2]={0,0};
929
930 for (Int_t iprong=0;iprong<2;iprong++){
931 if(fReadMC && labprong[iprong]>=0) {
932 mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
933 pdgProng[iprong]=mcprong->GetPdgCode();
0108fa62 934 }
a220990f 935 }
6b3e3c78 936
a220990f 937 if(fSys==0){
938 //no mass cut ditributions: ptbis
40445ada 939
a220990f 940 fillthispi="hptpiSnoMcut_";
941 fillthispi+=ptbin;
a8ce111e 942
a220990f 943 fillthisK="hptKSnoMcut_";
944 fillthisK+=ptbin;
a8ce111e 945
a220990f 946 if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
947 || (isSelectedPID==1 || isSelectedPID==3)){
948 ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
949 ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
950 }
a8ce111e 951
a220990f 952 if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
953 || (isSelectedPID==2 || isSelectedPID==3)){
954 ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
955 ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
a8ce111e 956 }
a220990f 957 }
a8ce111e 958
a220990f 959 //no mass cut ditributions: mass
960 fillthis="hMassS_";
961 fillthis+=ptbin;
40445ada 962
a220990f 963 if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
964 || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
965 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
966 }
967 else { //D0bar
968 if(fReadMC || (!fReadMC && isSelectedPID > 1))
40445ada 969 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
a220990f 970 }
4e61a020 971
a220990f 972 //apply cut on invariant mass on the pair
973 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
d7688946 974
a220990f 975 cosThetaStarD0 = part->CosThetaStarD0();
976 cosThetaStarD0bar = part->CosThetaStarD0bar();
977 cosPointingAngle = part->CosPointingAngle();
978 normalizedDecayLength2 = part->NormalizedDecayLength2();
979 decayLength2 = part->DecayLength2();
980 decayLengthxy = part->DecayLengthXY();
981 normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
a8ce111e 982
a220990f 983 ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
a220990f 984 d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
a8ce111e 985
a220990f 986 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
987 for (Int_t iprong=0; iprong<2; iprong++){
988 AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
989 if (fReadMC) labprong[iprong]=prong->GetLabel();
40445ada 990
a220990f 991 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
992 Int_t pdgprong=0;
993 if(fReadMC && labprong[iprong]>=0) {
994 mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
995 pdgprong=mcprong->GetPdgCode();
996 }
6b3e3c78 997
a220990f 998 Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
6b3e3c78 999
a220990f 1000 if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1001 //cout<<"pi"<<endl;
1002
1003 if(fSys==0){
40445ada 1004 fillthispi="hptpiS_";
1005 fillthispi+=ptbin;
d7688946 1006 ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
a220990f 1007 }
1008 fillthispi="hd0piS_";
1009 fillthispi+=ptbin;
1010 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1011 if(recalcvtx) {
4e61a020 1012
a220990f 1013 fillthispi="hd0vpiS_";
1014 fillthispi+=ptbin;
1015 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
b272aebf 1016 }
a220990f 1017
1018 }
40445ada 1019
a220990f 1020 if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1021 //cout<<"kappa"<<endl;
1022 if(fSys==0){
40445ada 1023 fillthisK="hptKS_";
1024 fillthisK+=ptbin;
d7688946 1025 ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
a220990f 1026 }
1027 fillthisK="hd0KS_";
1028 fillthisK+=ptbin;
1029 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1030 if (recalcvtx){
1031 fillthisK="hd0vKS_";
40445ada 1032 fillthisK+=ptbin;
a220990f 1033 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
b272aebf 1034 }
a220990f 1035 }
b272aebf 1036
a220990f 1037 if(fSys==0){
1038 fillthis="hcosthpointd0S_";
1039 fillthis+=ptbin;
1040 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1041 }
1042 } //end loop on prongs
b272aebf 1043
a220990f 1044 fillthis="hdcaS_";
1045 fillthis+=ptbin;
1046 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
4e61a020 1047
a220990f 1048 fillthis="hcosthetapointS_";
1049 fillthis+=ptbin;
1050 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
4e61a020 1051
a220990f 1052 fillthis="hcosthetapointxyS_";
1053 fillthis+=ptbin;
1054 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
a8ce111e 1055
4e61a020 1056
a220990f 1057 fillthis="hd0d0S_";
1058 fillthis+=ptbin;
1059 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
4e61a020 1060
a220990f 1061 fillthis="hdeclS_";
1062 fillthis+=ptbin;
1063 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
a8ce111e 1064
a220990f 1065 fillthis="hnormdeclS_";
1066 fillthis+=ptbin;
1067 ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1068
1069 fillthis="hdeclxyS_";
1070 fillthis+=ptbin;
1071 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1072
1073 fillthis="hnormdeclxyS_";
1074 fillthis+=ptbin;
1075 ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1076
1077 fillthis="hdeclxyd0d0S_";
1078 fillthis+=ptbin;
1079 ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1080
1081 fillthis="hnormdeclxyd0d0S_";
1082 fillthis+=ptbin;
1083 ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
a8ce111e 1084
a220990f 1085 if(recalcvtx) {
1086 fillthis="hdeclvS_";
a8ce111e 1087 fillthis+=ptbin;
a220990f 1088 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
a2121012 1089
a220990f 1090 fillthis="hnormdeclvS_";
1091 fillthis+=ptbin;
1092 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
a2121012 1093
a220990f 1094 fillthis="hd0d0vS_";
1095 fillthis+=ptbin;
1096 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1097 }
a8ce111e 1098
a220990f 1099 if(fSys==0){
1100 fillthis="hcosthetastarS_";
1101 fillthis+=ptbin;
1102 if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1103 else {
1104 if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1105 if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
a8ce111e 1106 }
1107
a220990f 1108 fillthis="hcosthpointd0d0S_";
1109 fillthis+=ptbin;
1110 ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1111 }
a8ce111e 1112
a220990f 1113 } //end mass cut
4e61a020 1114
a220990f 1115 } else{ //Background or LS
1116 //if(!fReadMC){
1117 //cout<<"is background"<<endl;
40445ada 1118
a220990f 1119 //no mass cut distributions: mass, ptbis
1120 fillthis="hMassB_";
1121 fillthis+=ptbin;
a8ce111e 1122
a220990f 1123 if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1124 if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1125 if(fSys==0){
1126 fillthis="hptB1prongnoMcut_";
1127 fillthis+=ptbin;
40445ada 1128
a220990f 1129 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
40445ada 1130
a220990f 1131 fillthis="hptB2prongsnoMcut_";
1132 fillthis+=ptbin;
1133 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
818c1271 1134 ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
a220990f 1135 }
1136 //apply cut on invariant mass on the pair
1137 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1138 if(fSys==0){
818c1271 1139 ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
d7688946 1140 cosThetaStarD0 = part->CosThetaStarD0();
1141 cosThetaStarD0bar = part->CosThetaStarD0bar();
a220990f 1142 }
1143
1144 cosPointingAngle = part->CosPointingAngle();
1145 normalizedDecayLength2 = part->NormalizedDecayLength2();
1146 decayLength2 = part->DecayLength2();
1147 decayLengthxy = part->DecayLengthXY();
1148 normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1149 d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1150
1151
1152 AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1153 if(!prongg) {
1154 if(fDebug>2) cout<<"No daughter found";
1155 return;
1156 }
1157 else{
1158 if(fArray==1){
1159 if(prongg->Charge()==1) {
1160 //fTotPosPairs[ptbin]++;
1161 ((TH1F*)fDistr->FindObject("hpospair"))->Fill(ptbin);
1162 } else {
1163 //fTotNegPairs[ptbin]++;
1164 ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(ptbin);
449b1302 1165 }
0108fa62 1166 }
a220990f 1167 }
40445ada 1168
a220990f 1169 fillthis="hd0B_";
1170 fillthis+=ptbin;
1171 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1172 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
0108fa62 1173
a220990f 1174 if(fReadMC){
1175 Int_t pdgMother[2]={0,0};
1176 Double_t factor[2]={1,1};
4e61a020 1177
a220990f 1178 for(Int_t iprong=0;iprong<2;iprong++){
1179 AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1180 lab=prong->GetLabel();
1181 if(lab>=0){
1182 AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1183 if(mcprong){
1184 Int_t labmom=mcprong->GetMother();
1185 if(labmom>=0){
1186 AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1187 if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
6b3e3c78 1188 }
1189 }
a220990f 1190 }
6b3e3c78 1191
a220990f 1192 if(fSys==0){
a8ce111e 1193
a220990f 1194 fillthis="hd0moresB_";
1195 fillthis+=ptbin;
6b3e3c78 1196
a220990f 1197 if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1198 if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1199 else factor[iprong]=1./.6;
1200 fNentries->Fill(11);
1201 }
6b3e3c78 1202
a220990f 1203 if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1204 factor[iprong]=1./0.25;
1205 fNentries->Fill(12);
1206 }
1207 fillthis="hd0moresB_";
1208 fillthis+=ptbin;
6b3e3c78 1209
a220990f 1210 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
6b3e3c78 1211
a220990f 1212 if(recalcvtx){
1213 fillthis="hd0vmoresB_";
1214 fillthis+=ptbin;
1215 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
9af24f46 1216 }
a220990f 1217 }
1218 } //loop on prongs
6b3e3c78 1219
a220990f 1220 if(fSys==0){
6b3e3c78 1221 fillthis="hd0d0moresB_";
1222 fillthis+=ptbin;
1223 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1224
a8ce111e 1225 fillthis="hcosthetapointmoresB_";
1226 fillthis+=ptbin;
1227 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1228
9af24f46 1229 if(recalcvtx){
1230 fillthis="hd0d0vmoresB_";
1231 fillthis+=ptbin;
1232 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1233 }
a220990f 1234 }
1235 } //readMC
ce39f0ac 1236
a220990f 1237 if(fSys==0){
1238 //normalise pt distr to half afterwards
1239 fillthis="hptB_";
40445ada 1240 fillthis+=ptbin;
a220990f 1241 ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1242 ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
34137226 1243
40445ada 1244 fillthis="hcosthetastarB_";
1245 fillthis+=ptbin;
d7688946 1246 if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1247 if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
40445ada 1248
a220990f 1249
1250 fillthis="hd0p0B_";
40445ada 1251 fillthis+=ptbin;
a220990f 1252 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1253 fillthis="hd0p1B_";
1254 fillthis+=ptbin;
1255 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1256
1257 fillthis="hcosthpointd0d0B_";
1258 fillthis+=ptbin;
1259 ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1260
1261 fillthis="hcosthpointd0B_";
1262 fillthis+=ptbin;
1263 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1264 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1265
40445ada 1266
9af24f46 1267 if(recalcvtx){
a220990f 1268
1269 fillthis="hd0vp0B_";
1270 fillthis+=ptbin;
1271 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1272 fillthis="hd0vp1B_";
1273 fillthis+=ptbin;
1274 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1275
1276 fillthis="hd0vB_";
a2121012 1277 fillthis+=ptbin;
a220990f 1278 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1279 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1280
a2121012 1281 }
1282
a220990f 1283 }
40445ada 1284
a220990f 1285 fillthis="hdcaB_";
1286 fillthis+=ptbin;
1287 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
a2121012 1288
a220990f 1289 fillthis="hd0d0B_";
1290 fillthis+=ptbin;
1291 ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
a8ce111e 1292
a220990f 1293 if(recalcvtx){
1294 fillthis="hd0d0vB_";
a8ce111e 1295 fillthis+=ptbin;
a220990f 1296 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1297 }
a8ce111e 1298
a220990f 1299 fillthis="hcosthetapointB_";
1300 fillthis+=ptbin;
1301 ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
a2121012 1302
a220990f 1303 fillthis="hcosthetapointxyB_";
1304 fillthis+=ptbin;
1305 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
9af24f46 1306
a220990f 1307 fillthis="hdeclB_";
1308 fillthis+=ptbin;
1309 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
a2121012 1310
a220990f 1311 fillthis="hnormdeclB_";
1312 fillthis+=ptbin;
1313 ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
9af24f46 1314
a220990f 1315 fillthis="hdeclxyB_";
1316 fillthis+=ptbin;
1317 ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
a8ce111e 1318
a220990f 1319 fillthis="hnormdeclxyB_";
1320 fillthis+=ptbin;
1321 ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1322
1323 fillthis="hdeclxyd0d0B_";
1324 fillthis+=ptbin;
1325 ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1326
1327 fillthis="hnormdeclxyd0d0B_";
1328 fillthis+=ptbin;
1329 ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1330
1331 if(recalcvtx) {
1332
1333 fillthis="hdeclvB_";
1334 fillthis+=ptbin;
1335 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1336
1337 fillthis="hnormdeclvB_";
1338 fillthis+=ptbin;
1339 ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1340
1341
1342 }
1343 }//mass cut
1344 }//else (background)
d7688946 1345
49061176 1346 return;
1347}
ea0d8716 1348//____________________________________________________________________________
1349
d7688946 1350void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
49061176 1351 //
40445ada 1352 // function used in UserExec to fill mass histograms:
49061176 1353 //
9de8c723 1354
1355
a8ce111e 1356 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
9de8c723 1357
d7688946 1358 //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
feb73eca 1359
ea0d8716 1360 //cout<<"check cuts = "<<endl;
1361 //cuts->PrintAll();
d7688946 1362 if (!fIsSelectedCandidate){
ea0d8716 1363 //cout<<"Not Selected"<<endl;
6b3e3c78 1364 //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
ea0d8716 1365 return;
1366 }
9de8c723 1367
6b3e3c78 1368
4e61a020 1369 if(fDebug>2) cout<<"Candidate selected"<<endl;
a41f6fad 1370
ea0d8716 1371 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1372 //printf("SELECTED\n");
1373 Int_t ptbin=cuts->PtBin(part->Pt());
9de8c723 1374
818c1271 1375 // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1376 // if(!prong) {
1377 // AliDebug(2,"No daughter found");
1378 // return;
1379 // }
1380 // else{
a8ce111e 1381 // if(prong->Charge()==1) {
1382 // ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
1383 // //fTotPosPairs[ptbin]++;
1384 // } else {
1385 // ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
1386 // //fTotNegPairs[ptbin]++;
1387 // }
818c1271 1388 // }
700e80e0 1389
9af24f46 1390 // for(Int_t it=0;it<2;it++){
1391
1392 // //request on spd points to be addes
1393 // if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
1394 // FILE *f=fopen("4display.txt","a");
1395 // 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());
1396 // fclose(f);
1397 // //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
1398 // }
1399 // }
6b3e3c78 1400
ea0d8716 1401 TString fillthis="";
1402 Int_t pdgDgD0toKpi[2]={321,211};
1403 Int_t labD0=-1;
1404 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)
1405
1406 //count candidates selected by cuts
1407 fNentries->Fill(1);
1408 //count true D0 selected by cuts
1409 if (fReadMC && labD0>=0) fNentries->Fill(2);
ea0d8716 1410
d7688946 1411 if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
a220990f 1412 if(fReadMC){
1413 if(labD0>=0) {
1414 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1415
1416 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1417 Int_t pdgD0 = partD0->GetPdgCode();
1418 //cout<<"pdg = "<<pdgD0<<endl;
1419 if (pdgD0==421){ //D0
1420 //cout<<"Fill S with D0"<<endl;
1421 fillthis="histSgn_";
1422 fillthis+=ptbin;
1423 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1424 if(fSys==0){
1425 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1426 fillthis="histSgn27_";
1427 fillthis+=ptbin;
1428 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1429 }
1430 }
1431 } else{ //it was a D0bar
1432 fillthis="histRfl_";
a8ce111e 1433 fillthis+=ptbin;
1434 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1435 }
a220990f 1436 } else {//background
1437 fillthis="histBkg_";
a4ae02cd 1438 fillthis+=ptbin;
1439 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1440 }
a220990f 1441
1442 }else{
1443 fillthis="histMass_";
a4ae02cd 1444 fillthis+=ptbin;
a220990f 1445 //cout<<"Filling "<<fillthis<<endl;
1446
1447 //printf("Fill mass with D0");
ea0d8716 1448 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1449 }
a220990f 1450
ea0d8716 1451 }
d7688946 1452 if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
a220990f 1453 if(fReadMC){
1454 if(labD0>=0) {
1455 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1456 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1457 Int_t pdgD0 = partD0->GetPdgCode();
1458 //cout<<" pdg = "<<pdgD0<<endl;
1459 if (pdgD0==-421){ //D0bar
1460 fillthis="histSgn_";
1461 fillthis+=ptbin;
1462 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1463 // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1464 // fillthis="histSgn27_";
1465 // fillthis+=ptbin;
1466 // ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1467 // }
9de8c723 1468
a4ae02cd 1469
a220990f 1470 } else{
1471 fillthis="histRfl_";
1472 fillthis+=ptbin;
1473 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1474 }
1475 } else {//background or LS
1476 fillthis="histBkg_";
a4ae02cd 1477 fillthis+=ptbin;
1478 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1479 }
a220990f 1480 }else{
1481 fillthis="histMass_";
ea0d8716 1482 fillthis+=ptbin;
a220990f 1483 //printf("Fill mass with D0bar");
1484 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1485
a4ae02cd 1486 }
ea0d8716 1487 }
a4ae02cd 1488
40445ada 1489 return;
49061176 1490}
4e61a020 1491
1492//__________________________________________________________________________
d7688946 1493AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
4e61a020 1494 //Calculate the primary vertex w/o the daughter tracks of the candidate
1495
4e61a020 1496 Int_t skipped[2];
1497 Int_t nTrksToSkip=2;
d7688946 1498 AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
4e61a020 1499 if(!dgTrack){
1500 AliDebug(2,"no daughter found!");
1501 return 0x0;
1502 }
1503 skipped[0]=dgTrack->GetID();
d7688946 1504 dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
4e61a020 1505 if(!dgTrack){
1506 AliDebug(2,"no daughter found!");
1507 return 0x0;
1508 }
1509 skipped[1]=dgTrack->GetID();
1510
6b3e3c78 1511 AliESDVertex *vertexESD=0x0;
1512 AliAODVertex *vertexAOD=0x0;
1513 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1514
4e61a020 1515 //
1516 vertexer->SetSkipTracks(nTrksToSkip,skipped);
4e61a020 1517 vertexer->SetMinClusters(4);
a2121012 1518 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
4e61a020 1519 if(!vertexESD) return vertexAOD;
1520 if(vertexESD->GetNContributors()<=0) {
1521 AliDebug(2,"vertexing failed");
1522 delete vertexESD; vertexESD=NULL;
1523 return vertexAOD;
1524 }
1525
1526 delete vertexer; vertexer=NULL;
1527
1528
1529 // convert to AliAODVertex
1530 Double_t pos[3],cov[6],chi2perNDF;
1531 vertexESD->GetXYZ(pos); // position
1532 vertexESD->GetCovMatrix(cov); //covariance matrix
1533 chi2perNDF = vertexESD->GetChi2toNDF();
1534 delete vertexESD; vertexESD=NULL;
1535
1536 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1537 return vertexAOD;
1538
1539}
1540
1541
49061176 1542//________________________________________________________________________
1543void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1544{
1545 // Terminate analysis
1546 //
1547 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1548
6321ee46 1549
ea0d8716 1550 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1551 if (!fOutputMass) {
1552 printf("ERROR: fOutputMass not available\n");
a4ae02cd 1553 return;
1554 }
a8ce111e 1555 if(fFillVarHists){
1556 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1557 if (!fDistr) {
1558 printf("ERROR: fDistr not available\n");
1559 return;
1560 }
a41f6fad 1561 }
a8ce111e 1562
40445ada 1563 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
5b2e5fae 1564
40445ada 1565 if(!fNentries){
1566 printf("ERROR: fNEntries not available\n");
1567 return;
1568 }
700e80e0 1569 fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
1570 if(!fCuts){
1571 printf("ERROR: fCuts not available\n");
34137226 1572 return;
1573 }
700e80e0 1574 fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
a96083b9 1575 if (!fCounter) {
1576 printf("ERROR: fCounter not available\n");
1577 return;
1578 }
700e80e0 1579 Int_t nptbins=fCuts->GetNPtBins();
1580 for(Int_t ipt=0;ipt<nptbins;ipt++){
4e61a020 1581
a8ce111e 1582 if(fArray==1 && fFillVarHists){
700e80e0 1583 fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fDistr->FindObject("hpospair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)*((TH1F*)fDistr->FindObject("hnegpair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)); //after cuts
feb73eca 1584
1585
ea0d8716 1586 if(fLsNormalization>1e-6) {
9de8c723 1587
feb73eca 1588 TString massName="histMass_";
ea0d8716 1589 massName+=ipt;
1590 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1591
feb73eca 1592 }
40445ada 1593
feb73eca 1594
700e80e0 1595 fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fDistr->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fDistr->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
1596 //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
40445ada 1597
ea0d8716 1598 if(fLsNormalization>1e-6) {
40445ada 1599
1600 TString nameDistr="hptB_";
ea0d8716 1601 nameDistr+=ipt;
40445ada 1602 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1603 nameDistr="hdcaB_";
ea0d8716 1604 nameDistr+=ipt;
40445ada 1605 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1606 nameDistr="hcosthetastarB_";
ea0d8716 1607 nameDistr+=ipt;
40445ada 1608 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1609 nameDistr="hd0B_";
ea0d8716 1610 nameDistr+=ipt;
40445ada 1611 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1612 nameDistr="hd0d0B_";
ea0d8716 1613 nameDistr+=ipt;
40445ada 1614 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1615 nameDistr="hcosthetapointB_";
ea0d8716 1616 nameDistr+=ipt;
40445ada 1617 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
a220990f 1618 if(fSys==0){
1619 nameDistr="hcosthpointd0d0B_";
1620 nameDistr+=ipt;
1621 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1622 }
40445ada 1623 }
feb73eca 1624 }
1625 }
4e61a020 1626 TString cvname,cstname;
527f330b 1627
1628 if (fArray==0){
1629 cvname="D0invmass";
4e61a020 1630 cstname="cstat0";
1631 } else {
1632 cvname="LSinvmass";
1633 cstname="cstat1";
1634 }
527f330b 1635
34137226 1636 TCanvas *cMass=new TCanvas(cvname,cvname);
1637 cMass->cd();
ea0d8716 1638 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
527f330b 1639
4e61a020 1640 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
34137226 1641 cStat->cd();
1642 cStat->SetGridy();
1643 fNentries->Draw("htext0");
527f330b 1644
4e61a020 1645 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1646 // ccheck->cd();
1647
49061176 1648 return;
1649}
1650