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