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.);
322 namedistr="hdeclvS_";
324 TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
326 namedistr="hdeclvB_";
328 TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
330 namedistr="hnormdeclvS_";
332 TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
334 namedistr="hnormdeclvB_";
336 TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
339 namedistr="hcosthetapointS_";
341 TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
342 namedistr="hcosthetapointB_";
344 TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
346 // costhetapoint vs d0 or d0d0
347 namedistr="hcosthpointd0S_";
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);
351 namedistr="hcosthpointd0B_";
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);
355 namedistr="hcosthpointd0d0S_";
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_";
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);
366 fDistr->Add(hptpiSnoMcut);
367 fDistr->Add(hptKSnoMcut);
368 fDistr->Add(hptB1pnoMcut);
369 fDistr->Add(hptB2pnoMcut);
380 fDistr->Add(hd0vpiS);
383 fDistr->Add(hd0vp0B);
384 fDistr->Add(hd0vp1B);
389 fDistr->Add(hd0d0vS);
390 fDistr->Add(hd0d0vB);
392 fDistr->Add(hcosthetastarS);
393 fDistr->Add(hcosthetastarB);
395 fDistr->Add(hcosthetapointS);
396 fDistr->Add(hcosthetapointB);
398 fDistr->Add(hdeclengthS);
399 fDistr->Add(hdeclengthB);
401 fDistr->Add(hnormdeclengthS);
402 fDistr->Add(hnormdeclengthB);
404 fDistr->Add(hdeclengthvS);
405 fDistr->Add(hdeclengthvB);
407 fDistr->Add(hnormdeclengthvS);
408 fDistr->Add(hnormdeclengthvB);
410 fDistr->Add(hcosthpointd0S);
411 fDistr->Add(hcosthpointd0B);
413 fDistr->Add(hcosthpointd0d0S);
414 fDistr->Add(hcosthpointd0d0B);
417 //histograms of invariant mass distributions
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();
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();
430 //distribution w/o cuts
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();
434 tmpMB->SetName(nameMassNocutsB.Data());
438 //MC signal and background
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();
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();
449 //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
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();
454 // printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
456 fOutputMass->Add(tmpMt);
457 fOutputMass->Add(tmpSt);
458 fOutputMass->Add(tmpS27t);
459 fOutputMass->Add(tmpBt);
460 fOutputMass->Add(tmpRt);
468 //histograms for vertex checking and TOF checking
469 TString checkname="hptGoodTr";
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";
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");
479 TH1F* hTOFsig=new TH1F(checkname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
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.);
484 checkname="hTOFtime";
485 TH1F* hTOFtime=new TH1F(checkname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
488 fChecks->Add(hptGoodTr);
489 fChecks->Add(hdistrGoodTr);
490 fChecks->Add(hTOFsig);
491 fChecks->Add(hTOFtimeKaonHyptime);
492 fChecks->Add(hTOFtime);
494 const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
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.);
498 fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
499 fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
500 fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
501 fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtx");
502 fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
503 fNentries->GetXaxis()->SetBinLabel(6,"ptbin = -1");
504 fNentries->GetXaxis()->SetBinLabel(7,"no daughter");
505 fNentries->GetXaxis()->SetBinLabel(8,"nCandSel(Tr)");
507 fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
511 PostData(1,fOutputMass);
513 PostData(3,fNentries);
519 //________________________________________________________________________
520 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
522 // Execute analysis for current event:
523 // heavy flavor candidates association to MC truth
524 //cout<<"I'm in UserExec"<<endl;
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]);
539 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
542 if(fArray==0){ //D0 candidates
543 // load D0->Kpi candidates
544 //cout<<"D0 candidates"<<endl;
546 } else { //LikeSign candidates
547 // load like sign candidates
548 //cout<<"LS candidates"<<endl;
549 bname="LikeSign2Prong";
552 TClonesArray *inputArray=0;
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());
563 if(aodHandler->GetExtensions()) {
564 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
565 AliAODEvent* aodFromExt = ext->GetAOD();
566 inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
569 inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
574 printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
579 // fix for temporary bug in ESDfilter
580 // the AODs with null vertex pointer didn't pass the PhysSel
581 if(!aod->GetPrimaryVertex()) return;
583 TClonesArray *mcArray = 0;
584 AliAODMCHeader *mcHeader = 0;
588 mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
590 printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
595 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
597 printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
602 //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
604 //histogram filled with 1 for every AOD
609 // AOD primary vertex
610 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
612 AliError("There is no primary vertex !");
616 Bool_t isGoodVtx=kFALSE;
619 TString primTitle = vtx1->GetTitle();
620 if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
625 //cout<<"Start checks"<<endl;
626 Int_t ntracks=0,isGoodTrack=0;
628 if(aod) ntracks=aod->GetNTracks();
630 //cout<<"ntracks = "<<ntracks<<endl;
632 //loop on tracks in the event
633 for (Int_t k=0;k<ntracks;k++){
634 AliAODTrack* track=aod->GetTrack(k);
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();
645 if(!pid) {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
648 pid->GetIntegratedTimes(times);
650 ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
651 ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
653 ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
654 if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
657 //check clusters of the tracks
658 Int_t nclsTot=0,nclsSPD=0;
660 for(Int_t l=0;l<6;l++) {
661 if(TESTBIT(track->GetITSClusterMap(),l)) {
662 nclsTot++; if(l<2) nclsSPD++;
666 if (track->Pt()>0.3 &&
667 track->GetStatus()&AliESDtrack::kTPCrefit &&
668 track->GetStatus()&AliESDtrack::kITSrefit &&
670 nclsSPD>0) {//fill hist good tracks
672 ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
675 //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
676 ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
678 //number of events with good vertex and at least 2 good tracks
679 if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
682 // loop over candidates
683 Int_t nInD0toKpi = inputArray->GetEntriesFast();
684 if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
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);
692 if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
693 AliDebug(1,"at least one daughter not found!");
698 Bool_t unsetvtx=kFALSE;
699 if(!d->GetOwnPrimaryVtx()) {
700 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
705 //check reco daughter in acceptance
707 Double_t eta0=d->EtaProng(0);
708 Double_t eta1=d->EtaProng(1);
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
713 Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
714 if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(1))->GetTPCNcls() < 70) isSelected=kFALSE;
715 if (!isSelected) return;
717 if(fDebug>2) cout<<"tracks selected"<<endl;
719 FillVarHists(aod,d,mcArray,fCuts,fDistr);
720 FillMassHists(d,mcArray,fCuts,fOutputMass);
723 if(unsetvtx) d->UnsetOwnPrimaryVtx();
728 PostData(1,fOutputMass);
730 PostData(3,fNentries);
736 //____________________________________________________________________________
737 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
739 // function used in UserExec to fill variable histograms:
745 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
748 Double_t minvD0 = part->InvMassD0();
751 Double_t minvD0bar = part->InvMassD0bar();
752 //cout<<"inside mass cut"<<endl;
754 Int_t pdgDgD0toKpi[2]={321,211};
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;
762 isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate);
764 //cout<<"Not Selected"<<endl;
769 Double_t invmasscut=0.03;
771 TString fillthispi="",fillthisK="",fillthis="";
773 Int_t ptbin=cuts->PtBin(part->Pt());
774 if(ptbin==-1) {fNentries->Fill(5); return;} //out of bounds
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;
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);
793 Double_t d0[2]={dz1[0],dz2[0]};
795 if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed
796 //printf("\nif no cuts or cuts passed\n");
797 if(lab>=0 && fReadMC){ //signal
799 //check pdg of the prongs
800 AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
801 AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
802 if(!prong0 || !prong1) {
806 labprong[0]=prong0->GetLabel();
807 labprong[1]=prong1->GetLabel();
808 AliAODMCParticle *mcprong=0;
809 Int_t pdgProng[2]={0,0};
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();
816 //no mass cut ditributions: ptbis
818 fillthispi="hptpiSnoMcut_";
821 fillthisK="hptKSnoMcut_";
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));
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));
834 //no mass cut ditributions: mass
838 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421){//D0
839 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
842 ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
845 //apply cut on invariant mass on the pair
846 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
848 if(fArray==1) cout<<"LS signal: ERROR"<<endl;
849 for (Int_t iprong=0; iprong<2; iprong++){
850 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
851 labprong[iprong]=prong->GetLabel();
853 //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
854 if(labprong[iprong]>=0) mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
855 Int_t pdgprong=mcprong->GetPdgCode();
857 if(TMath::Abs(pdgprong)==211) {
860 fillthispi="hptpiS_";
862 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
863 fillthispi="hd0piS_";
865 ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
866 if(d0[iprong] != -99) {
868 fillthispi="hd0vpiS_";
870 ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
875 if(TMath::Abs(pdgprong)==321) {
876 //cout<<"kappa"<<endl;
880 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
883 ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
884 if (d0[iprong] != -99){
887 ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
891 fillthis="hcosthpointd0S_";
893 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
895 } //end loop on prongs
899 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
901 fillthis="hcosthetapointS_";
903 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
905 fillthis="hcosthpointd0d0S_";
907 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
909 fillthis="hcosthetastarS_";
911 if (((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
912 else ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
916 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
918 if(d0[0] != -99 && d0[1] != -99 ){
921 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
926 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
928 fillthis="hnormdeclS_";
930 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
935 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
937 fillthis="hnormdeclvS_";
939 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
943 } else{ //Background or LS
944 //cout<<"is background"<<endl;
946 //no mass cut distributions: mass, ptbis
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);
953 fillthis="hptB1prongnoMcut_";
956 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
958 fillthis="hptB2prongsnoMcut_";
960 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
961 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
963 //apply cut on invariant mass on the pair
964 if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
967 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
969 if(fDebug>2) cout<<"No daughter found";
973 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
976 //normalise pt distr to half afterwards
979 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
980 ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
984 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
987 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
990 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
991 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
995 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
998 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1002 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1003 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1006 fillthis="hcosthpointd0B_";
1008 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
1009 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
1014 ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1016 fillthis="hcosthetastarB_";
1018 if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
1019 if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
1023 ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1025 if(d0[0] != -99 && d0[1] != -99 ){
1026 fillthis="hd0d0vB_";
1028 ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1031 fillthis="hcosthetapointB_";
1033 ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1035 fillthis="hcosthpointd0d0B_";
1037 ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
1041 ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1043 fillthis="hnormdeclB_";
1045 ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
1048 fillthis="hdeclvB_";
1050 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::DecayLength(vtxProp));
1052 fillthis="hnormdeclvB_";
1054 ((TH1F*)listout->FindObject(fillthis))->Fill(part->AliAODRecoDecay::NormalizedDecayLength(vtxProp));
1057 }//else (background)
1061 //____________________________________________________________________________
1063 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
1065 // function used in UserExec to fill mass histograms:
1069 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1071 Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1072 //cout<<"is selected = "<<isSelected<<endl;
1074 //cout<<"check cuts = "<<endl;
1077 //cout<<"Not Selected"<<endl;
1081 if(fDebug>2) cout<<"Candidate selected"<<endl;
1083 Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1084 //printf("SELECTED\n");
1085 Int_t ptbin=cuts->PtBin(part->Pt());
1087 AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
1089 AliDebug(2,"No daughter found");
1093 if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
1096 TString fillthis="";
1097 Int_t pdgDgD0toKpi[2]={321,211};
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)
1101 //count candidates selected by cuts
1103 //count true D0 selected by cuts
1104 if (fReadMC && labD0>=0) fNentries->Fill(2);
1105 //PostData(3,fNentries);
1107 if (isSelected==1 || isSelected==3) { //D0
1108 fillthis="histMass_";
1110 //cout<<"Filling "<<fillthis<<endl;
1112 //printf("Fill mass with D0");
1113 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1116 if(fArray==1) cout<<"LS signal ERROR"<<endl;
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_";
1125 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1126 if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1127 fillthis="histSgn27_";
1129 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1131 } else{ //it was a D0bar
1132 fillthis="histRfl_";
1134 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1136 } else {//background
1137 fillthis="histBkg_";
1139 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1143 if (isSelected>1) { //D0bar
1144 fillthis="histMass_";
1146 //printf("Fill mass with D0bar");
1147 ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
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_";
1157 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1158 if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1159 fillthis="histSgn27_";
1161 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1166 fillthis="histRfl_";
1168 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1170 } else {//background or LS
1171 fillthis="histBkg_";
1173 ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1180 //__________________________________________________________________________
1181 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1182 //Calculate the primary vertex w/o the daughter tracks of the candidate
1184 AliESDVertex *vertexESD=0x0;
1185 AliAODVertex *vertexAOD=0x0;
1186 AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1189 Int_t nTrksToSkip=2;
1190 AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1192 AliDebug(2,"no daughter found!");
1195 skipped[0]=dgTrack->GetID();
1196 dgTrack = (AliAODTrack*)d->GetDaughter(1);
1198 AliDebug(2,"no daughter found!");
1201 skipped[1]=dgTrack->GetID();
1205 vertexer->SetSkipTracks(nTrksToSkip,skipped);
1206 vertexer->SetMinClusters(4);
1207 vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
1208 if(!vertexESD) return vertexAOD;
1209 if(vertexESD->GetNContributors()<=0) {
1210 AliDebug(2,"vertexing failed");
1211 delete vertexESD; vertexESD=NULL;
1215 delete vertexer; vertexer=NULL;
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;
1225 vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1231 //________________________________________________________________________
1232 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1234 // Terminate analysis
1236 if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1239 fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1241 printf("ERROR: fOutputMass not available\n");
1244 fDistr = dynamic_cast<TList*> (GetOutputData(2));
1246 printf("ERROR: fDistr not available\n");
1250 fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
1253 printf("ERROR: fNEntries not available\n");
1257 fChecks = dynamic_cast<TList*> (GetOutputData(4));
1259 printf("ERROR: fChecks not available\n");
1264 for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1268 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]);
1271 if(fLsNormalization>1e-6) {
1273 TString massName="histMass_";
1275 ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1280 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1282 if(fLsNormalization>1e-6) {
1284 TString nameDistr="hptB_";
1286 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1289 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1290 nameDistr="hcosthetastarB_";
1292 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1295 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1296 nameDistr="hd0d0B_";
1298 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1299 nameDistr="hcosthetapointB_";
1301 ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1302 nameDistr="hcosthpointd0d0B_";
1304 ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1309 TString cvname,cstname;
1319 TCanvas *cMass=new TCanvas(cvname,cvname);
1321 ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
1323 TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
1326 fNentries->Draw("htext0");
1328 // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));