1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
19 // and comparison with the MC truth and cut variables distributions.
21 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
22 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
23 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
24 /////////////////////////////////////////////////////////////
26 #include <Riostream.h>
27 #include <TClonesArray.h>
33 #include <TDatabasePDG.h>
35 #include <AliAnalysisDataSlot.h>
36 #include <AliAnalysisDataContainer.h>
37 #include "AliAnalysisManager.h"
38 #include "AliESDtrack.h"
39 #include "AliVertexerTracks.h"
40 #include "AliAODHandler.h"
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"
53 ClassImp(AliAnalysisTaskSED0Mass)
56 //________________________________________________________________________
57 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
74 // Default constructor
77 //________________________________________________________________________
78 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
79 AliAnalysisTaskSE(name),
95 // Default constructor
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;}
104 // Output slot #1 writes into a TList container (mass with cuts)
105 DefineOutput(1,TList::Class()); //My private output
106 // Output slot #2 writes into a TList container (distributions)
107 DefineOutput(2,TList::Class()); //My private output
108 // Output slot #3 writes into a TH1F container (number of events)
109 DefineOutput(3,TH1F::Class()); //My private output
110 // Output slot #4 writes into a TList container (quality check)
111 DefineOutput(4,TList::Class()); //My private output
112 // Output slot #5 writes into a TList container (cuts)
113 DefineOutput(5,AliRDHFCutsD0toKpi::Class()); //My private output
117 //________________________________________________________________________
118 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
144 //________________________________________________________________________
145 void AliAnalysisTaskSED0Mass::Init()
149 if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
152 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
153 const char* nameoutput=GetOutputSlot(5)->GetContainer()->GetName();
154 copyfCuts->SetName(nameoutput);
156 PostData(5,copyfCuts);
162 //________________________________________________________________________
163 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
166 // Create the output container
168 if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
170 // Several histograms are more conveniently managed in a TList
171 fOutputMass = new TList();
172 fOutputMass->SetOwner();
173 fOutputMass->SetName("listMass");
175 fDistr = new TList();
177 fDistr->SetName("distributionslist");
179 fChecks = new TList();
181 fChecks->SetName("checkHistograms");
183 TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
185 for(Int_t i=0;i<fCuts->GetNPtBins();i++){
187 nameMass="histMass_";
189 nameSgn27="histSgn27_";
197 nameMassNocutsS="hMassS_";
199 nameMassNocutsB="hMassB_";
202 //histograms of cut variable distributions
207 TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
211 TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
215 TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
218 namedistr="hptpiSnoMcut_";
220 TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
222 namedistr="hptKSnoMcut_";
224 TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
226 namedistr="hptB1prongnoMcut_";
228 TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
230 namedistr="hptB2prongsnoMcut_";
232 TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
238 TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
241 TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
244 namedistr="hcosthetastarS_";
246 TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
247 namedistr="hcosthetastarB_";
249 TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
254 TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
258 TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
262 TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
266 TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
270 TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
272 namedistr="hd0vpiS_";
274 TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
278 TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
281 TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
283 namedistr="hd0vp0B_";
285 TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
287 namedistr="hd0vp1B_";
289 TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
293 TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
296 TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
298 namedistr="hd0d0vS_";
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_";
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);
308 TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.6);
312 TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.6);
314 namedistr="hnormdeclS_";
316 TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,10.);
318 namedistr="hnormdeclB_";
320 TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,10.);
323 namedistr="hcosthetapointS_";
325 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
326 namedistr="hcosthetapointB_";
328 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
330 // costhetapoint vs d0 or d0d0
331 namedistr="hcosthpointd0S_";
333 TH2F *hcosthpointd0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
335 namedistr="hcosthpointd0B_";
337 TH2F *hcosthpointd0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
339 namedistr="hcosthpointd0d0S_";
341 TH2F *hcosthpointd0d0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
342 namedistr="hcosthpointd0d0B_";
344 TH2F *hcosthpointd0d0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
350 fDistr->Add(hptpiSnoMcut);
351 fDistr->Add(hptKSnoMcut);
352 fDistr->Add(hptB1pnoMcut);
353 fDistr->Add(hptB2pnoMcut);
364 fDistr->Add(hd0vpiS);
367 fDistr->Add(hd0vp0B);
368 fDistr->Add(hd0vp1B);
373 fDistr->Add(hd0d0vS);
374 fDistr->Add(hd0d0vB);
376 fDistr->Add(hcosthetastarS);
377 fDistr->Add(hcosthetastarB);
379 fDistr->Add(hcosthetapointS);
380 fDistr->Add(hcosthetapointB);
382 fDistr->Add(hdeclengthS);
383 fDistr->Add(hdeclengthB);
385 fDistr->Add(hnormdeclengthS);
386 fDistr->Add(hnormdeclengthB);
388 fDistr->Add(hcosthpointd0S);
389 fDistr->Add(hcosthpointd0B);
391 fDistr->Add(hcosthpointd0d0S);
392 fDistr->Add(hcosthpointd0d0B);
395 //histograms of invariant mass distributions
397 TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
398 TH1F *tmpMl=(TH1F*)tmpMt->Clone();
402 //to compare with AliAnalysisTaskCharmFraction
403 TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.765,1.965);
404 TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
408 //distribution w/o cuts
409 // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
410 TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
411 TH1F *tmpMB=(TH1F*)tmpMS->Clone();
412 tmpMB->SetName(nameMassNocutsB.Data());
416 //MC signal and background
417 TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
418 TH1F *tmpSl=(TH1F*)tmpSt->Clone();
422 TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
423 TH1F *tmpBl=(TH1F*)tmpBt->Clone();
427 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
428 TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
429 TH1F *tmpRl=(TH1F*)tmpRt->Clone();
432 // printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
434 fOutputMass->Add(tmpMt);
435 fOutputMass->Add(tmpSt);
436 fOutputMass->Add(tmpS27t);
437 fOutputMass->Add(tmpBt);
438 fOutputMass->Add(tmpRt);
446 //histograms for vertex checking and TOF checking
447 TString checkname="hptGoodTr";
449 TH1F* hptGoodTr=new TH1F(checkname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Number",200,0.,8.);
450 hptGoodTr->SetTitleOffset(1.3,"Y");
451 checkname="hdistrGoodTr";
453 TH1F* hdistrGoodTr=new TH1F(checkname.Data(),"Distribution of number of good tracks per event;no.good-tracks/ev;Entries",31,0,31);
454 hdistrGoodTr->SetTitleOffset(1.3,"Y");
457 TH1F* hTOFsig=new TH1F(checkname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
459 checkname="hTOFtimeKaonHyptime";
460 TH2F* hTOFtimeKaonHyptime=new TH2F(checkname.Data(),"TOFtime - timeHypothesisForKaon;p_{t}[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
462 checkname="hTOFtime";
463 TH1F* hTOFtime=new TH1F(checkname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
466 fChecks->Add(hptGoodTr);
467 fChecks->Add(hdistrGoodTr);
468 fChecks->Add(hTOFsig);
469 fChecks->Add(hTOFtimeKaonHyptime);
470 fChecks->Add(hTOFtime);
472 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
474 fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex *** Integral(5,6) = pt out of bounds", 8,0.,8.);
476 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
477 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
478 fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
479 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtx");
480 fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
481 fNentries->GetXaxis()->SetBinLabel(6,"ptbin = -1");
482 fNentries->GetXaxis()->SetBinLabel(7,"no daughter");
483 fNentries->GetXaxis()->SetBinLabel(8,"nCandSel(Tr)");
485 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
489 PostData(1,fOutputMass);
491 PostData(3,fNentries);
497 //________________________________________________________________________
498 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
500 // Execute analysis for current event:
501 // heavy flavor candidates association to MC truth
502 //cout<<"I'm in UserExec"<<endl;
506 // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
507 // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
508 // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
509 // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
510 // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
511 // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
512 // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
513 // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
514 // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
517 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
520 if(fArray==0){ //D0 candidates
521 // load D0->Kpi candidates
522 //cout<<"D0 candidates"<<endl;
524 } else { //LikeSign candidates
525 // load like sign candidates
526 //cout<<"LS candidates"<<endl;
527 bname="LikeSign2Prong";
530 TClonesArray *inputArray=0;
532 if(!aod && AODEvent() && IsStandardAOD()) {
533 // In case there is an AOD handler writing a standard AOD, use the AOD
534 // event in memory rather than the input (ESD) event.
535 aod = dynamic_cast<AliAODEvent*> (AODEvent());
536 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
537 // have to taken from the AOD event hold by the AliAODExtension
538 AliAODHandler* aodHandler = (AliAODHandler*)
539 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
541 if(aodHandler->GetExtensions()) {
542 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
543 AliAODEvent* aodFromExt = ext->GetAOD();
544 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
547 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
552 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
557 TClonesArray *mcArray = 0;
558 AliAODMCHeader *mcHeader = 0;
562 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
564 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
569 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
571 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
576 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
578 //histogram filled with 1 for every AOD
583 // AOD primary vertex
584 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
586 Bool_t isGoodVtx=kFALSE;
589 TString primTitle = vtx1->GetTitle();
590 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
595 //cout<<"Start checks"<<endl;
596 Int_t ntracks=0,isGoodTrack=0;
598 if(aod) ntracks=aod->GetNTracks();
600 //cout<<"ntracks = "<<ntracks<<endl;
602 //loop on tracks in the event
603 for (Int_t k=0;k<ntracks;k++){
604 AliAODTrack* track=aod->GetTrack(k);
608 if(!(track->GetStatus()&AliESDtrack::kTPCrefit &&
609 track->GetStatus()&AliESDtrack::kITSrefit &&
610 track->GetTPCNcls() >=70 &&
611 track->GetStatus()&AliESDtrack::kTOFpid &&
612 track->GetStatus()&AliESDtrack::kTOFout &&
613 track->GetStatus()&AliESDtrack::kTIME)) continue;
614 AliAODPid *pid = track->GetDetPid();
615 if(!pid) {cout<<"No AliAODPid found"<<endl; continue;}
618 pid->GetIntegratedTimes(times);
620 ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
621 ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
623 ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
624 if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
627 //check clusters of the tracks
628 Int_t nclsTot=0,nclsSPD=0;
630 for(Int_t l=0;l<6;l++) {
631 if(TESTBIT(track->GetITSClusterMap(),l)) {
632 nclsTot++; if(l<2) nclsSPD++;
636 if (track->Pt()>0.3 &&
637 track->GetStatus()&AliESDtrack::kTPCrefit &&
638 track->GetStatus()&AliESDtrack::kITSrefit &&
640 nclsSPD>0) {//fill hist good tracks
641 //cout<<"in if"<<endl;
642 ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
645 //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
646 ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
648 //number of events with good vertex and at least 2 good tracks
649 if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
652 // loop over candidates
653 Int_t nInD0toKpi = inputArray->GetEntriesFast();
654 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
656 for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
657 //Int_t nPosPairs=0, nNegPairs=0;
658 //cout<<"inside the loop"<<endl;
659 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
662 if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
663 AliDebug(1,"at least one daughter not found!");
668 Bool_t unsetvtx=kFALSE;
669 if(!d->GetOwnPrimaryVtx()) {
670 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
675 //check reco daughter in acceptance
677 Double_t eta0=d->EtaProng(0);
678 Double_t eta1=d->EtaProng(1);
680 if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
681 //if( TMath::Abs(eta0)<0.9 && TMath::Abs(eta1)<0.9 ){
682 //apply cuts on tracks
683 Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
684 if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70) isSelected=kFALSE;
685 if (!isSelected) return;
687 if(fDebug>2) cout<<"tracks selected"<<endl;
689 FillVarHists(aod,d,mcArray,fCuts,fDistr);
690 FillMassHists(d,mcArray,fCuts,fOutputMass);
693 if(unsetvtx) d->UnsetOwnPrimaryVtx();
698 PostData(1,fOutputMass);
700 PostData(3,fNentries);
706 //____________________________________________________________________________
707 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
709 // function used in UserExec to fill variable histograms:
715 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
718 Double_t minvD0 = part->InvMassD0();
721 Double_t minvD0bar = part->InvMassD0bar();
722 //cout<<"inside mass cut"<<endl;
724 Int_t pdgDgD0toKpi[2]={321,211};
726 if(fReadMC) lab=part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
727 //Double_t pt = d->Pt(); //mother pt
728 Bool_t isSelected=kFALSE;
732 isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate);
734 //cout<<"Not Selected"<<endl;
739 Double_t invmasscut=0.03;
741 TString fillthispi="",fillthisK="",fillthis="";
743 Int_t ptbin=cuts->PtBin(part->Pt());
744 if(ptbin==-1) {fNentries->Fill(5); return;} //out of bounds
747 AliAODVertex *vtxProp=0x0;
748 vtxProp=GetPrimaryVtxSkipped(aod,part);
749 Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
750 dz1[0]=-99; dz2[0]=-99;
752 part->SetOwnPrimaryVtx(vtxProp);
753 //Bool_t unsetvtx=kTRUE;
754 //Calculate d0 for daughter tracks with Vtx Skipped
755 AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)part->GetDaughter(0));
756 AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)part->GetDaughter(1));
757 esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
758 esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
763 Double_t d0[2]={dz1[0],dz2[0]};
765 if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
766 //printf("\nif no cuts or cuts passed\n");
767 if(lab>=0 && fReadMC){ //signal
769 //check pdg of the prongs
770 AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
771 AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
772 if(!prong0 || !prong1) {
776 labprong[0]=prong0->GetLabel();
777 labprong[1]=prong1->GetLabel();
778 AliAODMCParticle *mcprong=0;
779 Int_t pdgProng[2]={0,0};
781 for (Int_t iprong=0;iprong<2;iprong++){
782 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
783 pdgProng[iprong]=mcprong->GetPdgCode();
786 //no mass cut ditributions: ptbis
788 fillthispi="hptpiSnoMcut_";
791 fillthisK="hptKSnoMcut_";
794 if (TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321){
795 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
796 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
798 if (TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211){
799 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
800 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
804 //no mass cut ditributions: mass
808 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
809 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
812 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
815 //apply cut on invariant mass on the pair
816 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
818 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
819 for (Int_t iprong=0; iprong<2; iprong++){
820 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
821 labprong[iprong]=prong->GetLabel();
823 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
824 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
825 Int_t pdgprong=mcprong->GetPdgCode();
827 if(TMath::Abs(pdgprong)==211) {
830 fillthispi="hptpiS_";
832 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
833 fillthispi="hd0piS_";
835 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
836 if(d0[iprong] != -99) {
838 fillthispi="hd0vpiS_";
840 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
845 if(TMath::Abs(pdgprong)==321) {
846 //cout<<"kappa"<<endl;
850 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
853 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
854 if (d0[iprong] != -99){
857 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
861 fillthis="hcosthpointd0S_";
863 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
865 } //end loop on prongs
869 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
871 fillthis="hcosthetapointS_";
873 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
875 fillthis="hcosthpointd0d0S_";
877 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
879 fillthis="hcosthetastarS_";
881 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
882 else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
886 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
888 if(d0[0] != -99 && d0[1] != -99 ){
891 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
896 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
898 fillthis="hnormdeclS_";
900 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
904 } else{ //Background or LS
905 //cout<<"is background"<<endl;
907 //no mass cut distributions: mass, ptbis
911 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
912 if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
914 fillthis="hptB1prongnoMcut_";
917 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
919 fillthis="hptB2prongsnoMcut_";
921 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
922 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
924 //apply cut on invariant mass on the pair
925 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
928 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
930 if(fDebug>2) cout<<"No daughter found";
934 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
937 //normalise pt distr to half afterwards
940 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
941 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
945 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
948 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
952 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
955 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
957 fillthis="hcosthpointd0B_";
959 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
960 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
965 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
967 fillthis="hcosthetastarB_";
969 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
970 if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
974 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
976 fillthis="hcosthetapointB_";
978 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
980 fillthis="hcosthpointd0d0B_";
982 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
989 //____________________________________________________________________________
991 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
993 // function used in UserExec to fill mass histograms:
997 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
999 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1000 //cout<<"is selected = "<<isSelected<<endl;
1002 //cout<<"check cuts = "<<endl;
1005 //cout<<"Not Selected"<<endl;
1009 if(fDebug>2) cout<<"Candidate selected"<<endl;
1011 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1012 //printf("SELECTED\n");
1013 Int_t ptbin=cuts->PtBin(part->Pt());
1015 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
1017 AliDebug(2,"No daughter found");
1021 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
1024 TString fillthis="";
1025 Int_t pdgDgD0toKpi[2]={321,211};
1027 if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
1029 //count candidates selected by cuts
1031 //count true D0 selected by cuts
1032 if (fReadMC && labD0>=0) fNentries->Fill(2);
1033 //PostData(3,fNentries);
1035 if (isSelected==1 || isSelected==3) { //D0
1036 fillthis="histMass_";
1038 //cout<<"Filling "<<fillthis<<endl;
1040 //printf("Fill mass with D0");
1041 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1044 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1046 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1047 Int_t pdgD0 = partD0->GetPdgCode();
1048 //cout<<"pdg = "<<pdgD0<<endl;
1049 if (pdgD0==421){ //D0
1050 //cout<<"Fill S with D0"<<endl;
1051 fillthis="histSgn_";
1053 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1054 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1055 fillthis="histSgn27_";
1057 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1059 } else{ //it was a D0bar
1060 fillthis="histRfl_";
1062 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1064 } else {//background
1065 fillthis="histBkg_";
1067 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1071 if (isSelected>1) { //D0bar
1072 fillthis="histMass_";
1074 //printf("Fill mass with D0bar");
1075 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1078 if(fArray==1) cout<<"LS signal ERROR"<<endl;
1079 AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1080 Int_t pdgD0 = partD0->GetPdgCode();
1081 //cout<<" pdg = "<<pdgD0<<endl;
1082 if (pdgD0==-421){ //D0bar
1083 fillthis="histSgn_";
1085 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1086 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1087 fillthis="histSgn27_";
1089 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1094 fillthis="histRfl_";
1096 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1098 } else {//background or LS
1099 fillthis="histBkg_";
1101 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1108 //__________________________________________________________________________
1109 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1110 //Calculate the primary vertex w/o the daughter tracks of the candidate
1112 AliESDVertex *vertexESD=0x0;
1113 AliAODVertex *vertexAOD=0x0;
1114 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1117 Int_t nTrksToSkip=2;
1118 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1120 AliDebug(2,"no daughter found!");
1123 skipped[0]=dgTrack->GetID();
1124 dgTrack = (AliAODTrack*)d->GetDaughter(1);
1126 AliDebug(2,"no daughter found!");
1129 skipped[1]=dgTrack->GetID();
1133 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1134 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
1135 vertexer->SetMinClusters(4);
1136 if(!vertexESD) return vertexAOD;
1137 if(vertexESD->GetNContributors()<=0) {
1138 AliDebug(2,"vertexing failed");
1139 delete vertexESD; vertexESD=NULL;
1143 delete vertexer; vertexer=NULL;
1146 // convert to AliAODVertex
1147 Double_t pos[3],cov[6],chi2perNDF;
1148 vertexESD->GetXYZ(pos); // position
1149 vertexESD->GetCovMatrix(cov); //covariance matrix
1150 chi2perNDF = vertexESD->GetChi2toNDF();
1151 delete vertexESD; vertexESD=NULL;
1153 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1159 //________________________________________________________________________
1160 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1162 // Terminate analysis
1164 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1167 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1169 printf("ERROR: fOutputMass not available\n");
1172 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1174 printf("ERROR: fDistr not available\n");
1178 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
1181 printf("ERROR: fNEntries not available\n");
1185 fChecks = dynamic_cast<TList*> (GetOutputData(4));
1187 printf("ERROR: fChecks not available\n");
1192 for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1195 TString d0Bname=Form("hd0B_%d",ipt),d0p0Bname=Form("hd0p0B_%d",ipt),d0p1Bname=Form("hd0p1B_%d",ipt);
1196 ((TH1F*)fDistr->FindObject(d0Bname))->Add( ((TH1F*)fDistr->FindObject(d0p0Bname)), ((TH1F*)fDistr->FindObject(d0p1Bname)) );
1200 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1203 if(fLsNormalization>1e-6) {
1205 TString massName="histMass_";
1207 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1212 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1214 if(fLsNormalization>1e-6) {
1216 TString nameDistr="hptB_";
1218 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1221 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1222 nameDistr="hcosthetastarB_";
1224 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1227 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1228 nameDistr="hd0d0B_";
1230 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1231 nameDistr="hcosthetapointB_";
1233 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1234 nameDistr="hcosthpointd0d0B_";
1236 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1241 TString cvname,cstname;
1251 TCanvas *cMass=new TCanvas(cvname,cvname);
1253 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
1255 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
1258 fNentries->Draw("htext0");
1260 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));