]> 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
fixed bugs with the rule checker.
[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
a2121012 322 namedistr="hdeclvS_";
323 namedistr+=i;
324 TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
325
326 namedistr="hdeclvB_";
327 namedistr+=i;
328 TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
329
330 namedistr="hnormdeclvS_";
331 namedistr+=i;
332 TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
333
334 namedistr="hnormdeclvB_";
335 namedistr+=i;
336 TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
337
4e61a020 338 // costhetapoint
b272aebf 339 namedistr="hcosthetapointS_";
ea0d8716 340 namedistr+=i;
b272aebf 341 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
342 namedistr="hcosthetapointB_";
ea0d8716 343 namedistr+=i;
b272aebf 344 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
345
4e61a020 346 // costhetapoint vs d0 or d0d0
347 namedistr="hcosthpointd0S_";
348 namedistr+=i;
349 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);
350
351 namedistr="hcosthpointd0B_";
352 namedistr+=i;
353 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);
354
b272aebf 355 namedistr="hcosthpointd0d0S_";
ea0d8716 356 namedistr+=i;
b272aebf 357 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);
358 namedistr="hcosthpointd0d0B_";
ea0d8716 359 namedistr+=i;
b272aebf 360 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);
361
362 fDistr->Add(hptpiS);
363 fDistr->Add(hptKS);
364 fDistr->Add(hptB);
365
366 fDistr->Add(hptpiSnoMcut);
367 fDistr->Add(hptKSnoMcut);
368 fDistr->Add(hptB1pnoMcut);
369 fDistr->Add(hptB2pnoMcut);
370
371 fDistr->Add(hdcaS);
372 fDistr->Add(hdcaB);
373
374 fDistr->Add(hd0piS);
375 fDistr->Add(hd0KS);
376 fDistr->Add(hd0B);
4e61a020 377 fDistr->Add(hd0p0B);
378 fDistr->Add(hd0p1B);
379
380 fDistr->Add(hd0vpiS);
381 fDistr->Add(hd0vKS);
382 fDistr->Add(hd0vB);
383 fDistr->Add(hd0vp0B);
384 fDistr->Add(hd0vp1B);
b272aebf 385
386 fDistr->Add(hd0d0S);
387 fDistr->Add(hd0d0B);
388
4e61a020 389 fDistr->Add(hd0d0vS);
390 fDistr->Add(hd0d0vB);
391
b272aebf 392 fDistr->Add(hcosthetastarS);
393 fDistr->Add(hcosthetastarB);
394
395 fDistr->Add(hcosthetapointS);
396 fDistr->Add(hcosthetapointB);
397
4e61a020 398 fDistr->Add(hdeclengthS);
399 fDistr->Add(hdeclengthB);
400
401 fDistr->Add(hnormdeclengthS);
402 fDistr->Add(hnormdeclengthB);
403
a2121012 404 fDistr->Add(hdeclengthvS);
405 fDistr->Add(hdeclengthvB);
406
407 fDistr->Add(hnormdeclengthvS);
408 fDistr->Add(hnormdeclengthvB);
409
4e61a020 410 fDistr->Add(hcosthpointd0S);
411 fDistr->Add(hcosthpointd0B);
412
b272aebf 413 fDistr->Add(hcosthpointd0d0S);
414 fDistr->Add(hcosthpointd0d0B);
415
416
a41f6fad 417 //histograms of invariant mass distributions
418
46c96ce5 419 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
420 TH1F *tmpMl=(TH1F*)tmpMt->Clone();
421 tmpMt->Sumw2();
422 tmpMl->Sumw2();
a4ae02cd 423
9de8c723 424 //to compare with AliAnalysisTaskCharmFraction
425 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);
426 TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
427 tmpS27t->Sumw2();
428 tmpS27l->Sumw2();
0108fa62 429
430 //distribution w/o cuts
6321ee46 431 // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
432 TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
433 TH1F *tmpMB=(TH1F*)tmpMS->Clone();
0108fa62 434 tmpMB->SetName(nameMassNocutsB.Data());
435 tmpMS->Sumw2();
436 tmpMB->Sumw2();
437
438 //MC signal and background
46c96ce5 439 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
440 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
441 tmpSt->Sumw2();
442 tmpSl->Sumw2();
a4ae02cd 443
46c96ce5 444 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
445 TH1F *tmpBl=(TH1F*)tmpBt->Clone();
446 tmpBt->Sumw2();
447 tmpBl->Sumw2();
a4ae02cd 448
feb73eca 449 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
46c96ce5 450 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
451 TH1F *tmpRl=(TH1F*)tmpRt->Clone();
452 tmpRt->Sumw2();
453 tmpRl->Sumw2();
a4ae02cd 454 // printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
455
ea0d8716 456 fOutputMass->Add(tmpMt);
457 fOutputMass->Add(tmpSt);
458 fOutputMass->Add(tmpS27t);
459 fOutputMass->Add(tmpBt);
460 fOutputMass->Add(tmpRt);
461
0108fa62 462 fDistr->Add(tmpMS);
463 fDistr->Add(tmpMB);
464
34137226 465
7646d6da 466 }
a4ae02cd 467
4e61a020 468 //histograms for vertex checking and TOF checking
34137226 469 TString checkname="hptGoodTr";
470
471 TH1F* hptGoodTr=new TH1F(checkname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Number",200,0.,8.);
472 hptGoodTr->SetTitleOffset(1.3,"Y");
473 checkname="hdistrGoodTr";
a41f6fad 474
34137226 475 TH1F* hdistrGoodTr=new TH1F(checkname.Data(),"Distribution of number of good tracks per event;no.good-tracks/ev;Entries",31,0,31);
476 hdistrGoodTr->SetTitleOffset(1.3,"Y");
477
4e61a020 478 checkname="hTOFsig";
479 TH1F* hTOFsig=new TH1F(checkname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
480
481 checkname="hTOFtimeKaonHyptime";
482 TH2F* hTOFtimeKaonHyptime=new TH2F(checkname.Data(),"TOFtime - timeHypothesisForKaon;p_{t}[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
483
484 checkname="hTOFtime";
485 TH1F* hTOFtime=new TH1F(checkname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
486
487
34137226 488 fChecks->Add(hptGoodTr);
489 fChecks->Add(hdistrGoodTr);
4e61a020 490 fChecks->Add(hTOFsig);
491 fChecks->Add(hTOFtimeKaonHyptime);
492 fChecks->Add(hTOFtime);
34137226 493
5b2e5fae 494 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
4e61a020 495
496 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 497
498 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
4e61a020 499 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
34137226 500 fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
501 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtx");
502 fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
4e61a020 503 fNentries->GetXaxis()->SetBinLabel(6,"ptbin = -1");
504 fNentries->GetXaxis()->SetBinLabel(7,"no daughter");
505 fNentries->GetXaxis()->SetBinLabel(8,"nCandSel(Tr)");
506
34137226 507 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
a4ae02cd 508
ea0d8716 509
40445ada 510 // Post the data
ea0d8716 511 PostData(1,fOutputMass);
5b2e5fae 512 PostData(2,fDistr);
40445ada 513 PostData(3,fNentries);
5b2e5fae 514 PostData(4,fChecks);
ea0d8716 515
49061176 516 return;
517}
518
519//________________________________________________________________________
520void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
521{
522 // Execute analysis for current event:
523 // heavy flavor candidates association to MC truth
a4ae02cd 524 //cout<<"I'm in UserExec"<<endl;
ea0d8716 525
526
527 //cuts order
528 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
529 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
530 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
531 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
532 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
533 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
534 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
535 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
536 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
537
538
49061176 539 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
feb73eca 540
b557eb43 541 TString bname;
feb73eca 542 if(fArray==0){ //D0 candidates
b557eb43 543 // load D0->Kpi candidates
feb73eca 544 //cout<<"D0 candidates"<<endl;
b557eb43 545 bname="D0toKpi";
feb73eca 546 } else { //LikeSign candidates
feb73eca 547 // load like sign candidates
b557eb43 548 //cout<<"LS candidates"<<endl;
549 bname="LikeSign2Prong";
550 }
551
552 TClonesArray *inputArray=0;
34137226 553
b557eb43 554 if(!aod && AODEvent() && IsStandardAOD()) {
555 // In case there is an AOD handler writing a standard AOD, use the AOD
556 // event in memory rather than the input (ESD) event.
557 aod = dynamic_cast<AliAODEvent*> (AODEvent());
558 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
559 // have to taken from the AOD event hold by the AliAODExtension
560 AliAODHandler* aodHandler = (AliAODHandler*)
561 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
34137226 562
b557eb43 563 if(aodHandler->GetExtensions()) {
564 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
34137226 565 AliAODEvent* aodFromExt = ext->GetAOD();
b557eb43 566 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
feb73eca 567 }
b557eb43 568 } else {
569 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
570 }
feb73eca 571
b557eb43 572
573 if(!inputArray) {
574 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
575 return;
49061176 576 }
577
34137226 578
ce39f0ac 579 TClonesArray *mcArray = 0;
580 AliAODMCHeader *mcHeader = 0;
581
582 if(fReadMC) {
583 // load MC particles
584 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
585 if(!mcArray) {
586 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
587 return;
588 }
40445ada 589
ce39f0ac 590 // load MC header
591 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
592 if(!mcHeader) {
593 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
594 return;
595 }
49061176 596 }
597
b557eb43 598 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
40445ada 599
a4ae02cd 600 //histogram filled with 1 for every AOD
34137226 601 fNentries->Fill(0);
49061176 602
0108fa62 603
40445ada 604
605 // AOD primary vertex
606 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
c59e1898 607 if(!vtx1) {
608 AliError("There is no primary vertex !");
609 return;
610 }
b272aebf 611
40445ada 612 Bool_t isGoodVtx=kFALSE;
b272aebf 613
40445ada 614 //vtx1->Print();
615 TString primTitle = vtx1->GetTitle();
616 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
617 isGoodVtx=kTRUE;
618 fNentries->Fill(3);
619 }
620
621 //cout<<"Start checks"<<endl;
622 Int_t ntracks=0,isGoodTrack=0;
623
624 if(aod) ntracks=aod->GetNTracks();
625
75638da0 626 //cout<<"ntracks = "<<ntracks<<endl;
ea0d8716 627
75638da0 628 //loop on tracks in the event
629 for (Int_t k=0;k<ntracks;k++){
630 AliAODTrack* track=aod->GetTrack(k);
4e61a020 631
632 //check TOF
633
634 if(!(track->GetStatus()&AliESDtrack::kTPCrefit &&
635 track->GetStatus()&AliESDtrack::kITSrefit &&
636 track->GetTPCNcls() >=70 &&
637 track->GetStatus()&AliESDtrack::kTOFpid &&
638 track->GetStatus()&AliESDtrack::kTOFout &&
639 track->GetStatus()&AliESDtrack::kTIME)) continue;
640 AliAODPid *pid = track->GetDetPid();
a2121012 641 if(!pid) {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
4e61a020 642
643 Double_t times[5];
644 pid->GetIntegratedTimes(times);
645
646 ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
647 ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
648
649 ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
650 if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
651
652
75638da0 653 //check clusters of the tracks
654 Int_t nclsTot=0,nclsSPD=0;
40445ada 655
75638da0 656 for(Int_t l=0;l<6;l++) {
657 if(TESTBIT(track->GetITSClusterMap(),l)) {
658 nclsTot++; if(l<2) nclsSPD++;
b272aebf 659 }
75638da0 660 }
40445ada 661
75638da0 662 if (track->Pt()>0.3 &&
663 track->GetStatus()&AliESDtrack::kTPCrefit &&
664 track->GetStatus()&AliESDtrack::kITSrefit &&
665 nclsTot>3 &&
666 nclsSPD>0) {//fill hist good tracks
a2121012 667
75638da0 668 ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
669 isGoodTrack++;
40445ada 670 }
75638da0 671 //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
672 ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
673 }
674 //number of events with good vertex and at least 2 good tracks
675 if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
ea0d8716 676
677
75638da0 678 // loop over candidates
679 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4e61a020 680 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
75638da0 681
75638da0 682 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
683 //Int_t nPosPairs=0, nNegPairs=0;
684 //cout<<"inside the loop"<<endl;
685 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
4e61a020 686
687 //check daughters
688 if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
689 AliDebug(1,"at least one daughter not found!");
690 fNentries->Fill(6);
691 return;
692 }
693
75638da0 694 Bool_t unsetvtx=kFALSE;
695 if(!d->GetOwnPrimaryVtx()) {
696 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
697 unsetvtx=kTRUE;
698 }
ea0d8716 699
4e61a020 700
75638da0 701 //check reco daughter in acceptance
4e61a020 702 /*
75638da0 703 Double_t eta0=d->EtaProng(0);
704 Double_t eta1=d->EtaProng(1);
4e61a020 705 */
706 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
707 //if( TMath::Abs(eta0)<0.9 && TMath::Abs(eta1)<0.9 ){
708 //apply cuts on tracks
3cc4604b 709 Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
a2121012 710 if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(1))->GetTPCNcls() < 70) isSelected=kFALSE;
3cc4604b 711 if (!isSelected) return;
4e61a020 712 fNentries->Fill(7);
713 if(fDebug>2) cout<<"tracks selected"<<endl;
3cc4604b 714
4e61a020 715 FillVarHists(aod,d,mcArray,fCuts,fDistr);
ea0d8716 716 FillMassHists(d,mcArray,fCuts,fOutputMass);
75638da0 717 }
ea0d8716 718
75638da0 719 if(unsetvtx) d->UnsetOwnPrimaryVtx();
720 } //end for prongs
721
ea0d8716 722
40445ada 723 // Post the data
ea0d8716 724 PostData(1,fOutputMass);
725 PostData(2,fDistr);
40445ada 726 PostData(3,fNentries);
ea0d8716 727 PostData(4,fChecks);
728
40445ada 729 return;
730}
b272aebf 731
40445ada 732//____________________________________________________________________________
4e61a020 733void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
40445ada 734 //
735 // function used in UserExec to fill variable histograms:
736 //
b272aebf 737
40445ada 738 //add distr here
739 UInt_t pdgs[2];
740
741 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
742 pdgs[0]=211;
743 pdgs[1]=321;
744 Double_t minvD0 = part->InvMassD0();
745 pdgs[1]=211;
746 pdgs[0]=321;
747 Double_t minvD0bar = part->InvMassD0bar();
748 //cout<<"inside mass cut"<<endl;
b272aebf 749
40445ada 750 Int_t pdgDgD0toKpi[2]={321,211};
751 Int_t lab=-9999;
752 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)
753 //Double_t pt = d->Pt(); //mother pt
754 Bool_t isSelected=kFALSE;
b272aebf 755
3cc4604b 756
ea0d8716 757 if(fCutOnDistr){
3cc4604b 758 isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate);
ea0d8716 759 if (!isSelected){
760 //cout<<"Not Selected"<<endl;
761 return;
762 }
763 }
0108fa62 764
4e61a020 765 Double_t invmasscut=0.03;
766
40445ada 767 TString fillthispi="",fillthisK="",fillthis="";
b272aebf 768
ea0d8716 769 Int_t ptbin=cuts->PtBin(part->Pt());
4e61a020 770 if(ptbin==-1) {fNentries->Fill(5); return;} //out of bounds
771
772 //recalculate vertex
773 AliAODVertex *vtxProp=0x0;
774 vtxProp=GetPrimaryVtxSkipped(aod,part);
775 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
776 dz1[0]=-99; dz2[0]=-99;
777 if(vtxProp) {
778 part->SetOwnPrimaryVtx(vtxProp);
779 //Bool_t unsetvtx=kTRUE;
780 //Calculate d0 for daughter tracks with Vtx Skipped
781 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)part->GetDaughter(0));
782 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)part->GetDaughter(1));
783 esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
784 esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
785 delete esdtrack1;
786 delete esdtrack2;
787 }
788
789 Double_t d0[2]={dz1[0],dz2[0]};
790
40445ada 791 if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
ea0d8716 792 //printf("\nif no cuts or cuts passed\n");
40445ada 793 if(lab>=0 && fReadMC){ //signal
0108fa62 794
40445ada 795 //check pdg of the prongs
796 AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
797 AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
4e61a020 798 if(!prong0 || !prong1) {
799 return;
800 }
40445ada 801 Int_t labprong[2];
802 labprong[0]=prong0->GetLabel();
803 labprong[1]=prong1->GetLabel();
804 AliAODMCParticle *mcprong=0;
805 Int_t pdgProng[2]={0,0};
4e61a020 806
40445ada 807 for (Int_t iprong=0;iprong<2;iprong++){
808 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
809 pdgProng[iprong]=mcprong->GetPdgCode();
0108fa62 810 }
811
40445ada 812 //no mass cut ditributions: ptbis
813
814 fillthispi="hptpiSnoMcut_";
815 fillthispi+=ptbin;
816
817 fillthisK="hptKSnoMcut_";
818 fillthisK+=ptbin;
819
820 if (TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321){
821 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
822 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
823 }else {
824 if (TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211){
825 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
826 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
827 }
828 }
829
830 //no mass cut ditributions: mass
831 fillthis="hMassS_";
832 fillthis+=ptbin;
833
834 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
835 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
836 }
837 else { //D0bar
838 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
839 }
4e61a020 840
0108fa62 841 //apply cut on invariant mass on the pair
4e61a020 842 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
40445ada 843
0108fa62 844 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
845 for (Int_t iprong=0; iprong<2; iprong++){
40445ada 846 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
04e09ffc 847 labprong[iprong]=prong->GetLabel();
40445ada 848
0108fa62 849 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
40445ada 850 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
0108fa62 851 Int_t pdgprong=mcprong->GetPdgCode();
40445ada 852
0108fa62 853 if(TMath::Abs(pdgprong)==211) {
854 //cout<<"pi"<<endl;
40445ada 855
856 fillthispi="hptpiS_";
857 fillthispi+=ptbin;
858 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
4e61a020 859 fillthispi="hd0piS_";
860 fillthispi+=ptbin;
861 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
862 if(d0[iprong] != -99) {
863
864 fillthispi="hd0vpiS_";
865 fillthispi+=ptbin;
866 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
867 }
868
b272aebf 869 }
40445ada 870
871 if(TMath::Abs(pdgprong)==321) {
872 //cout<<"kappa"<<endl;
873
874 fillthisK="hptKS_";
875 fillthisK+=ptbin;
876 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
40445ada 877 fillthisK="hd0KS_";
878 fillthisK+=ptbin;
879 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
4e61a020 880 if (d0[iprong] != -99){
881 fillthisK="hd0vKS_";
882 fillthisK+=ptbin;
883 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
884 }
b272aebf 885 }
b272aebf 886
4e61a020 887 fillthis="hcosthpointd0S_";
40445ada 888 fillthis+=ptbin;
4e61a020 889 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
b272aebf 890
4e61a020 891 } //end loop on prongs
b272aebf 892
4e61a020 893 fillthis="hdcaS_";
894 fillthis+=ptbin;
895 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
896
897 fillthis="hcosthetapointS_";
898 fillthis+=ptbin;
899 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
900
901 fillthis="hcosthpointd0d0S_";
902 fillthis+=ptbin;
903 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
b272aebf 904
4e61a020 905 fillthis="hcosthetastarS_";
906 fillthis+=ptbin;
907 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
908 else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
0108fa62 909
4e61a020 910 fillthis="hd0d0S_";
911 fillthis+=ptbin;
912 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
913
914 if(d0[0] != -99 && d0[1] != -99 ){
915 fillthis="hd0d0vS_";
40445ada 916 fillthis+=ptbin;
4e61a020 917 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
b272aebf 918 }
4e61a020 919
920 fillthis="hdeclS_";
921 fillthis+=ptbin;
922 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
923
924 fillthis="hnormdeclS_";
925 fillthis+=ptbin;
926 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
a2121012 927
928 if(vtxProp) {
929 fillthis="hdeclvS_";
930 fillthis+=ptbin;
931 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
932
933 fillthis="hnormdeclvS_";
934 fillthis+=ptbin;
935 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
936 }
937 } //end mass cut
4e61a020 938
0108fa62 939 } else{ //Background or LS
940 //cout<<"is background"<<endl;
40445ada 941
b272aebf 942 //no mass cut distributions: mass, ptbis
40445ada 943 fillthis="hMassB_";
944 fillthis+=ptbin;
ea0d8716 945
946 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
947 if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
b272aebf 948
40445ada 949 fillthis="hptB1prongnoMcut_";
950 fillthis+=ptbin;
951
952 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
953
954 fillthis="hptB2prongsnoMcut_";
955 fillthis+=ptbin;
956 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
957 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 958
959 //apply cut on invariant mass on the pair
4e61a020 960 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
0108fa62 961
962
40445ada 963 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 964 if(!prong) {
965 if(fDebug>2) cout<<"No daughter found";
966 return;
967 }
0108fa62 968 else{
99012cda 969 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
0108fa62 970 }
40445ada 971
972 //normalise pt distr to half afterwards
973 fillthis="hptB_";
974 fillthis+=ptbin;
4e61a020 975 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
976 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 977
4e61a020 978 fillthis="hd0p0B_";
40445ada 979 fillthis+=ptbin;
980 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
4e61a020 981 fillthis="hd0p1B_";
982 fillthis+=ptbin;
983 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
a2121012 984 fillthis="hd0B_";
985 fillthis+=ptbin;
986 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
987 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
4e61a020 988
989 fillthis="hd0vp0B_";
990 fillthis+=ptbin;
991 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
992 fillthis="hd0vp1B_";
993 fillthis+=ptbin;
994 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
995
a2121012 996 fillthis="hd0vB_";
997 fillthis+=ptbin;
998 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
999 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1000
1001
4e61a020 1002 fillthis="hcosthpointd0B_";
1003 fillthis+=ptbin;
1004 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
1005 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
1006
ce39f0ac 1007
40445ada 1008 fillthis="hdcaB_";
1009 fillthis+=ptbin;
1010 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
34137226 1011
40445ada 1012 fillthis="hcosthetastarB_";
1013 fillthis+=ptbin;
ea0d8716 1014 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
4e61a020 1015 if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
40445ada 1016
1017 fillthis="hd0d0B_";
1018 fillthis+=ptbin;
1019 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1020
a2121012 1021 if(d0[0] != -99 && d0[1] != -99 ){
1022 fillthis="hd0d0vB_";
1023 fillthis+=ptbin;
1024 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1025 }
1026
40445ada 1027 fillthis="hcosthetapointB_";
1028 fillthis+=ptbin;
1029 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1030
1031 fillthis="hcosthpointd0d0B_";
1032 fillthis+=ptbin;
1033 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
a2121012 1034
1035 fillthis="hdeclB_";
1036 fillthis+=ptbin;
1037 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1038
1039 fillthis="hnormdeclB_";
1040 fillthis+=ptbin;
1041 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
1042
1043 if(vtxProp) {
1044 fillthis="hdeclvB_";
1045 fillthis+=ptbin;
1046 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
1047
1048 fillthis="hnormdeclvB_";
1049 fillthis+=ptbin;
1050 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
1051 }
1052 }//mass cut
4e61a020 1053 }//else (background)
40445ada 1054 }
49061176 1055 return;
1056}
ea0d8716 1057//____________________________________________________________________________
1058
1059void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
49061176 1060 //
40445ada 1061 // function used in UserExec to fill mass histograms:
49061176 1062 //
9de8c723 1063
1064
1065 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1066
3cc4604b 1067 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
ea0d8716 1068 //cout<<"is selected = "<<isSelected<<endl;
feb73eca 1069
ea0d8716 1070 //cout<<"check cuts = "<<endl;
1071 //cuts->PrintAll();
1072 if (!isSelected){
1073 //cout<<"Not Selected"<<endl;
1074 return;
1075 }
9de8c723 1076
4e61a020 1077 if(fDebug>2) cout<<"Candidate selected"<<endl;
a41f6fad 1078
ea0d8716 1079 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1080 //printf("SELECTED\n");
1081 Int_t ptbin=cuts->PtBin(part->Pt());
9de8c723 1082
ea0d8716 1083 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 1084 if(!prong) {
1085 AliDebug(2,"No daughter found");
1086 return;
1087 }
ea0d8716 1088 else{
4b4a1d25 1089 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
ea0d8716 1090 }
7646d6da 1091
ea0d8716 1092 TString fillthis="";
1093 Int_t pdgDgD0toKpi[2]={321,211};
1094 Int_t labD0=-1;
1095 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)
1096
1097 //count candidates selected by cuts
1098 fNentries->Fill(1);
1099 //count true D0 selected by cuts
1100 if (fReadMC && labD0>=0) fNentries->Fill(2);
1101 //PostData(3,fNentries);
1102
1103 if (isSelected==1 || isSelected==3) { //D0
1104 fillthis="histMass_";
1105 fillthis+=ptbin;
1106 //cout<<"Filling "<<fillthis<<endl;
1107
1108 //printf("Fill mass with D0");
1109 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
feb73eca 1110
ea0d8716 1111 if(labD0>=0) {
1112 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1113
1114 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1115 Int_t pdgD0 = partD0->GetPdgCode();
1116 //cout<<"pdg = "<<pdgD0<<endl;
1117 if (pdgD0==421){ //D0
1118 //cout<<"Fill S with D0"<<endl;
1119 fillthis="histSgn_";
1120 fillthis+=ptbin;
1121 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1122 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1123 fillthis="histSgn27_";
46c96ce5 1124 fillthis+=ptbin;
1125 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
a4ae02cd 1126 }
ea0d8716 1127 } else{ //it was a D0bar
1128 fillthis="histRfl_";
a4ae02cd 1129 fillthis+=ptbin;
1130 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1131 }
ea0d8716 1132 } else {//background
1133 fillthis="histBkg_";
a4ae02cd 1134 fillthis+=ptbin;
ea0d8716 1135 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1136 }
1137
1138 }
1139 if (isSelected>1) { //D0bar
1140 fillthis="histMass_";
1141 fillthis+=ptbin;
1142 //printf("Fill mass with D0bar");
1143 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
feb73eca 1144
ea0d8716 1145 if(labD0>=0) {
1146 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1147 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1148 Int_t pdgD0 = partD0->GetPdgCode();
1149 //cout<<" pdg = "<<pdgD0<<endl;
1150 if (pdgD0==-421){ //D0bar
1151 fillthis="histSgn_";
1152 fillthis+=ptbin;
1153 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1154 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1155 fillthis="histSgn27_";
a4ae02cd 1156 fillthis+=ptbin;
1157 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
ea0d8716 1158 }
9de8c723 1159
a4ae02cd 1160
ea0d8716 1161 } else{
1162 fillthis="histRfl_";
a4ae02cd 1163 fillthis+=ptbin;
1164 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1165 }
ea0d8716 1166 } else {//background or LS
1167 fillthis="histBkg_";
1168 fillthis+=ptbin;
1169 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
a4ae02cd 1170 }
ea0d8716 1171 }
a4ae02cd 1172
40445ada 1173 return;
49061176 1174}
4e61a020 1175
1176//__________________________________________________________________________
1177AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1178 //Calculate the primary vertex w/o the daughter tracks of the candidate
1179
1180 AliESDVertex *vertexESD=0x0;
1181 AliAODVertex *vertexAOD=0x0;
1182 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1183
1184 Int_t skipped[2];
1185 Int_t nTrksToSkip=2;
1186 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1187 if(!dgTrack){
1188 AliDebug(2,"no daughter found!");
1189 return 0x0;
1190 }
1191 skipped[0]=dgTrack->GetID();
1192 dgTrack = (AliAODTrack*)d->GetDaughter(1);
1193 if(!dgTrack){
1194 AliDebug(2,"no daughter found!");
1195 return 0x0;
1196 }
1197 skipped[1]=dgTrack->GetID();
1198
1199
1200 //
1201 vertexer->SetSkipTracks(nTrksToSkip,skipped);
4e61a020 1202 vertexer->SetMinClusters(4);
a2121012 1203 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
4e61a020 1204 if(!vertexESD) return vertexAOD;
1205 if(vertexESD->GetNContributors()<=0) {
1206 AliDebug(2,"vertexing failed");
1207 delete vertexESD; vertexESD=NULL;
1208 return vertexAOD;
1209 }
1210
1211 delete vertexer; vertexer=NULL;
1212
1213
1214 // convert to AliAODVertex
1215 Double_t pos[3],cov[6],chi2perNDF;
1216 vertexESD->GetXYZ(pos); // position
1217 vertexESD->GetCovMatrix(cov); //covariance matrix
1218 chi2perNDF = vertexESD->GetChi2toNDF();
1219 delete vertexESD; vertexESD=NULL;
1220
1221 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1222 return vertexAOD;
1223
1224}
1225
1226
49061176 1227//________________________________________________________________________
1228void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1229{
1230 // Terminate analysis
1231 //
1232 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1233
6321ee46 1234
ea0d8716 1235 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1236 if (!fOutputMass) {
1237 printf("ERROR: fOutputMass not available\n");
a4ae02cd 1238 return;
1239 }
ea0d8716 1240 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1241 if (!fDistr) {
1242 printf("ERROR: fDistr not available\n");
a41f6fad 1243 return;
1244 }
5b2e5fae 1245
40445ada 1246 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
5b2e5fae 1247
40445ada 1248 if(!fNentries){
1249 printf("ERROR: fNEntries not available\n");
1250 return;
1251 }
5b2e5fae 1252
ea0d8716 1253 fChecks = dynamic_cast<TList*> (GetOutputData(4));
34137226 1254 if (!fChecks) {
1255 printf("ERROR: fChecks not available\n");
1256 return;
1257 }
ea0d8716 1258
4e61a020 1259
1260 for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1261
4e61a020 1262
1263 if(fArray==1){
feb73eca 1264 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1265
1266
ea0d8716 1267 if(fLsNormalization>1e-6) {
9de8c723 1268
feb73eca 1269 TString massName="histMass_";
ea0d8716 1270 massName+=ipt;
1271 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1272
feb73eca 1273 }
40445ada 1274
feb73eca 1275
40445ada 1276 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1277
ea0d8716 1278 if(fLsNormalization>1e-6) {
40445ada 1279
1280 TString nameDistr="hptB_";
ea0d8716 1281 nameDistr+=ipt;
40445ada 1282 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1283 nameDistr="hdcaB_";
ea0d8716 1284 nameDistr+=ipt;
40445ada 1285 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1286 nameDistr="hcosthetastarB_";
ea0d8716 1287 nameDistr+=ipt;
40445ada 1288 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1289 nameDistr="hd0B_";
ea0d8716 1290 nameDistr+=ipt;
40445ada 1291 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1292 nameDistr="hd0d0B_";
ea0d8716 1293 nameDistr+=ipt;
40445ada 1294 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1295 nameDistr="hcosthetapointB_";
ea0d8716 1296 nameDistr+=ipt;
40445ada 1297 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1298 nameDistr="hcosthpointd0d0B_";
ea0d8716 1299 nameDistr+=ipt;
40445ada 1300 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
feb73eca 1301
40445ada 1302 }
feb73eca 1303 }
1304 }
4e61a020 1305 TString cvname,cstname;
527f330b 1306
1307 if (fArray==0){
1308 cvname="D0invmass";
4e61a020 1309 cstname="cstat0";
1310 } else {
1311 cvname="LSinvmass";
1312 cstname="cstat1";
1313 }
527f330b 1314
34137226 1315 TCanvas *cMass=new TCanvas(cvname,cvname);
1316 cMass->cd();
ea0d8716 1317 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
527f330b 1318
4e61a020 1319 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
34137226 1320 cStat->cd();
1321 cStat->SetGridy();
1322 fNentries->Draw("htext0");
527f330b 1323
4e61a020 1324 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1325 // ccheck->cd();
1326
49061176 1327 return;
1328}
1329