]> 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
Histogram range changed
[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
16/////////////////////////////////////////////////////////////
17//
18// AliAnalysisTaskSE for D0 candidates invariant mass histogram
a41f6fad 19// and comparison with the MC truth and cut variables distributions.
49061176 20//
21// Authors: A.Dainese, andrea.dainese@lnl.infn.it
feb73eca 22// Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
23// Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
49061176 24/////////////////////////////////////////////////////////////
25
26#include <Riostream.h>
27#include <TClonesArray.h>
527f330b 28#include <TCanvas.h>
49061176 29#include <TNtuple.h>
30#include <TList.h>
31#include <TH1F.h>
a41f6fad 32#include <TH2F.h>
33#include <TDatabasePDG.h>
49061176 34
5b2e5fae 35#include <AliAnalysisDataSlot.h>
36#include <AliAnalysisDataContainer.h>
b557eb43 37#include "AliAnalysisManager.h"
34137226 38#include "AliESDtrack.h"
4e61a020 39#include "AliVertexerTracks.h"
b557eb43 40#include "AliAODHandler.h"
49061176 41#include "AliAODEvent.h"
42#include "AliAODVertex.h"
43#include "AliAODTrack.h"
44#include "AliAODMCHeader.h"
45#include "AliAODMCParticle.h"
46#include "AliAODRecoDecayHF2Prong.h"
47#include "AliAODRecoCascadeHF.h"
48#include "AliAnalysisVertexingHF.h"
49#include "AliAnalysisTaskSE.h"
50#include "AliAnalysisTaskSED0Mass.h"
51
b557eb43 52
49061176 53ClassImp(AliAnalysisTaskSED0Mass)
54
55
56//________________________________________________________________________
57AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
58AliAnalysisTaskSE(),
ea0d8716 59fOutputMass(0),
60fDistr(0),
61fNentries(0),
34137226 62fChecks(0),
ea0d8716 63fCuts(0),
feb73eca 64fArray(0),
ce39f0ac 65fReadMC(0),
40445ada 66fCutOnDistr(0),
87020237 67fNPtBins(1),
4b4a1d25 68fTotPosPairs(0),
69fTotNegPairs(0),
feb73eca 70fLsNormalization(1.)
71
ea0d8716 72
49061176 73{
74 // Default constructor
75}
76
77//________________________________________________________________________
ea0d8716 78AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
49061176 79AliAnalysisTaskSE(name),
ea0d8716 80fOutputMass(0),
34137226 81fDistr(0),
ea0d8716 82fNentries(0),
34137226 83fChecks(0),
ea0d8716 84fCuts(0),
feb73eca 85fArray(0),
ce39f0ac 86fReadMC(0),
40445ada 87fCutOnDistr(0),
87020237 88fNPtBins(1),
4b4a1d25 89fTotPosPairs(0),
90fTotNegPairs(0),
feb73eca 91fLsNormalization(1.)
6321ee46 92
ea0d8716 93
49061176 94{
95 // Default constructor
87020237 96
97 fNPtBins=cuts->GetNPtBins();
98 fTotPosPairs=new Int_t[fNPtBins];
99 fTotNegPairs=new Int_t[fNPtBins];
100 for(Int_t i=0;i<fNPtBins;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
4b4a1d25 101
ea0d8716 102 fCuts=cuts;
103
104 // Output slot #1 writes into a TList container (mass with cuts)
49061176 105 DefineOutput(1,TList::Class()); //My private output
ea0d8716 106 // Output slot #2 writes into a TList container (distributions)
a4ae02cd 107 DefineOutput(2,TList::Class()); //My private output
ea0d8716 108 // Output slot #3 writes into a TH1F container (number of events)
a4ae02cd 109 DefineOutput(3,TH1F::Class()); //My private output
ea0d8716 110 // Output slot #4 writes into a TList container (quality check)
a41f6fad 111 DefineOutput(4,TList::Class()); //My private output
ea0d8716 112 // Output slot #5 writes into a TList container (cuts)
113 DefineOutput(5,AliRDHFCutsD0toKpi::Class()); //My private output
ea0d8716 114
49061176 115}
116
117//________________________________________________________________________
118AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
119{
ea0d8716 120 if (fOutputMass) {
121 delete fOutputMass;
122 fOutputMass = 0;
a4ae02cd 123 }
a41f6fad 124 if (fDistr) {
125 delete fDistr;
126 fDistr = 0;
127 }
34137226 128 if (fChecks) {
129 delete fChecks;
130 fChecks = 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 }
ea0d8716 140
ea0d8716 141
49061176 142}
143
144//________________________________________________________________________
145void AliAnalysisTaskSED0Mass::Init()
146{
147 // Initialization
148
149 if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
150
ea0d8716 151
ea0d8716 152 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
4e61a020 153 const char* nameoutput=GetOutputSlot(5)->GetContainer()->GetName();
154 copyfCuts->SetName(nameoutput);
ea0d8716 155 // Post the data
156 PostData(5,copyfCuts);
4e61a020 157
49061176 158
49061176 159 return;
160}
161
162//________________________________________________________________________
163void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
164{
165
166 // Create the output container
167 //
168 if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
169
170 // Several histograms are more conveniently managed in a TList
ea0d8716 171 fOutputMass = new TList();
172 fOutputMass->SetOwner();
173 fOutputMass->SetName("listMass");
49061176 174
a41f6fad 175 fDistr = new TList();
176 fDistr->SetOwner();
177 fDistr->SetName("distributionslist");
178
34137226 179 fChecks = new TList();
180 fChecks->SetOwner();
181 fChecks->SetName("checkHistograms");
182
0108fa62 183 TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
49061176 184
ea0d8716 185 for(Int_t i=0;i<fCuts->GetNPtBins();i++){
b272aebf 186
7646d6da 187 nameMass="histMass_";
ea0d8716 188 nameMass+=i;
9de8c723 189 nameSgn27="histSgn27_";
ea0d8716 190 nameSgn27+=i;
7646d6da 191 nameSgn="histSgn_";
ea0d8716 192 nameSgn+=i;
7646d6da 193 nameBkg="histBkg_";
ea0d8716 194 nameBkg+=i;
a4ae02cd 195 nameRfl="histRfl_";
ea0d8716 196 nameRfl+=i;
0108fa62 197 nameMassNocutsS="hMassS_";
ea0d8716 198 nameMassNocutsS+=i;
0108fa62 199 nameMassNocutsB="hMassB_";
ea0d8716 200 nameMassNocutsB+=i;
7646d6da 201
b272aebf 202 //histograms of cut variable distributions
203
204 // pT
205 namedistr="hptpiS_";
ea0d8716 206 namedistr+=i;
b272aebf 207 TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
208
209 namedistr="hptKS_";
ea0d8716 210 namedistr+=i;
b272aebf 211 TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
212
213 namedistr="hptB_";
ea0d8716 214 namedistr+=i;
b272aebf 215 TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
216
217 // pT no mass cut
218 namedistr="hptpiSnoMcut_";
ea0d8716 219 namedistr+=i;
b272aebf 220 TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
221
222 namedistr="hptKSnoMcut_";
ea0d8716 223 namedistr+=i;
b272aebf 224 TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
225
226 namedistr="hptB1prongnoMcut_";
ea0d8716 227 namedistr+=i;
b272aebf 228 TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
229
230 namedistr="hptB2prongsnoMcut_";
ea0d8716 231 namedistr+=i;
b272aebf 232 TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
233
234
235 // dca
236 namedistr="hdcaS_";
ea0d8716 237 namedistr+=i;
b272aebf 238 TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
239 namedistr="hdcaB_";
ea0d8716 240 namedistr+=i;
b272aebf 241 TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
242
243 // costhetastar
244 namedistr="hcosthetastarS_";
ea0d8716 245 namedistr+=i;
b272aebf 246 TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
247 namedistr="hcosthetastarB_";
ea0d8716 248 namedistr+=i;
b272aebf 249 TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
250
251 // impact parameter
252 namedistr="hd0piS_";
ea0d8716 253 namedistr+=i;
b272aebf 254 TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
255
256 namedistr="hd0KS_";
ea0d8716 257 namedistr+=i;
b272aebf 258 TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
4e61a020 259
b272aebf 260 namedistr="hd0B_";
ea0d8716 261 namedistr+=i;
4e61a020 262 TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
263
264 namedistr="hd0p0B_";
265 namedistr+=i;
266 TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
267
268 namedistr="hd0p1B_";
269 namedistr+=i;
270 TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
271
272 namedistr="hd0vpiS_";
273 namedistr+=i;
274 TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
275
276 namedistr="hd0vKS_";
277 namedistr+=i;
278 TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
279 namedistr="hd0vB_";
280 namedistr+=i;
281 TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
282
283 namedistr="hd0vp0B_";
284 namedistr+=i;
285 TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
286
287 namedistr="hd0vp1B_";
288 namedistr+=i;
289 TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
b272aebf 290
291 namedistr="hd0d0S_";
ea0d8716 292 namedistr+=i;
b272aebf 293 TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
294 namedistr="hd0d0B_";
ea0d8716 295 namedistr+=i;
b272aebf 296 TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
297
4e61a020 298 namedistr="hd0d0vS_";
299 namedistr+=i;
300 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);
301 namedistr="hd0d0vB_";
302 namedistr+=i;
303 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);
304
305 //decay lenght
306 namedistr="hdeclS_";
307 namedistr+=i;
308 TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.6);
309
310 namedistr="hdeclB_";
311 namedistr+=i;
312 TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.6);
313
314 namedistr="hnormdeclS_";
315 namedistr+=i;
316 TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,10.);
317
318 namedistr="hnormdeclB_";
319 namedistr+=i;
320 TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,10.);
321
322 // costhetapoint
b272aebf 323 namedistr="hcosthetapointS_";
ea0d8716 324 namedistr+=i;
b272aebf 325 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
326 namedistr="hcosthetapointB_";
ea0d8716 327 namedistr+=i;
b272aebf 328 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
329
4e61a020 330 // costhetapoint vs d0 or d0d0
331 namedistr="hcosthpointd0S_";
332 namedistr+=i;
333 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);
334
335 namedistr="hcosthpointd0B_";
336 namedistr+=i;
337 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);
338
b272aebf 339 namedistr="hcosthpointd0d0S_";
ea0d8716 340 namedistr+=i;
b272aebf 341 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);
342 namedistr="hcosthpointd0d0B_";
ea0d8716 343 namedistr+=i;
b272aebf 344 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);
345
346 fDistr->Add(hptpiS);
347 fDistr->Add(hptKS);
348 fDistr->Add(hptB);
349
350 fDistr->Add(hptpiSnoMcut);
351 fDistr->Add(hptKSnoMcut);
352 fDistr->Add(hptB1pnoMcut);
353 fDistr->Add(hptB2pnoMcut);
354
355 fDistr->Add(hdcaS);
356 fDistr->Add(hdcaB);
357
358 fDistr->Add(hd0piS);
359 fDistr->Add(hd0KS);
360 fDistr->Add(hd0B);
4e61a020 361 fDistr->Add(hd0p0B);
362 fDistr->Add(hd0p1B);
363
364 fDistr->Add(hd0vpiS);
365 fDistr->Add(hd0vKS);
366 fDistr->Add(hd0vB);
367 fDistr->Add(hd0vp0B);
368 fDistr->Add(hd0vp1B);
b272aebf 369
370 fDistr->Add(hd0d0S);
371 fDistr->Add(hd0d0B);
372
4e61a020 373 fDistr->Add(hd0d0vS);
374 fDistr->Add(hd0d0vB);
375
b272aebf 376 fDistr->Add(hcosthetastarS);
377 fDistr->Add(hcosthetastarB);
378
379 fDistr->Add(hcosthetapointS);
380 fDistr->Add(hcosthetapointB);
381
4e61a020 382 fDistr->Add(hdeclengthS);
383 fDistr->Add(hdeclengthB);
384
385 fDistr->Add(hnormdeclengthS);
386 fDistr->Add(hnormdeclengthB);
387
388 fDistr->Add(hcosthpointd0S);
389 fDistr->Add(hcosthpointd0B);
390
b272aebf 391 fDistr->Add(hcosthpointd0d0S);
392 fDistr->Add(hcosthpointd0d0B);
393
394
a41f6fad 395 //histograms of invariant mass distributions
396
46c96ce5 397 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
398 TH1F *tmpMl=(TH1F*)tmpMt->Clone();
399 tmpMt->Sumw2();
400 tmpMl->Sumw2();
a4ae02cd 401
9de8c723 402 //to compare with AliAnalysisTaskCharmFraction
403 TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.765,1.965);
404 TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
405 tmpS27t->Sumw2();
406 tmpS27l->Sumw2();
0108fa62 407
408 //distribution w/o cuts
6321ee46 409 // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
410 TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
411 TH1F *tmpMB=(TH1F*)tmpMS->Clone();
0108fa62 412 tmpMB->SetName(nameMassNocutsB.Data());
413 tmpMS->Sumw2();
414 tmpMB->Sumw2();
415
416 //MC signal and background
46c96ce5 417 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
418 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
419 tmpSt->Sumw2();
420 tmpSl->Sumw2();
a4ae02cd 421
46c96ce5 422 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
423 TH1F *tmpBl=(TH1F*)tmpBt->Clone();
424 tmpBt->Sumw2();
425 tmpBl->Sumw2();
a4ae02cd 426
feb73eca 427 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
46c96ce5 428 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
429 TH1F *tmpRl=(TH1F*)tmpRt->Clone();
430 tmpRt->Sumw2();
431 tmpRl->Sumw2();
a4ae02cd 432 // printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
433
ea0d8716 434 fOutputMass->Add(tmpMt);
435 fOutputMass->Add(tmpSt);
436 fOutputMass->Add(tmpS27t);
437 fOutputMass->Add(tmpBt);
438 fOutputMass->Add(tmpRt);
439
0108fa62 440 fDistr->Add(tmpMS);
441 fDistr->Add(tmpMB);
442
34137226 443
7646d6da 444 }
a4ae02cd 445
4e61a020 446 //histograms for vertex checking and TOF checking
34137226 447 TString checkname="hptGoodTr";
448
449 TH1F* hptGoodTr=new TH1F(checkname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Number",200,0.,8.);
450 hptGoodTr->SetTitleOffset(1.3,"Y");
451 checkname="hdistrGoodTr";
a41f6fad 452
34137226 453 TH1F* hdistrGoodTr=new TH1F(checkname.Data(),"Distribution of number of good tracks per event;no.good-tracks/ev;Entries",31,0,31);
454 hdistrGoodTr->SetTitleOffset(1.3,"Y");
455
4e61a020 456 checkname="hTOFsig";
457 TH1F* hTOFsig=new TH1F(checkname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
458
459 checkname="hTOFtimeKaonHyptime";
460 TH2F* hTOFtimeKaonHyptime=new TH2F(checkname.Data(),"TOFtime - timeHypothesisForKaon;p_{t}[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
461
462 checkname="hTOFtime";
463 TH1F* hTOFtime=new TH1F(checkname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
464
465
34137226 466 fChecks->Add(hptGoodTr);
467 fChecks->Add(hdistrGoodTr);
4e61a020 468 fChecks->Add(hTOFsig);
469 fChecks->Add(hTOFtimeKaonHyptime);
470 fChecks->Add(hTOFtime);
34137226 471
5b2e5fae 472 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
4e61a020 473
474 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", 8,0.,8.);
34137226 475
476 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
4e61a020 477 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
34137226 478 fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
479 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtx");
480 fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
4e61a020 481 fNentries->GetXaxis()->SetBinLabel(6,"ptbin = -1");
482 fNentries->GetXaxis()->SetBinLabel(7,"no daughter");
483 fNentries->GetXaxis()->SetBinLabel(8,"nCandSel(Tr)");
484
34137226 485 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
a4ae02cd 486
ea0d8716 487
40445ada 488 // Post the data
ea0d8716 489 PostData(1,fOutputMass);
5b2e5fae 490 PostData(2,fDistr);
40445ada 491 PostData(3,fNentries);
5b2e5fae 492 PostData(4,fChecks);
ea0d8716 493
49061176 494 return;
495}
496
497//________________________________________________________________________
498void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
499{
500 // Execute analysis for current event:
501 // heavy flavor candidates association to MC truth
a4ae02cd 502 //cout<<"I'm in UserExec"<<endl;
ea0d8716 503
504
505 //cuts order
506 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
507 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
508 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
509 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
510 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
511 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
512 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
513 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
514 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
515
516
49061176 517 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
feb73eca 518
b557eb43 519 TString bname;
feb73eca 520 if(fArray==0){ //D0 candidates
b557eb43 521 // load D0->Kpi candidates
feb73eca 522 //cout<<"D0 candidates"<<endl;
b557eb43 523 bname="D0toKpi";
feb73eca 524 } else { //LikeSign candidates
feb73eca 525 // load like sign candidates
b557eb43 526 //cout<<"LS candidates"<<endl;
527 bname="LikeSign2Prong";
528 }
529
530 TClonesArray *inputArray=0;
34137226 531
b557eb43 532 if(!aod && AODEvent() && IsStandardAOD()) {
533 // In case there is an AOD handler writing a standard AOD, use the AOD
534 // event in memory rather than the input (ESD) event.
535 aod = dynamic_cast<AliAODEvent*> (AODEvent());
536 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
537 // have to taken from the AOD event hold by the AliAODExtension
538 AliAODHandler* aodHandler = (AliAODHandler*)
539 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
34137226 540
b557eb43 541 if(aodHandler->GetExtensions()) {
542 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
34137226 543 AliAODEvent* aodFromExt = ext->GetAOD();
b557eb43 544 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
feb73eca 545 }
b557eb43 546 } else {
547 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
548 }
feb73eca 549
b557eb43 550
551 if(!inputArray) {
552 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
553 return;
49061176 554 }
555
34137226 556
ce39f0ac 557 TClonesArray *mcArray = 0;
558 AliAODMCHeader *mcHeader = 0;
559
560 if(fReadMC) {
561 // load MC particles
562 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
563 if(!mcArray) {
564 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
565 return;
566 }
40445ada 567
ce39f0ac 568 // load MC header
569 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
570 if(!mcHeader) {
571 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
572 return;
573 }
49061176 574 }
575
b557eb43 576 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
40445ada 577
a4ae02cd 578 //histogram filled with 1 for every AOD
34137226 579 fNentries->Fill(0);
49061176 580
0108fa62 581
40445ada 582
583 // AOD primary vertex
584 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
b272aebf 585
40445ada 586 Bool_t isGoodVtx=kFALSE;
b272aebf 587
40445ada 588 //vtx1->Print();
589 TString primTitle = vtx1->GetTitle();
590 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
591 isGoodVtx=kTRUE;
592 fNentries->Fill(3);
593 }
594
595 //cout<<"Start checks"<<endl;
596 Int_t ntracks=0,isGoodTrack=0;
597
598 if(aod) ntracks=aod->GetNTracks();
599
75638da0 600 //cout<<"ntracks = "<<ntracks<<endl;
ea0d8716 601
75638da0 602 //loop on tracks in the event
603 for (Int_t k=0;k<ntracks;k++){
604 AliAODTrack* track=aod->GetTrack(k);
4e61a020 605
606 //check TOF
607
608 if(!(track->GetStatus()&AliESDtrack::kTPCrefit &&
609 track->GetStatus()&AliESDtrack::kITSrefit &&
610 track->GetTPCNcls() >=70 &&
611 track->GetStatus()&AliESDtrack::kTOFpid &&
612 track->GetStatus()&AliESDtrack::kTOFout &&
613 track->GetStatus()&AliESDtrack::kTIME)) continue;
614 AliAODPid *pid = track->GetDetPid();
615 if(!pid) {cout<<"No AliAODPid found"<<endl; continue;}
616
617 Double_t times[5];
618 pid->GetIntegratedTimes(times);
619
620 ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
621 ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
622
623 ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
624 if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
625
626
75638da0 627 //check clusters of the tracks
628 Int_t nclsTot=0,nclsSPD=0;
40445ada 629
75638da0 630 for(Int_t l=0;l<6;l++) {
631 if(TESTBIT(track->GetITSClusterMap(),l)) {
632 nclsTot++; if(l<2) nclsSPD++;
b272aebf 633 }
75638da0 634 }
40445ada 635
75638da0 636 if (track->Pt()>0.3 &&
637 track->GetStatus()&AliESDtrack::kTPCrefit &&
638 track->GetStatus()&AliESDtrack::kITSrefit &&
639 nclsTot>3 &&
640 nclsSPD>0) {//fill hist good tracks
641 //cout<<"in if"<<endl;
642 ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
643 isGoodTrack++;
40445ada 644 }
75638da0 645 //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
646 ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
647 }
648 //number of events with good vertex and at least 2 good tracks
649 if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
ea0d8716 650
651
75638da0 652 // loop over candidates
653 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4e61a020 654 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
75638da0 655
75638da0 656 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
657 //Int_t nPosPairs=0, nNegPairs=0;
658 //cout<<"inside the loop"<<endl;
659 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
4e61a020 660
661 //check daughters
662 if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
663 AliDebug(1,"at least one daughter not found!");
664 fNentries->Fill(6);
665 return;
666 }
667
75638da0 668 Bool_t unsetvtx=kFALSE;
669 if(!d->GetOwnPrimaryVtx()) {
670 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
671 unsetvtx=kTRUE;
672 }
ea0d8716 673
4e61a020 674
75638da0 675 //check reco daughter in acceptance
4e61a020 676 /*
75638da0 677 Double_t eta0=d->EtaProng(0);
678 Double_t eta1=d->EtaProng(1);
4e61a020 679 */
680 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
681 //if( TMath::Abs(eta0)<0.9 && TMath::Abs(eta1)<0.9 ){
682 //apply cuts on tracks
3cc4604b 683 Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
4e61a020 684 if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70) isSelected=kFALSE;
3cc4604b 685 if (!isSelected) return;
4e61a020 686 fNentries->Fill(7);
687 if(fDebug>2) cout<<"tracks selected"<<endl;
3cc4604b 688
4e61a020 689 FillVarHists(aod,d,mcArray,fCuts,fDistr);
ea0d8716 690 FillMassHists(d,mcArray,fCuts,fOutputMass);
75638da0 691 }
ea0d8716 692
75638da0 693 if(unsetvtx) d->UnsetOwnPrimaryVtx();
694 } //end for prongs
695
ea0d8716 696
40445ada 697 // Post the data
ea0d8716 698 PostData(1,fOutputMass);
699 PostData(2,fDistr);
40445ada 700 PostData(3,fNentries);
ea0d8716 701 PostData(4,fChecks);
702
40445ada 703 return;
704}
b272aebf 705
40445ada 706//____________________________________________________________________________
4e61a020 707void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
40445ada 708 //
709 // function used in UserExec to fill variable histograms:
710 //
b272aebf 711
40445ada 712 //add distr here
713 UInt_t pdgs[2];
714
715 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
716 pdgs[0]=211;
717 pdgs[1]=321;
718 Double_t minvD0 = part->InvMassD0();
719 pdgs[1]=211;
720 pdgs[0]=321;
721 Double_t minvD0bar = part->InvMassD0bar();
722 //cout<<"inside mass cut"<<endl;
b272aebf 723
40445ada 724 Int_t pdgDgD0toKpi[2]={321,211};
725 Int_t lab=-9999;
726 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)
727 //Double_t pt = d->Pt(); //mother pt
728 Bool_t isSelected=kFALSE;
b272aebf 729
3cc4604b 730
ea0d8716 731 if(fCutOnDistr){
3cc4604b 732 isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate);
ea0d8716 733 if (!isSelected){
734 //cout<<"Not Selected"<<endl;
735 return;
736 }
737 }
0108fa62 738
4e61a020 739 Double_t invmasscut=0.03;
740
40445ada 741 TString fillthispi="",fillthisK="",fillthis="";
b272aebf 742
ea0d8716 743 Int_t ptbin=cuts->PtBin(part->Pt());
4e61a020 744 if(ptbin==-1) {fNentries->Fill(5); return;} //out of bounds
745
746 //recalculate vertex
747 AliAODVertex *vtxProp=0x0;
748 vtxProp=GetPrimaryVtxSkipped(aod,part);
749 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
750 dz1[0]=-99; dz2[0]=-99;
751 if(vtxProp) {
752 part->SetOwnPrimaryVtx(vtxProp);
753 //Bool_t unsetvtx=kTRUE;
754 //Calculate d0 for daughter tracks with Vtx Skipped
755 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)part->GetDaughter(0));
756 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)part->GetDaughter(1));
757 esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
758 esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
759 delete esdtrack1;
760 delete esdtrack2;
761 }
762
763 Double_t d0[2]={dz1[0],dz2[0]};
764
40445ada 765 if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
ea0d8716 766 //printf("\nif no cuts or cuts passed\n");
40445ada 767 if(lab>=0 && fReadMC){ //signal
0108fa62 768
40445ada 769 //check pdg of the prongs
770 AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
771 AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
4e61a020 772 if(!prong0 || !prong1) {
773 return;
774 }
40445ada 775 Int_t labprong[2];
776 labprong[0]=prong0->GetLabel();
777 labprong[1]=prong1->GetLabel();
778 AliAODMCParticle *mcprong=0;
779 Int_t pdgProng[2]={0,0};
4e61a020 780
40445ada 781 for (Int_t iprong=0;iprong<2;iprong++){
782 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
783 pdgProng[iprong]=mcprong->GetPdgCode();
0108fa62 784 }
785
40445ada 786 //no mass cut ditributions: ptbis
787
788 fillthispi="hptpiSnoMcut_";
789 fillthispi+=ptbin;
790
791 fillthisK="hptKSnoMcut_";
792 fillthisK+=ptbin;
793
794 if (TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321){
795 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
796 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
797 }else {
798 if (TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211){
799 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
800 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
801 }
802 }
803
804 //no mass cut ditributions: mass
805 fillthis="hMassS_";
806 fillthis+=ptbin;
807
808 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
809 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
810 }
811 else { //D0bar
812 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
813 }
4e61a020 814
0108fa62 815 //apply cut on invariant mass on the pair
4e61a020 816 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
40445ada 817
0108fa62 818 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
819 for (Int_t iprong=0; iprong<2; iprong++){
40445ada 820 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
04e09ffc 821 labprong[iprong]=prong->GetLabel();
40445ada 822
0108fa62 823 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
40445ada 824 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
0108fa62 825 Int_t pdgprong=mcprong->GetPdgCode();
40445ada 826
0108fa62 827 if(TMath::Abs(pdgprong)==211) {
828 //cout<<"pi"<<endl;
40445ada 829
830 fillthispi="hptpiS_";
831 fillthispi+=ptbin;
832 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
4e61a020 833 fillthispi="hd0piS_";
834 fillthispi+=ptbin;
835 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
836 if(d0[iprong] != -99) {
837
838 fillthispi="hd0vpiS_";
839 fillthispi+=ptbin;
840 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
841 }
842
b272aebf 843 }
40445ada 844
845 if(TMath::Abs(pdgprong)==321) {
846 //cout<<"kappa"<<endl;
847
848 fillthisK="hptKS_";
849 fillthisK+=ptbin;
850 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
40445ada 851 fillthisK="hd0KS_";
852 fillthisK+=ptbin;
853 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
4e61a020 854 if (d0[iprong] != -99){
855 fillthisK="hd0vKS_";
856 fillthisK+=ptbin;
857 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
858 }
b272aebf 859 }
b272aebf 860
4e61a020 861 fillthis="hcosthpointd0S_";
40445ada 862 fillthis+=ptbin;
4e61a020 863 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
b272aebf 864
4e61a020 865 } //end loop on prongs
b272aebf 866
4e61a020 867 fillthis="hdcaS_";
868 fillthis+=ptbin;
869 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
870
871 fillthis="hcosthetapointS_";
872 fillthis+=ptbin;
873 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
874
875 fillthis="hcosthpointd0d0S_";
876 fillthis+=ptbin;
877 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
b272aebf 878
4e61a020 879 fillthis="hcosthetastarS_";
880 fillthis+=ptbin;
881 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
882 else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
0108fa62 883
4e61a020 884 fillthis="hd0d0S_";
885 fillthis+=ptbin;
886 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
887
888 if(d0[0] != -99 && d0[1] != -99 ){
889 fillthis="hd0d0vS_";
40445ada 890 fillthis+=ptbin;
4e61a020 891 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
b272aebf 892 }
4e61a020 893
894 fillthis="hdeclS_";
895 fillthis+=ptbin;
896 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
897
898 fillthis="hnormdeclS_";
899 fillthis+=ptbin;
900 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
901
902 } //end mass cut
903
0108fa62 904 } else{ //Background or LS
905 //cout<<"is background"<<endl;
40445ada 906
b272aebf 907 //no mass cut distributions: mass, ptbis
40445ada 908 fillthis="hMassB_";
909 fillthis+=ptbin;
ea0d8716 910
911 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
912 if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
b272aebf 913
40445ada 914 fillthis="hptB1prongnoMcut_";
915 fillthis+=ptbin;
916
917 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
918
919 fillthis="hptB2prongsnoMcut_";
920 fillthis+=ptbin;
921 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
922 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 923
924 //apply cut on invariant mass on the pair
4e61a020 925 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
0108fa62 926
927
40445ada 928 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 929 if(!prong) {
930 if(fDebug>2) cout<<"No daughter found";
931 return;
932 }
0108fa62 933 else{
99012cda 934 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
0108fa62 935 }
40445ada 936
937 //normalise pt distr to half afterwards
938 fillthis="hptB_";
939 fillthis+=ptbin;
4e61a020 940 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
941 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 942
4e61a020 943 fillthis="hd0p0B_";
40445ada 944 fillthis+=ptbin;
945 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
4e61a020 946 fillthis="hd0p1B_";
947 fillthis+=ptbin;
948 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
949
950 fillthis="hd0vp0B_";
951 fillthis+=ptbin;
952 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
953 fillthis="hd0vp1B_";
954 fillthis+=ptbin;
955 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
956
957 fillthis="hcosthpointd0B_";
958 fillthis+=ptbin;
959 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
960 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
961
ce39f0ac 962
40445ada 963 fillthis="hdcaB_";
964 fillthis+=ptbin;
965 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
34137226 966
40445ada 967 fillthis="hcosthetastarB_";
968 fillthis+=ptbin;
ea0d8716 969 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
4e61a020 970 if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
40445ada 971
972 fillthis="hd0d0B_";
973 fillthis+=ptbin;
974 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
975
976 fillthis="hcosthetapointB_";
977 fillthis+=ptbin;
978 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
979
980 fillthis="hcosthpointd0d0B_";
981 fillthis+=ptbin;
982 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
983
4e61a020 984 }//mass cut
985 }//else (background)
40445ada 986 }
49061176 987 return;
988}
ea0d8716 989//____________________________________________________________________________
990
991void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
49061176 992 //
40445ada 993 // function used in UserExec to fill mass histograms:
49061176 994 //
9de8c723 995
996
997 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
998
3cc4604b 999 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
ea0d8716 1000 //cout<<"is selected = "<<isSelected<<endl;
feb73eca 1001
ea0d8716 1002 //cout<<"check cuts = "<<endl;
1003 //cuts->PrintAll();
1004 if (!isSelected){
1005 //cout<<"Not Selected"<<endl;
1006 return;
1007 }
9de8c723 1008
4e61a020 1009 if(fDebug>2) cout<<"Candidate selected"<<endl;
a41f6fad 1010
ea0d8716 1011 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1012 //printf("SELECTED\n");
1013 Int_t ptbin=cuts->PtBin(part->Pt());
9de8c723 1014
ea0d8716 1015 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 1016 if(!prong) {
1017 AliDebug(2,"No daughter found");
1018 return;
1019 }
ea0d8716 1020 else{
4b4a1d25 1021 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
ea0d8716 1022 }
7646d6da 1023
ea0d8716 1024 TString fillthis="";
1025 Int_t pdgDgD0toKpi[2]={321,211};
1026 Int_t labD0=-1;
1027 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)
1028
1029 //count candidates selected by cuts
1030 fNentries->Fill(1);
1031 //count true D0 selected by cuts
1032 if (fReadMC && labD0>=0) fNentries->Fill(2);
1033 //PostData(3,fNentries);
1034
1035 if (isSelected==1 || isSelected==3) { //D0
1036 fillthis="histMass_";
1037 fillthis+=ptbin;
1038 //cout<<"Filling "<<fillthis<<endl;
1039
1040 //printf("Fill mass with D0");
1041 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
feb73eca 1042
ea0d8716 1043 if(labD0>=0) {
1044 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1045
1046 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1047 Int_t pdgD0 = partD0->GetPdgCode();
1048 //cout<<"pdg = "<<pdgD0<<endl;
1049 if (pdgD0==421){ //D0
1050 //cout<<"Fill S with D0"<<endl;
1051 fillthis="histSgn_";
1052 fillthis+=ptbin;
1053 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1054 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1055 fillthis="histSgn27_";
46c96ce5 1056 fillthis+=ptbin;
1057 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
a4ae02cd 1058 }
ea0d8716 1059 } else{ //it was a D0bar
1060 fillthis="histRfl_";
a4ae02cd 1061 fillthis+=ptbin;
1062 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1063 }
ea0d8716 1064 } else {//background
1065 fillthis="histBkg_";
a4ae02cd 1066 fillthis+=ptbin;
ea0d8716 1067 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1068 }
1069
1070 }
1071 if (isSelected>1) { //D0bar
1072 fillthis="histMass_";
1073 fillthis+=ptbin;
1074 //printf("Fill mass with D0bar");
1075 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
feb73eca 1076
ea0d8716 1077 if(labD0>=0) {
1078 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1079 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1080 Int_t pdgD0 = partD0->GetPdgCode();
1081 //cout<<" pdg = "<<pdgD0<<endl;
1082 if (pdgD0==-421){ //D0bar
1083 fillthis="histSgn_";
1084 fillthis+=ptbin;
1085 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1086 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1087 fillthis="histSgn27_";
a4ae02cd 1088 fillthis+=ptbin;
1089 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
ea0d8716 1090 }
9de8c723 1091
a4ae02cd 1092
ea0d8716 1093 } else{
1094 fillthis="histRfl_";
a4ae02cd 1095 fillthis+=ptbin;
1096 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1097 }
ea0d8716 1098 } else {//background or LS
1099 fillthis="histBkg_";
1100 fillthis+=ptbin;
1101 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
a4ae02cd 1102 }
ea0d8716 1103 }
a4ae02cd 1104
40445ada 1105 return;
49061176 1106}
4e61a020 1107
1108//__________________________________________________________________________
1109AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1110 //Calculate the primary vertex w/o the daughter tracks of the candidate
1111
1112 AliESDVertex *vertexESD=0x0;
1113 AliAODVertex *vertexAOD=0x0;
1114 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1115
1116 Int_t skipped[2];
1117 Int_t nTrksToSkip=2;
1118 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1119 if(!dgTrack){
1120 AliDebug(2,"no daughter found!");
1121 return 0x0;
1122 }
1123 skipped[0]=dgTrack->GetID();
1124 dgTrack = (AliAODTrack*)d->GetDaughter(1);
1125 if(!dgTrack){
1126 AliDebug(2,"no daughter found!");
1127 return 0x0;
1128 }
1129 skipped[1]=dgTrack->GetID();
1130
1131
1132 //
1133 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1134 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
1135 vertexer->SetMinClusters(4);
1136 if(!vertexESD) return vertexAOD;
1137 if(vertexESD->GetNContributors()<=0) {
1138 AliDebug(2,"vertexing failed");
1139 delete vertexESD; vertexESD=NULL;
1140 return vertexAOD;
1141 }
1142
1143 delete vertexer; vertexer=NULL;
1144
1145
1146 // convert to AliAODVertex
1147 Double_t pos[3],cov[6],chi2perNDF;
1148 vertexESD->GetXYZ(pos); // position
1149 vertexESD->GetCovMatrix(cov); //covariance matrix
1150 chi2perNDF = vertexESD->GetChi2toNDF();
1151 delete vertexESD; vertexESD=NULL;
1152
1153 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1154 return vertexAOD;
1155
1156}
1157
1158
49061176 1159//________________________________________________________________________
1160void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1161{
1162 // Terminate analysis
1163 //
1164 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1165
6321ee46 1166
ea0d8716 1167 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1168 if (!fOutputMass) {
1169 printf("ERROR: fOutputMass not available\n");
a4ae02cd 1170 return;
1171 }
ea0d8716 1172 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1173 if (!fDistr) {
1174 printf("ERROR: fDistr not available\n");
a41f6fad 1175 return;
1176 }
5b2e5fae 1177
40445ada 1178 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
5b2e5fae 1179
40445ada 1180 if(!fNentries){
1181 printf("ERROR: fNEntries not available\n");
1182 return;
1183 }
5b2e5fae 1184
ea0d8716 1185 fChecks = dynamic_cast<TList*> (GetOutputData(4));
34137226 1186 if (!fChecks) {
1187 printf("ERROR: fChecks not available\n");
1188 return;
1189 }
ea0d8716 1190
4e61a020 1191
1192 for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1193
1194 //sum d0p0 and d0p1
1195 TString d0Bname=Form("hd0B_%d",ipt),d0p0Bname=Form("hd0p0B_%d",ipt),d0p1Bname=Form("hd0p1B_%d",ipt);
1196 ((TH1F*)fDistr->FindObject(d0Bname))->Add( ((TH1F*)fDistr->FindObject(d0p0Bname)), ((TH1F*)fDistr->FindObject(d0p1Bname)) );
1197
1198
1199 if(fArray==1){
feb73eca 1200 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1201
1202
ea0d8716 1203 if(fLsNormalization>1e-6) {
9de8c723 1204
feb73eca 1205 TString massName="histMass_";
ea0d8716 1206 massName+=ipt;
1207 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1208
feb73eca 1209 }
40445ada 1210
feb73eca 1211
40445ada 1212 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1213
ea0d8716 1214 if(fLsNormalization>1e-6) {
40445ada 1215
1216 TString nameDistr="hptB_";
ea0d8716 1217 nameDistr+=ipt;
40445ada 1218 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1219 nameDistr="hdcaB_";
ea0d8716 1220 nameDistr+=ipt;
40445ada 1221 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1222 nameDistr="hcosthetastarB_";
ea0d8716 1223 nameDistr+=ipt;
40445ada 1224 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1225 nameDistr="hd0B_";
ea0d8716 1226 nameDistr+=ipt;
40445ada 1227 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1228 nameDistr="hd0d0B_";
ea0d8716 1229 nameDistr+=ipt;
40445ada 1230 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1231 nameDistr="hcosthetapointB_";
ea0d8716 1232 nameDistr+=ipt;
40445ada 1233 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1234 nameDistr="hcosthpointd0d0B_";
ea0d8716 1235 nameDistr+=ipt;
40445ada 1236 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
feb73eca 1237
40445ada 1238 }
feb73eca 1239 }
1240 }
4e61a020 1241 TString cvname,cstname;
527f330b 1242
1243 if (fArray==0){
1244 cvname="D0invmass";
4e61a020 1245 cstname="cstat0";
1246 } else {
1247 cvname="LSinvmass";
1248 cstname="cstat1";
1249 }
527f330b 1250
34137226 1251 TCanvas *cMass=new TCanvas(cvname,cvname);
1252 cMass->cd();
ea0d8716 1253 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
527f330b 1254
4e61a020 1255 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
34137226 1256 cStat->cd();
1257 cStat->SetGridy();
1258 fNentries->Draw("htext0");
527f330b 1259
4e61a020 1260 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1261 // ccheck->cd();
1262
49061176 1263 return;
1264}
1265