]> 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
Update and addition of LS analysis (Renu, Giacomo, Francesco)
[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>
28#include <TNtuple.h>
29#include <TList.h>
30#include <TH1F.h>
a41f6fad 31#include <TH2F.h>
32#include <TDatabasePDG.h>
49061176 33
b557eb43 34#include "AliAnalysisManager.h"
35#include "AliAODHandler.h"
49061176 36#include "AliAODEvent.h"
37#include "AliAODVertex.h"
38#include "AliAODTrack.h"
39#include "AliAODMCHeader.h"
40#include "AliAODMCParticle.h"
41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliAODRecoCascadeHF.h"
43#include "AliAnalysisVertexingHF.h"
44#include "AliAnalysisTaskSE.h"
45#include "AliAnalysisTaskSED0Mass.h"
46
b557eb43 47
49061176 48ClassImp(AliAnalysisTaskSED0Mass)
49
50
51//________________________________________________________________________
52AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
53AliAnalysisTaskSE(),
9de8c723 54fOutputPPR(0),
a41f6fad 55fOutputloose(0),
56fDistr(0),
feb73eca 57fNentries(0),
9de8c723 58fVHFPPR(0),
6ce5994f 59fVHFloose(0),
feb73eca 60fArray(0),
ce39f0ac 61fReadMC(0),
feb73eca 62fLsNormalization(1.)
63
49061176 64{
65 // Default constructor
9de8c723 66 for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
49061176 67}
68
69//________________________________________________________________________
70AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
71AliAnalysisTaskSE(name),
9de8c723 72fOutputPPR(0),
a4ae02cd 73fOutputloose(0),
a41f6fad 74fDistr(0),
feb73eca 75fNentries(0),
9de8c723 76fVHFPPR(0),
6ce5994f 77fVHFloose(0),
feb73eca 78fArray(0),
ce39f0ac 79fReadMC(0),
feb73eca 80fLsNormalization(1.)
49061176 81{
82 // Default constructor
9de8c723 83 for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
49061176 84
85 // Output slot #1 writes into a TList container
86 DefineOutput(1,TList::Class()); //My private output
a4ae02cd 87 // Output slot #2 writes into a TList container
88 DefineOutput(2,TList::Class()); //My private output
89 // Output slot #3 writes into a TH1F container
90 DefineOutput(3,TH1F::Class()); //My private output
a41f6fad 91 // Output slot #4 writes into a TList container
92 DefineOutput(4,TList::Class()); //My private output
49061176 93}
94
95//________________________________________________________________________
96AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
97{
9de8c723 98 if (fOutputPPR) {
99 delete fOutputPPR;
100 fOutputPPR = 0;
a4ae02cd 101 }
9de8c723 102 if (fVHFPPR) {
103 delete fVHFPPR;
104 fVHFPPR = 0;
a4ae02cd 105 }
106 if (fOutputloose) {
107 delete fOutputloose;
108 fOutputloose = 0;
109 }
a41f6fad 110
111 if (fDistr) {
112 delete fDistr;
113 fDistr = 0;
114 }
115
a4ae02cd 116 if (fVHFloose) {
117 delete fVHFloose;
118 fVHFloose = 0;
49061176 119 }
46c96ce5 120 if (fNentries){
121 delete fNentries;
122 fNentries = 0;
123 }
a41f6fad 124
125
49061176 126}
127
128//________________________________________________________________________
129void AliAnalysisTaskSED0Mass::Init()
130{
131 // Initialization
132
133 if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
134
135 gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
136
a4ae02cd 137 // 2 sets of dedidcated cuts -- defined in UserExec
9de8c723 138 fVHFPPR = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
a4ae02cd 139 fVHFloose = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
140
49061176 141 return;
142}
143
144//________________________________________________________________________
145void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
146{
147
148 // Create the output container
149 //
150 if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
151
152 // Several histograms are more conveniently managed in a TList
9de8c723 153 fOutputPPR = new TList();
154 fOutputPPR->SetOwner();
155 fOutputPPR->SetName("listPPR");
a4ae02cd 156
157 fOutputloose = new TList();
158 fOutputloose->SetOwner();
159 fOutputloose->SetName("listloose");
49061176 160
a41f6fad 161 fDistr = new TList();
162 fDistr->SetOwner();
163 fDistr->SetName("distributionslist");
164
9de8c723 165 const Int_t nhist=5;
a4ae02cd 166
9de8c723 167 TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ", namedistr=" ";
49061176 168
169 for(Int_t i=0;i<nhist;i++){
7646d6da 170 nameMass="histMass_";
49061176 171 nameMass+=i+1;
9de8c723 172 nameSgn27="histSgn27_";
173 nameSgn27+=i+1;
7646d6da 174 nameSgn="histSgn_";
49061176 175 nameSgn+=i+1;
7646d6da 176 nameBkg="histBkg_";
49061176 177 nameBkg+=i+1;
a4ae02cd 178 nameRfl="histRfl_";
179 nameRfl+=i+1;
7646d6da 180
a41f6fad 181 //histograms of invariant mass distributions
182
46c96ce5 183 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
184 TH1F *tmpMl=(TH1F*)tmpMt->Clone();
185 tmpMt->Sumw2();
186 tmpMl->Sumw2();
a4ae02cd 187
9de8c723 188 //to compare with AliAnalysisTaskCharmFraction
189 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);
190 TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
191 tmpS27t->Sumw2();
192 tmpS27l->Sumw2();
193
46c96ce5 194 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
195 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
196 tmpSt->Sumw2();
197 tmpSl->Sumw2();
a4ae02cd 198
46c96ce5 199 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
200 TH1F *tmpBl=(TH1F*)tmpBt->Clone();
201 tmpBt->Sumw2();
202 tmpBl->Sumw2();
a4ae02cd 203
feb73eca 204 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
46c96ce5 205 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
206 TH1F *tmpRl=(TH1F*)tmpRt->Clone();
207 tmpRt->Sumw2();
208 tmpRl->Sumw2();
a4ae02cd 209 // printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
210
9de8c723 211 fOutputPPR->Add(tmpMt);
212 fOutputPPR->Add(tmpSt);
213 fOutputPPR->Add(tmpS27t);
214 fOutputPPR->Add(tmpBt);
215 fOutputPPR->Add(tmpRt);
a4ae02cd 216
46c96ce5 217 fOutputloose->Add(tmpMl);
218 fOutputloose->Add(tmpSl);
9de8c723 219 fOutputloose->Add(tmpS27l);
46c96ce5 220 fOutputloose->Add(tmpBl);
221 fOutputloose->Add(tmpRl);
a4ae02cd 222
7646d6da 223 }
a4ae02cd 224
a41f6fad 225 //histograms of cut variable distributions
226 // pT
227 namedistr="hptpiS";
228 TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
a41f6fad 229
230 namedistr="hptKS";
231 TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
232
233 namedistr="hptB";
234 TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
a41f6fad 235
236 // dca
237 namedistr="hdcaS";
238 TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
239 namedistr="hdcaB";
240 TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
a41f6fad 241
242 // costhetastar
feb73eca 243 namedistr="hcosthetastarS";
a41f6fad 244 TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
feb73eca 245 namedistr="hcosthetastarB";
a41f6fad 246 TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
a41f6fad 247
248 // impact parameter
249 namedistr="hd0piS";
250 TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
a41f6fad 251
252 namedistr="hd0KS";
253 TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
254 namedistr="hd0B";
255 TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution;d0 [cm]",200,-0.1,0.1);
a41f6fad 256
257 namedistr="hd0d0S";
258 TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
259 namedistr="hd0d0B";
260 TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
a41f6fad 261
262 // costhetapoint
263 namedistr="hcosthetapointS";
264 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
265 namedistr="hcosthetapointB";
266 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
a41f6fad 267
268 namedistr="hcosthpointd0d0S";
269 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);
270 namedistr="hcosthpointd0d0B";
271 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);
272
273 fDistr->Add(hptpiS);
a41f6fad 274 fDistr->Add(hptKS);
275 fDistr->Add(hptB);
a41f6fad 276
277 fDistr->Add(hdcaS);
278 fDistr->Add(hdcaB);
a41f6fad 279
280 fDistr->Add(hd0piS);
a41f6fad 281 fDistr->Add(hd0KS);
282 fDistr->Add(hd0B);
a41f6fad 283
284 fDistr->Add(hd0d0S);
285 fDistr->Add(hd0d0B);
a41f6fad 286
287 fDistr->Add(hcosthetastarS);
288 fDistr->Add(hcosthetastarB);
a41f6fad 289
290 fDistr->Add(hcosthetapointS);
291 fDistr->Add(hcosthetapointB);
a41f6fad 292
293 fDistr->Add(hcosthpointd0d0S);
294 fDistr->Add(hcosthpointd0d0B);
295
9de8c723 296 fNentries=new TH1F("nentriesD0", "nentriesD0->Integral(1,2) = number of AODs *** nentriesD0->Integral(3,4) = number of candidates selected with cuts *** nentriesD0->Integral(5,6) = number of D0 selected with cuts", 6,1.,4.);
a4ae02cd 297
49061176 298 return;
299}
300
301//________________________________________________________________________
302void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
303{
304 // Execute analysis for current event:
305 // heavy flavor candidates association to MC truth
a4ae02cd 306 //cout<<"I'm in UserExec"<<endl;
49061176 307 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
feb73eca 308
b557eb43 309 TString bname;
feb73eca 310 if(fArray==0){ //D0 candidates
b557eb43 311 // load D0->Kpi candidates
feb73eca 312 //cout<<"D0 candidates"<<endl;
b557eb43 313 bname="D0toKpi";
feb73eca 314 } else { //LikeSign candidates
feb73eca 315 // load like sign candidates
b557eb43 316 //cout<<"LS candidates"<<endl;
317 bname="LikeSign2Prong";
318 }
319
320 TClonesArray *inputArray=0;
321
322 if(!aod && AODEvent() && IsStandardAOD()) {
323 // In case there is an AOD handler writing a standard AOD, use the AOD
324 // event in memory rather than the input (ESD) event.
325 aod = dynamic_cast<AliAODEvent*> (AODEvent());
326 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
327 // have to taken from the AOD event hold by the AliAODExtension
328 AliAODHandler* aodHandler = (AliAODHandler*)
329 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
330 if(aodHandler->GetExtensions()) {
331 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
332 AliAODEvent *aodFromExt = ext->GetAOD();
333 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
feb73eca 334 }
b557eb43 335 } else {
336 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
337 }
feb73eca 338
b557eb43 339
340 if(!inputArray) {
341 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
342 return;
49061176 343 }
344
345 // AOD primary vertex
346 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
347 //vtx1->Print();
feb73eca 348
ce39f0ac 349 TClonesArray *mcArray = 0;
350 AliAODMCHeader *mcHeader = 0;
351
352 if(fReadMC) {
353 // load MC particles
354 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
355 if(!mcArray) {
356 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
357 return;
358 }
49061176 359
ce39f0ac 360 // load MC header
361 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
362 if(!mcHeader) {
363 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
364 return;
365 }
49061176 366 }
367
b557eb43 368 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
a4ae02cd 369
370 //histogram filled with 1 for every AOD
371 fNentries->Fill(1);
372 PostData(3,fNentries);
9de8c723 373
feb73eca 374 // loop over candidates
375 Int_t nInD0toKpi = inputArray->GetEntriesFast();
4464ce7e 376 if(fDebug>1) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
49061176 377
378 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
feb73eca 379 //Int_t nPosPairs=0, nNegPairs=0;
a41f6fad 380 //cout<<"inside the loop"<<endl;
feb73eca 381 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
49061176 382 Bool_t unsetvtx=kFALSE;
383 if(!d->GetOwnPrimaryVtx()) {
384 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
385 unsetvtx=kTRUE;
386 }
387
9de8c723 388 //check reco daughter in acceptance
389 Double_t eta0=d->EtaProng(0);
390 Double_t eta1=d->EtaProng(1);
391 Bool_t prongsinacc=kFALSE;
392 if (TMath::Abs(eta0) < 0.9 && TMath::Abs(eta1) < 0.9) prongsinacc=kTRUE;
a41f6fad 393
394 //add distr here
395 UInt_t pdgs[2];
49061176 396
a41f6fad 397 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
398 pdgs[0]=211;
399 pdgs[1]=321;
400 Double_t minvD0 = d->InvMassD0();
401 pdgs[1]=211;
402 pdgs[0]=321;
403 Double_t minvD0bar = d->InvMassD0bar();
404 //apply cut on invariant mass on the pair
ce39f0ac 405 if(fReadMC){
406 if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
407 //cout<<"inside mass cut"<<endl;
408 Int_t pdgDgD0toKpi[2]={321,211};
409 Int_t lab=d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
410
411 if(lab>=0){ //signal
412 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
413 for (Int_t iprong=0; iprong<2; iprong++){
414 AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(iprong);
415 Int_t labprong=prong->GetLabel();
416
417 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
418 AliAODMCParticle *mcprong=0;
419 if(labprong>=0) mcprong= (AliAODMCParticle*)mcArray->At(labprong);
420 Int_t pdgprong=mcprong->GetPdgCode();
421 if(TMath::Abs(pdgprong)==211) {
422 //cout<<"pi"<<endl;
423 ((TH1F*)fDistr->FindObject("hptpiS"))->Fill(d->PtProng(iprong));
424 ((TH1F*)fDistr->FindObject("hd0piS"))->Fill(d->Getd0Prong(iprong));
425 }
426
427 if(TMath::Abs(pdgprong)==321) {
428 //cout<<"kappa"<<endl;
429 ((TH1F*)fDistr->FindObject("hptKS"))->Fill(d->PtProng(iprong));
430 ((TH1F*)fDistr->FindObject("hd0KS"))->Fill(d->Getd0Prong(iprong));
431 }
432 ((TH1F*)fDistr->FindObject("hdcaS"))->Fill(d->GetDCA());
433
434 }
435
436 if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
437 ((TH1F*)fDistr->FindObject("hcosthetastarS"))->Fill(d->CosThetaStarD0());
438 else ((TH1F*)fDistr->FindObject("hcosthetastarS"))->Fill(d->CosThetaStarD0bar());
439
440 ((TH1F*)fDistr->FindObject("hd0d0S"))->Fill(d->Prodd0d0());
441
442 ((TH1F*)fDistr->FindObject("hcosthetapointS"))->Fill(d->CosPointingAngle());
443 ((TH1F*)fDistr->FindObject("hcosthpointd0d0S"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
444
445 //cout<<"ok point"<<endl;
a41f6fad 446
feb73eca 447
ce39f0ac 448 } else{ //Background or LS
449 //cout<<"is background"<<endl;
450 AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(0);
451 if(!prong) cout<<"No daughter found";
452 else{
453 if(prong->Charge()==1) {fTotPosPairs[4]++;} else {fTotNegPairs[4]++;}
454 }
455 ((TH1F*)fDistr->FindObject("hptB"))->Fill(d->PtProng(0));
456 //cout<<"ptok"<<endl;
457 ((TH1F*)fDistr->FindObject("hd0B"))->Fill(d->Getd0Prong(0));
458 //cout<<"d0ok"<<endl;
459 ((TH1F*)fDistr->FindObject("hdcaB"))->Fill(d->GetDCA());
460 //cout<<"dcaok"<<endl;
461 ((TH1F*)fDistr->FindObject("hcosthetastarB"))->Fill(d->CosThetaStarD0());
462 ((TH1F*)fDistr->FindObject("hcosthetastarB"))->Fill(d->CosThetaStarD0bar());
463 ((TH1F*)fDistr->FindObject("hd0d0B"))->Fill(d->Prodd0d0());
464 //cout<<"d0d0ok"<<endl;
465 ((TH1F*)fDistr->FindObject("hcosthetapointB"))->Fill(d->CosPointingAngle());
466 ((TH1F*)fDistr->FindObject("hcosthpointd0d0B"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
467
468 //cout<<"pointok"<<endl;
469
470
471 }
472 } //inv mass cut
473 }
a41f6fad 474
9de8c723 475 //cuts order
476// printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
477// printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
478// printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
479// printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
480// printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
481// printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
482// printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
483// printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
484// printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
a41f6fad 485
49061176 486 Double_t pt = d->Pt();
a41f6fad 487
49061176 488 Int_t ptbin=0;
489
490 //cout<<"P_t = "<<pt<<endl;
491 if (pt>0. && pt<=1.) {
a4ae02cd 492 ptbin=1;
9de8c723 493 fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.5);
6ce5994f 494 fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.00025,0.7);
9de8c723 495// fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.7);
6ce5994f 496// fVHFloose->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,1,1,-0.00015,0.5);
49061176 497 //printf("I'm in the bin %d\n",ptbin);
a4ae02cd 498
49061176 499 }
46c96ce5 500
501 if(pt>1. && pt<=3.) {
9de8c723 502 if(pt>1. && pt<=2.) ptbin=2;
503 if(pt>2. && pt<=3.) ptbin=3;
504 fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0002,0.8);
46c96ce5 505 fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
506 //printf("I'm in the bin %d\n",ptbin);
507 }
508
509 if(pt>3. && pt<=5.){
9de8c723 510 ptbin=4;
511 fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0001,0.8);
6ce5994f 512 fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.8);
49061176 513 //printf("I'm in the bin %d\n",ptbin);
49061176 514 }
46c96ce5 515 if(pt>5.){
9de8c723 516 ptbin=5;
517 fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00005,0.8);
46c96ce5 518 fVHFloose->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.9);
519 }//if(pt>5)
520
a4ae02cd 521 //printf("I'm in the bin %d\n",ptbin);
522 //old
523 //fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed
9de8c723 524 if(prongsinacc){
525 FillHists(ptbin,d,mcArray,fVHFPPR,fOutputPPR);
526 FillHists(ptbin,d,mcArray,fVHFloose,fOutputloose);
527 }
49061176 528 if(unsetvtx) d->UnsetOwnPrimaryVtx();
49061176 529 }
46c96ce5 530
49061176 531
532
533
534 // Post the data
9de8c723 535 PostData(1,fOutputPPR);
a4ae02cd 536 PostData(2,fOutputloose);
a41f6fad 537 PostData(4,fDistr);
49061176 538 return;
539}
540//____________________________________________________________________________*
a4ae02cd 541void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAnalysisVertexingHF *vhf, TList *listout){
49061176 542 //
543 // function used in UserExec:
544 //
9de8c723 545
546
547 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
548
49061176 549 Int_t okD0=0,okD0bar=0;
feb73eca 550 //cout<<"inside FillHist"<<endl;
a4ae02cd 551 if(part->SelectD0(vhf->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
49061176 552 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
553 //printf("SELECTED\n");
feb73eca 554
9de8c723 555
feb73eca 556 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
557 if(!prong) cout<<"No daughter found";
558 else{
9de8c723 559 if(prong->Charge()==1) {fTotPosPairs[ptbin-1]++;} else {fTotNegPairs[ptbin-1]++;}
feb73eca 560 }
561
562 TString fillthis="";
c8112d2f 563 Int_t pdgDgD0toKpi[2]={321,211};
ce39f0ac 564 Int_t labD0=-1;
565 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)
a41f6fad 566
9de8c723 567 //count candidates selected by cuts
568 fNentries->Fill(2);
569 //count true D0 selected by cuts
570 if (labD0>=0) fNentries->Fill(3);
571 PostData(3,fNentries);
572
a4ae02cd 573 if (okD0==1) {
574 fillthis="histMass_";
7646d6da 575 fillthis+=ptbin;
576 //cout<<"Filling "<<fillthis<<endl;
577
49061176 578 //printf("Fill mass with D0");
a4ae02cd 579 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
feb73eca 580
a4ae02cd 581 if(labD0>=0) {
feb73eca 582 if(fArray==1) cout<<"LS signal ERROR"<<endl;
583
a4ae02cd 584 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
585 Int_t pdgD0 = partD0->GetPdgCode();
586 //cout<<"pdg = "<<pdgD0<<endl;
a4ae02cd 587 if (pdgD0==421){ //D0
588 //cout<<"Fill S with D0"<<endl;
46c96ce5 589 fillthis="histSgn_";
590 fillthis+=ptbin;
a4ae02cd 591 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
9de8c723 592 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
593 fillthis="histSgn27_";
594 fillthis+=ptbin;
595 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
596 }
a4ae02cd 597 } else{ //it was a D0bar
46c96ce5 598 fillthis="histRfl_";
599 fillthis+=ptbin;
600 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
a4ae02cd 601 }
602 } else {//background
603 fillthis="histBkg_";
604 fillthis+=ptbin;
605 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
606 }
feb73eca 607
49061176 608 }
609 if (okD0bar==1) {
a4ae02cd 610 fillthis="histMass_";
611 fillthis+=ptbin;
49061176 612 //printf("Fill mass with D0bar");
a4ae02cd 613 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
feb73eca 614
a4ae02cd 615 if(labD0>=0) {
feb73eca 616 if(fArray==1) cout<<"LS signal ERROR"<<endl;
a4ae02cd 617 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
618 Int_t pdgD0 = partD0->GetPdgCode();
619 //cout<<" pdg = "<<pdgD0<<endl;
620 if (pdgD0==-421){ //D0bar
621 fillthis="histSgn_";
622 fillthis+=ptbin;
623 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
9de8c723 624 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
625 fillthis="histSgn27_";
626 fillthis+=ptbin;
627 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
628 }
629
a4ae02cd 630
631 } else{
46c96ce5 632 fillthis="histRfl_";
633 fillthis+=ptbin;
634 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
a4ae02cd 635 }
feb73eca 636 } else {//background or LS
a4ae02cd 637 fillthis="histBkg_";
638 fillthis+=ptbin;
639 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
640 }
641 }
642
643
49061176 644 } //else cout<<"NOT SELECTED"<<endl;
645
646}
647//________________________________________________________________________
648void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
649{
650 // Terminate analysis
651 //
652 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
653
9de8c723 654 fOutputPPR = dynamic_cast<TList*> (GetOutputData(1));
655 if (!fOutputPPR) {
a41f6fad 656 printf("ERROR: fOutputthight not available\n");
a4ae02cd 657 return;
658 }
659 fOutputloose = dynamic_cast<TList*> (GetOutputData(2));
660 if (!fOutputloose) {
a41f6fad 661 printf("ERROR: fOutputloose not available\n");
662 return;
663 }
9de8c723 664 fDistr = dynamic_cast<TList*> (GetOutputData(4));
a41f6fad 665 if (!fDistr) {
666 printf("ERROR: fDistr not available\n");
49061176 667 return;
668 }
669
feb73eca 670 if(fArray==1){
671 for(Int_t ipt=0;ipt<4;ipt++){
672 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
673
674
675 if(fLsNormalization>0) {
9de8c723 676
feb73eca 677 TString massName="histMass_";
678 massName+=ipt+1;
9de8c723 679 ((TH1F*)fOutputPPR->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputPPR->FindObject(massName))->GetEntries());
feb73eca 680 ((TH1F*)fOutputloose->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputloose->FindObject(massName))->GetEntries());
681 }
682 }
683
684 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
685
686 if(fLsNormalization>0) {
687
688 TString nameDistr="hptB";
689 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
690 nameDistr="hdcaB";
691 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
692 nameDistr="hcosthetastarB";
693 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
694 nameDistr="hd0B";
695 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
696 nameDistr="hd0d0B";
697 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
698 nameDistr="hcosthetapointB";
699 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
700 nameDistr="hcosthpointd0d0B";
701 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
702
703 }
704 }
6ce5994f 705
49061176 706 return;
707}
708