]> 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
In case the input is MC Kinematics, need some patches in the vertex and extraction...
[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
7c23877d 579 // fix for temporary bug in ESDfilter
580 // the AODs with null vertex pointer didn't pass the PhysSel
581 if(!aod->GetPrimaryVertex()) return;
582
ce39f0ac 583 TClonesArray *mcArray = 0;
584 AliAODMCHeader *mcHeader = 0;
585
586 if(fReadMC) {
587 // load MC particles
588 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
589 if(!mcArray) {
590 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
591 return;
592 }
40445ada 593
ce39f0ac 594 // load MC header
595 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
596 if(!mcHeader) {
597 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
598 return;
599 }
49061176 600 }
601
b557eb43 602 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
40445ada 603
a4ae02cd 604 //histogram filled with 1 for every AOD
34137226 605 fNentries->Fill(0);
49061176 606
0108fa62 607
40445ada 608
609 // AOD primary vertex
610 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
c59e1898 611 if(!vtx1) {
612 AliError("There is no primary vertex !");
613 return;
614 }
b272aebf 615
40445ada 616 Bool_t isGoodVtx=kFALSE;
b272aebf 617
40445ada 618 //vtx1->Print();
619 TString primTitle = vtx1->GetTitle();
620 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
621 isGoodVtx=kTRUE;
622 fNentries->Fill(3);
623 }
624
625 //cout<<"Start checks"<<endl;
626 Int_t ntracks=0,isGoodTrack=0;
627
628 if(aod) ntracks=aod->GetNTracks();
629
75638da0 630 //cout<<"ntracks = "<<ntracks<<endl;
ea0d8716 631
75638da0 632 //loop on tracks in the event
633 for (Int_t k=0;k<ntracks;k++){
634 AliAODTrack* track=aod->GetTrack(k);
4e61a020 635
636 //check TOF
637
638 if(!(track->GetStatus()&AliESDtrack::kTPCrefit &&
639 track->GetStatus()&AliESDtrack::kITSrefit &&
640 track->GetTPCNcls() >=70 &&
641 track->GetStatus()&AliESDtrack::kTOFpid &&
642 track->GetStatus()&AliESDtrack::kTOFout &&
643 track->GetStatus()&AliESDtrack::kTIME)) continue;
644 AliAODPid *pid = track->GetDetPid();
a2121012 645 if(!pid) {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
4e61a020 646
647 Double_t times[5];
648 pid->GetIntegratedTimes(times);
649
650 ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
651 ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
652
653 ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
654 if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
655
656
75638da0 657 //check clusters of the tracks
658 Int_t nclsTot=0,nclsSPD=0;
40445ada 659
75638da0 660 for(Int_t l=0;l<6;l++) {
661 if(TESTBIT(track->GetITSClusterMap(),l)) {
662 nclsTot++; if(l<2) nclsSPD++;
b272aebf 663 }
75638da0 664 }
40445ada 665
75638da0 666 if (track->Pt()>0.3 &&
667 track->GetStatus()&AliESDtrack::kTPCrefit &&
668 track->GetStatus()&AliESDtrack::kITSrefit &&
669 nclsTot>3 &&
670 nclsSPD>0) {//fill hist good tracks
a2121012 671
75638da0 672 ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
673 isGoodTrack++;
40445ada 674 }
75638da0 675 //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
676 ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
677 }
678 //number of events with good vertex and at least 2 good tracks
679 if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
ea0d8716 680
681
75638da0 682 // loop over candidates
683 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4e61a020 684 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
75638da0 685
75638da0 686 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
687 //Int_t nPosPairs=0, nNegPairs=0;
688 //cout<<"inside the loop"<<endl;
689 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
4e61a020 690
691 //check daughters
692 if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
693 AliDebug(1,"at least one daughter not found!");
694 fNentries->Fill(6);
695 return;
696 }
697
75638da0 698 Bool_t unsetvtx=kFALSE;
699 if(!d->GetOwnPrimaryVtx()) {
700 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
701 unsetvtx=kTRUE;
702 }
ea0d8716 703
4e61a020 704
75638da0 705 //check reco daughter in acceptance
4e61a020 706 /*
75638da0 707 Double_t eta0=d->EtaProng(0);
708 Double_t eta1=d->EtaProng(1);
4e61a020 709 */
710 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
711 //if( TMath::Abs(eta0)<0.9 && TMath::Abs(eta1)<0.9 ){
712 //apply cuts on tracks
3cc4604b 713 Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
a2121012 714 if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(1))->GetTPCNcls() < 70) isSelected=kFALSE;
3cc4604b 715 if (!isSelected) return;
4e61a020 716 fNentries->Fill(7);
717 if(fDebug>2) cout<<"tracks selected"<<endl;
3cc4604b 718
4e61a020 719 FillVarHists(aod,d,mcArray,fCuts,fDistr);
ea0d8716 720 FillMassHists(d,mcArray,fCuts,fOutputMass);
75638da0 721 }
ea0d8716 722
75638da0 723 if(unsetvtx) d->UnsetOwnPrimaryVtx();
724 } //end for prongs
725
ea0d8716 726
40445ada 727 // Post the data
ea0d8716 728 PostData(1,fOutputMass);
729 PostData(2,fDistr);
40445ada 730 PostData(3,fNentries);
ea0d8716 731 PostData(4,fChecks);
732
40445ada 733 return;
734}
b272aebf 735
40445ada 736//____________________________________________________________________________
4e61a020 737void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
40445ada 738 //
739 // function used in UserExec to fill variable histograms:
740 //
b272aebf 741
40445ada 742 //add distr here
743 UInt_t pdgs[2];
744
745 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
746 pdgs[0]=211;
747 pdgs[1]=321;
748 Double_t minvD0 = part->InvMassD0();
749 pdgs[1]=211;
750 pdgs[0]=321;
751 Double_t minvD0bar = part->InvMassD0bar();
752 //cout<<"inside mass cut"<<endl;
b272aebf 753
40445ada 754 Int_t pdgDgD0toKpi[2]={321,211};
755 Int_t lab=-9999;
756 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)
757 //Double_t pt = d->Pt(); //mother pt
758 Bool_t isSelected=kFALSE;
b272aebf 759
3cc4604b 760
ea0d8716 761 if(fCutOnDistr){
3cc4604b 762 isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate);
ea0d8716 763 if (!isSelected){
764 //cout<<"Not Selected"<<endl;
765 return;
766 }
767 }
0108fa62 768
4e61a020 769 Double_t invmasscut=0.03;
770
40445ada 771 TString fillthispi="",fillthisK="",fillthis="";
b272aebf 772
ea0d8716 773 Int_t ptbin=cuts->PtBin(part->Pt());
4e61a020 774 if(ptbin==-1) {fNentries->Fill(5); return;} //out of bounds
775
776 //recalculate vertex
777 AliAODVertex *vtxProp=0x0;
778 vtxProp=GetPrimaryVtxSkipped(aod,part);
779 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
780 dz1[0]=-99; dz2[0]=-99;
781 if(vtxProp) {
782 part->SetOwnPrimaryVtx(vtxProp);
783 //Bool_t unsetvtx=kTRUE;
784 //Calculate d0 for daughter tracks with Vtx Skipped
785 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)part->GetDaughter(0));
786 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)part->GetDaughter(1));
787 esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
788 esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
789 delete esdtrack1;
790 delete esdtrack2;
791 }
792
793 Double_t d0[2]={dz1[0],dz2[0]};
794
40445ada 795 if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
ea0d8716 796 //printf("\nif no cuts or cuts passed\n");
40445ada 797 if(lab>=0 && fReadMC){ //signal
0108fa62 798
40445ada 799 //check pdg of the prongs
800 AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
801 AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
4e61a020 802 if(!prong0 || !prong1) {
803 return;
804 }
40445ada 805 Int_t labprong[2];
806 labprong[0]=prong0->GetLabel();
807 labprong[1]=prong1->GetLabel();
808 AliAODMCParticle *mcprong=0;
809 Int_t pdgProng[2]={0,0};
4e61a020 810
40445ada 811 for (Int_t iprong=0;iprong<2;iprong++){
812 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
813 pdgProng[iprong]=mcprong->GetPdgCode();
0108fa62 814 }
815
40445ada 816 //no mass cut ditributions: ptbis
817
818 fillthispi="hptpiSnoMcut_";
819 fillthispi+=ptbin;
820
821 fillthisK="hptKSnoMcut_";
822 fillthisK+=ptbin;
823
824 if (TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321){
825 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
826 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
827 }else {
828 if (TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211){
829 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
830 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
831 }
832 }
833
834 //no mass cut ditributions: mass
835 fillthis="hMassS_";
836 fillthis+=ptbin;
837
838 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
839 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
840 }
841 else { //D0bar
842 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
843 }
4e61a020 844
0108fa62 845 //apply cut on invariant mass on the pair
4e61a020 846 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
40445ada 847
0108fa62 848 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
849 for (Int_t iprong=0; iprong<2; iprong++){
40445ada 850 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
04e09ffc 851 labprong[iprong]=prong->GetLabel();
40445ada 852
0108fa62 853 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
40445ada 854 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
0108fa62 855 Int_t pdgprong=mcprong->GetPdgCode();
40445ada 856
0108fa62 857 if(TMath::Abs(pdgprong)==211) {
858 //cout<<"pi"<<endl;
40445ada 859
860 fillthispi="hptpiS_";
861 fillthispi+=ptbin;
862 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
4e61a020 863 fillthispi="hd0piS_";
864 fillthispi+=ptbin;
865 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
866 if(d0[iprong] != -99) {
867
868 fillthispi="hd0vpiS_";
869 fillthispi+=ptbin;
870 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
871 }
872
b272aebf 873 }
40445ada 874
875 if(TMath::Abs(pdgprong)==321) {
876 //cout<<"kappa"<<endl;
877
878 fillthisK="hptKS_";
879 fillthisK+=ptbin;
880 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
40445ada 881 fillthisK="hd0KS_";
882 fillthisK+=ptbin;
883 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
4e61a020 884 if (d0[iprong] != -99){
885 fillthisK="hd0vKS_";
886 fillthisK+=ptbin;
887 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
888 }
b272aebf 889 }
b272aebf 890
4e61a020 891 fillthis="hcosthpointd0S_";
40445ada 892 fillthis+=ptbin;
4e61a020 893 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
b272aebf 894
4e61a020 895 } //end loop on prongs
b272aebf 896
4e61a020 897 fillthis="hdcaS_";
898 fillthis+=ptbin;
899 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
900
901 fillthis="hcosthetapointS_";
902 fillthis+=ptbin;
903 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
904
905 fillthis="hcosthpointd0d0S_";
906 fillthis+=ptbin;
907 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
b272aebf 908
4e61a020 909 fillthis="hcosthetastarS_";
910 fillthis+=ptbin;
911 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
912 else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
0108fa62 913
4e61a020 914 fillthis="hd0d0S_";
915 fillthis+=ptbin;
916 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
917
918 if(d0[0] != -99 && d0[1] != -99 ){
919 fillthis="hd0d0vS_";
40445ada 920 fillthis+=ptbin;
4e61a020 921 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
b272aebf 922 }
4e61a020 923
924 fillthis="hdeclS_";
925 fillthis+=ptbin;
926 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
927
928 fillthis="hnormdeclS_";
929 fillthis+=ptbin;
930 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
a2121012 931
932 if(vtxProp) {
933 fillthis="hdeclvS_";
934 fillthis+=ptbin;
935 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
936
937 fillthis="hnormdeclvS_";
938 fillthis+=ptbin;
939 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
940 }
941 } //end mass cut
4e61a020 942
0108fa62 943 } else{ //Background or LS
944 //cout<<"is background"<<endl;
40445ada 945
b272aebf 946 //no mass cut distributions: mass, ptbis
40445ada 947 fillthis="hMassB_";
948 fillthis+=ptbin;
ea0d8716 949
950 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
951 if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
b272aebf 952
40445ada 953 fillthis="hptB1prongnoMcut_";
954 fillthis+=ptbin;
955
956 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
957
958 fillthis="hptB2prongsnoMcut_";
959 fillthis+=ptbin;
960 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
961 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 962
963 //apply cut on invariant mass on the pair
4e61a020 964 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
0108fa62 965
966
40445ada 967 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 968 if(!prong) {
969 if(fDebug>2) cout<<"No daughter found";
970 return;
971 }
0108fa62 972 else{
99012cda 973 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
0108fa62 974 }
40445ada 975
976 //normalise pt distr to half afterwards
977 fillthis="hptB_";
978 fillthis+=ptbin;
4e61a020 979 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
980 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
0108fa62 981
4e61a020 982 fillthis="hd0p0B_";
40445ada 983 fillthis+=ptbin;
984 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
4e61a020 985 fillthis="hd0p1B_";
986 fillthis+=ptbin;
987 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
a2121012 988 fillthis="hd0B_";
989 fillthis+=ptbin;
990 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
991 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
4e61a020 992
993 fillthis="hd0vp0B_";
994 fillthis+=ptbin;
995 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
996 fillthis="hd0vp1B_";
997 fillthis+=ptbin;
998 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
999
a2121012 1000 fillthis="hd0vB_";
1001 fillthis+=ptbin;
1002 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1003 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1004
1005
4e61a020 1006 fillthis="hcosthpointd0B_";
1007 fillthis+=ptbin;
1008 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
1009 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
1010
ce39f0ac 1011
40445ada 1012 fillthis="hdcaB_";
1013 fillthis+=ptbin;
1014 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
34137226 1015
40445ada 1016 fillthis="hcosthetastarB_";
1017 fillthis+=ptbin;
ea0d8716 1018 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
4e61a020 1019 if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
40445ada 1020
1021 fillthis="hd0d0B_";
1022 fillthis+=ptbin;
1023 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1024
a2121012 1025 if(d0[0] != -99 && d0[1] != -99 ){
1026 fillthis="hd0d0vB_";
1027 fillthis+=ptbin;
1028 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1029 }
1030
40445ada 1031 fillthis="hcosthetapointB_";
1032 fillthis+=ptbin;
1033 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1034
1035 fillthis="hcosthpointd0d0B_";
1036 fillthis+=ptbin;
1037 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
a2121012 1038
1039 fillthis="hdeclB_";
1040 fillthis+=ptbin;
1041 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1042
1043 fillthis="hnormdeclB_";
1044 fillthis+=ptbin;
1045 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
1046
1047 if(vtxProp) {
1048 fillthis="hdeclvB_";
1049 fillthis+=ptbin;
1050 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
1051
1052 fillthis="hnormdeclvB_";
1053 fillthis+=ptbin;
1054 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
1055 }
1056 }//mass cut
4e61a020 1057 }//else (background)
40445ada 1058 }
49061176 1059 return;
1060}
ea0d8716 1061//____________________________________________________________________________
1062
1063void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
49061176 1064 //
40445ada 1065 // function used in UserExec to fill mass histograms:
49061176 1066 //
9de8c723 1067
1068
1069 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1070
3cc4604b 1071 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
ea0d8716 1072 //cout<<"is selected = "<<isSelected<<endl;
feb73eca 1073
ea0d8716 1074 //cout<<"check cuts = "<<endl;
1075 //cuts->PrintAll();
1076 if (!isSelected){
1077 //cout<<"Not Selected"<<endl;
1078 return;
1079 }
9de8c723 1080
4e61a020 1081 if(fDebug>2) cout<<"Candidate selected"<<endl;
a41f6fad 1082
ea0d8716 1083 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1084 //printf("SELECTED\n");
1085 Int_t ptbin=cuts->PtBin(part->Pt());
9de8c723 1086
ea0d8716 1087 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
4e61a020 1088 if(!prong) {
1089 AliDebug(2,"No daughter found");
1090 return;
1091 }
ea0d8716 1092 else{
4b4a1d25 1093 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
ea0d8716 1094 }
7646d6da 1095
ea0d8716 1096 TString fillthis="";
1097 Int_t pdgDgD0toKpi[2]={321,211};
1098 Int_t labD0=-1;
1099 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)
1100
1101 //count candidates selected by cuts
1102 fNentries->Fill(1);
1103 //count true D0 selected by cuts
1104 if (fReadMC && labD0>=0) fNentries->Fill(2);
1105 //PostData(3,fNentries);
1106
1107 if (isSelected==1 || isSelected==3) { //D0
1108 fillthis="histMass_";
1109 fillthis+=ptbin;
1110 //cout<<"Filling "<<fillthis<<endl;
1111
1112 //printf("Fill mass with D0");
1113 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
feb73eca 1114
ea0d8716 1115 if(labD0>=0) {
1116 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1117
1118 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1119 Int_t pdgD0 = partD0->GetPdgCode();
1120 //cout<<"pdg = "<<pdgD0<<endl;
1121 if (pdgD0==421){ //D0
1122 //cout<<"Fill S with D0"<<endl;
1123 fillthis="histSgn_";
1124 fillthis+=ptbin;
1125 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1126 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1127 fillthis="histSgn27_";
46c96ce5 1128 fillthis+=ptbin;
1129 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
a4ae02cd 1130 }
ea0d8716 1131 } else{ //it was a D0bar
1132 fillthis="histRfl_";
a4ae02cd 1133 fillthis+=ptbin;
1134 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1135 }
ea0d8716 1136 } else {//background
1137 fillthis="histBkg_";
a4ae02cd 1138 fillthis+=ptbin;
ea0d8716 1139 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1140 }
1141
1142 }
1143 if (isSelected>1) { //D0bar
1144 fillthis="histMass_";
1145 fillthis+=ptbin;
1146 //printf("Fill mass with D0bar");
1147 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
feb73eca 1148
ea0d8716 1149 if(labD0>=0) {
1150 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1151 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1152 Int_t pdgD0 = partD0->GetPdgCode();
1153 //cout<<" pdg = "<<pdgD0<<endl;
1154 if (pdgD0==-421){ //D0bar
1155 fillthis="histSgn_";
1156 fillthis+=ptbin;
1157 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1158 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1159 fillthis="histSgn27_";
a4ae02cd 1160 fillthis+=ptbin;
1161 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
ea0d8716 1162 }
9de8c723 1163
a4ae02cd 1164
ea0d8716 1165 } else{
1166 fillthis="histRfl_";
a4ae02cd 1167 fillthis+=ptbin;
1168 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1169 }
ea0d8716 1170 } else {//background or LS
1171 fillthis="histBkg_";
1172 fillthis+=ptbin;
1173 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
a4ae02cd 1174 }
ea0d8716 1175 }
a4ae02cd 1176
40445ada 1177 return;
49061176 1178}
4e61a020 1179
1180//__________________________________________________________________________
1181AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1182 //Calculate the primary vertex w/o the daughter tracks of the candidate
1183
1184 AliESDVertex *vertexESD=0x0;
1185 AliAODVertex *vertexAOD=0x0;
1186 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1187
1188 Int_t skipped[2];
1189 Int_t nTrksToSkip=2;
1190 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1191 if(!dgTrack){
1192 AliDebug(2,"no daughter found!");
1193 return 0x0;
1194 }
1195 skipped[0]=dgTrack->GetID();
1196 dgTrack = (AliAODTrack*)d->GetDaughter(1);
1197 if(!dgTrack){
1198 AliDebug(2,"no daughter found!");
1199 return 0x0;
1200 }
1201 skipped[1]=dgTrack->GetID();
1202
1203
1204 //
1205 vertexer->SetSkipTracks(nTrksToSkip,skipped);
4e61a020 1206 vertexer->SetMinClusters(4);
a2121012 1207 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
4e61a020 1208 if(!vertexESD) return vertexAOD;
1209 if(vertexESD->GetNContributors()<=0) {
1210 AliDebug(2,"vertexing failed");
1211 delete vertexESD; vertexESD=NULL;
1212 return vertexAOD;
1213 }
1214
1215 delete vertexer; vertexer=NULL;
1216
1217
1218 // convert to AliAODVertex
1219 Double_t pos[3],cov[6],chi2perNDF;
1220 vertexESD->GetXYZ(pos); // position
1221 vertexESD->GetCovMatrix(cov); //covariance matrix
1222 chi2perNDF = vertexESD->GetChi2toNDF();
1223 delete vertexESD; vertexESD=NULL;
1224
1225 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1226 return vertexAOD;
1227
1228}
1229
1230
49061176 1231//________________________________________________________________________
1232void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1233{
1234 // Terminate analysis
1235 //
1236 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1237
6321ee46 1238
ea0d8716 1239 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1240 if (!fOutputMass) {
1241 printf("ERROR: fOutputMass not available\n");
a4ae02cd 1242 return;
1243 }
ea0d8716 1244 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1245 if (!fDistr) {
1246 printf("ERROR: fDistr not available\n");
a41f6fad 1247 return;
1248 }
5b2e5fae 1249
40445ada 1250 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
5b2e5fae 1251
40445ada 1252 if(!fNentries){
1253 printf("ERROR: fNEntries not available\n");
1254 return;
1255 }
5b2e5fae 1256
ea0d8716 1257 fChecks = dynamic_cast<TList*> (GetOutputData(4));
34137226 1258 if (!fChecks) {
1259 printf("ERROR: fChecks not available\n");
1260 return;
1261 }
ea0d8716 1262
4e61a020 1263
1264 for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1265
4e61a020 1266
1267 if(fArray==1){
feb73eca 1268 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1269
1270
ea0d8716 1271 if(fLsNormalization>1e-6) {
9de8c723 1272
feb73eca 1273 TString massName="histMass_";
ea0d8716 1274 massName+=ipt;
1275 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1276
feb73eca 1277 }
40445ada 1278
feb73eca 1279
40445ada 1280 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1281
ea0d8716 1282 if(fLsNormalization>1e-6) {
40445ada 1283
1284 TString nameDistr="hptB_";
ea0d8716 1285 nameDistr+=ipt;
40445ada 1286 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1287 nameDistr="hdcaB_";
ea0d8716 1288 nameDistr+=ipt;
40445ada 1289 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1290 nameDistr="hcosthetastarB_";
ea0d8716 1291 nameDistr+=ipt;
40445ada 1292 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1293 nameDistr="hd0B_";
ea0d8716 1294 nameDistr+=ipt;
40445ada 1295 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1296 nameDistr="hd0d0B_";
ea0d8716 1297 nameDistr+=ipt;
40445ada 1298 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1299 nameDistr="hcosthetapointB_";
ea0d8716 1300 nameDistr+=ipt;
40445ada 1301 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1302 nameDistr="hcosthpointd0d0B_";
ea0d8716 1303 nameDistr+=ipt;
40445ada 1304 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
feb73eca 1305
40445ada 1306 }
feb73eca 1307 }
1308 }
4e61a020 1309 TString cvname,cstname;
527f330b 1310
1311 if (fArray==0){
1312 cvname="D0invmass";
4e61a020 1313 cstname="cstat0";
1314 } else {
1315 cvname="LSinvmass";
1316 cstname="cstat1";
1317 }
527f330b 1318
34137226 1319 TCanvas *cMass=new TCanvas(cvname,cvname);
1320 cMass->cd();
ea0d8716 1321 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
527f330b 1322
4e61a020 1323 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
34137226 1324 cStat->cd();
1325 cStat->SetGridy();
1326 fNentries->Draw("htext0");
527f330b 1327
4e61a020 1328 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1329 // ccheck->cd();
1330
49061176 1331 return;
1332}
1333