[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / ReadAODVertexingHF.C

void ReadAODVertexingHF(const char *aodFileName="AliAOD.root",
			const char *aodHFFileName="AliAOD.VertexingHF.root") 
{
  //
  // Example macro to read D0->Kpi candidates from AOD (having the
  // standard AOD + a friend heavy-flavour AOD) and apply cuts
  // Origin: A.Dainese
  //

  Bool_t useParFiles=kFALSE;
  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
  LoadLibraries(useParFiles);

  // create a test histogram
  TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
  hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
  hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
  TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
  hMass->SetXTitle("Invariant mass [GeV]");
  hMass->SetYTitle("Entries");
  TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
  hSecVtxZ->SetXTitle("z of decay vertex [cm]");
  hSecVtxZ->SetYTitle("Entries");
  TH1F *hDeltaMassDstar = new TH1F("hDeltaMassDstar","D* delta mass plot",100,0,0.3);
  hDeltaMassDstar->SetXTitle("M(Kpipi)-M(Kpi) [GeV]");
  hDeltaMassDstar->SetYTitle("Entries");


  // open input file and get the TTree
  TFile inFile(aodFileName,"READ");
  if (!inFile.IsOpen()) return;

  TTree *aodTree = (TTree*)inFile.Get("aodTree");
  aodTree->AddFriend("aodTree",aodHFFileName);

  AliAODEvent *aod = new AliAODEvent();

  aod->ReadFromTree(aodTree);

  // load heavy flavour vertices
  TClonesArray *arrayVerticesHF = 
    (TClonesArray*)aod->GetList()->FindObject("VerticesHF"); 
  
  // load D0->Kpi candidates
  TClonesArray *arrayD0toKpi = 
    (TClonesArray*)aod->GetList()->FindObject("D0toKpi"); 
     
  // load 3prong candidates
  TClonesArray *array3Prong = 
    (TClonesArray*)aod->GetList()->FindObject("Charm3Prong"); 
    
  // load D* candidates
  TClonesArray *arrayDstar = 
    (TClonesArray*)aod->GetList()->FindObject("Dstar"); 
    
  // load cascade (V0+track) candidates
  TClonesArray *arrayCascades = 
    (TClonesArray*)aod->GetList()->FindObject("CascadesHF"); 
    

  Double_t cutsD0[9]=
  // cutsD0[0] = inv. mass half width [GeV]
  // cutsD0[1] = dca [cm]
  // cutsD0[2] = cosThetaStar
  // cutsD0[3] = pTK [GeV/c]
  // cutsD0[4] = pTPi [GeV/c]
  // cutsD0[5] = d0K [cm]   upper limit!
  // cutsD0[6] = d0Pi [cm]  upper limit!
  // cutsD0[7] = d0d0 [cm^2]
  // cutsD0[8] = cosThetaPoint
                     {1000.,
		      100000.,
		      1.1,
		      0.,
		      0.,
		      100000.,
		      100000.,
		      100000000.,
		      -1.1}; 
  Double_t cutsDstar[5]=
    // (to be passed to AliAODRecoCascadeHF::SelectDstar())
    // 0 = inv. mass half width of D* [GeV]
    // 1 = half width of (M_Kpipi-M_Kpi) [GeV]
    // 2 = PtMin of pi_s [GeV/c]
    // 3 = PtMax of pi_s [GeV/c]
    // 4 = theta, angle between the trace of pi_s and D0 decay plane [rad]
                     {999999.,
		      999999.,
		      0.1, 
		      1.0, 
		      0.5};

  Double_t cutsLctoV0[9]=// cuts on Lambdac candidates to V0+bachelor
                        // (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
                        // 0 = inv. mass half width in K0s hypothesis [GeV]   
                        // 1 = inv. mass half width in Lambda hypothesis [GeV]   
                        // 2 = inv. mass V0 in K0s hypothesis half width [GeV]   
                        // 3 = inv. mass V0 in Lambda hypothesis half width [GeV]   
                        // 4 = pT min Bachelor track [GeV/c]
                        // 5 = pT min V0-Positive track [GeV/c]
                        // 6 = pT min V0-Negative track [GeV/c]
                        // 7 = dca cut on the V0 (cm)
                        // 8 = dca cut on the cascade (cm)
    {2.0,2.0,1.0,1.0,0.0,0.0,0.0,1000.,1000.};

  Int_t nTotHF=0,nTotD0toKpi=0,nTotDstar=0,nTot3Prong=0,nTotCasc=0;
  AliAODVertex *vtx1=0;

  // loop over events
  Int_t nEvents = aodTree->GetEntries();
  for (Int_t nEv = 0; nEv < nEvents; nEv++) {
    cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;

    // read event
    aodTree->GetEvent(nEv);
    //aodTree->BranchRef();

    //print event info
    aod->GetHeader()->Print();

    // primary vertex
    vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
    vtx1->Print();

    // make trkIDtoEntry register (temporary)
    Int_t trkIDtoEntry[100000];
    for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
      AliAODTrack *track = aod->GetTrack(it);
      trkIDtoEntry[track->GetID()]=it;
    }
    

    // loop over D0->Kpi candidates
    Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
    nTotD0toKpi += nD0toKpi;
    cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
    
    for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
      AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);

      Bool_t unsetvtx=kFALSE;
      if(!d->GetOwnPrimaryVtx()) {
	d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	unsetvtx=kTRUE;
      }
      Int_t okD0=0,okD0bar=0; 
      if(d->SelectD0(cutsD0,okD0,okD0bar)) {
	//cout<<1e8*d->Prodd0d0()<<endl;
	hMass->Fill(d->InvMassD0(),0.5);
	hMass->Fill(d->InvMassD0bar(),0.5);
	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
	hSecVtxZ->Fill(d->GetSecVtxZ());
	//cout<<d->GetSecVtxX() <<endl;

	// get daughter AOD tracks
	AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
	AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
	if(!trk0 || !trk1) {
	  trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
	  trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
	  cout<<"references to standard AOD not available"<<endl;
	}
	cout<<"pt of positive track: "<<trk0->Pt()<<endl;

	// make a AliNeutralTrackParam from the D0 
	// and calculate impact parameters to primary vertex
	AliNeutralTrackParam trackD0(d);
	cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
	//trackD0.Print();
	Double_t dz[2],covdz[3];
	trackD0.PropagateToDCA(vtx1,aod->GetMagneticField(),1000.,dz,covdz);
	cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
      }
      if(unsetvtx) d->UnsetOwnPrimaryVtx();
    }


    // loop over D* candidates
    Int_t nDstar = arrayDstar->GetEntriesFast();
    nTotDstar += nDstar;
    cout<<"Number of D*->D0pi: "<<nDstar<<endl;
    
    for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
      AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar);
      Bool_t unsetvtx=kFALSE;
      if(!c->GetOwnPrimaryVtx()) {
	c->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
	c->Get2Prong()->SetOwnPrimaryVtx(vtx1);
	unsetvtx=kTRUE;
      }
      if(c->SelectDstar(cutsDstar,cutsD0)) {
	hDeltaMassDstar->Fill(c->DeltaInvMass());
	// get daughters
	AliAODTrack *trk = (AliAODTrack*)c->GetBachelor();
	AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)c->Get2Prong();
	cout<<"pt of soft pi: "<<trk->Pt()<<endl;
	cout<<"pt of D0: "<<d->Pt()<<endl;
      }

      if(unsetvtx) {
	c->UnsetOwnPrimaryVtx();
	c->Get2Prong()->UnsetOwnPrimaryVtx();
      }
    }

    
    // count 3prong candidates
    Int_t n3Prong = array3Prong->GetEntriesFast();
    nTot3Prong += n3Prong;
    cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
    
    // loop over HF vertices
    Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
    nTotHF += nVtxsHF;
    cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
    for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
      AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
      // print info
      //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
    }


    // count cascade candidates
    if (arrayCascades){
      Int_t nCasc = arrayCascades->GetEntriesFast();
      nTotCasc+=nCasc;
      cout << "Number of Cascades: "<<nCasc<<endl;
    }
    
  }
  
  printf("\n Total HF vertices: %d\n",nTotHF);
  printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
  printf("\n Total D*->D0pi: %d\n",nTotDstar);
  printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
  if (arrayCascades) printf("\n Total Cascades: %d\n",nTotCasc);

  TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
  c1->Divide(2,2);
  c1->cd(1);
  hCPtaVSd0d0->Draw("colz");
  c1->cd(2);
  hMass->SetFillColor(4);
  hMass->Draw();
  c1->cd(3);
  hSecVtxZ->SetFillColor(2);
  hSecVtxZ->Draw();
  c1->cd(4);
  hDeltaMassDstar->SetFillColor(3);
  hDeltaMassDstar->Draw();

  return;
}
//------------------------------------------------------------------------
void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root") 
{
  //
  // Example macro to read D0->Kpi candidates from a stand-alone
  // heavy-flavour AOD (i.e. without standard AOD) and apply cuts
  // Origin: A.Dainese
  //
  Bool_t useParFiles=kFALSE;
  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
  LoadLibraries(useParFiles);

  // create a test histogram
  TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
  hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
  hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
  TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
  hMass->SetXTitle("Invariant mass [GeV]");
  hMass->SetYTitle("Entries");
  TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
  hSecVtxZ->SetXTitle("z of decay vertex [cm]");
  hSecVtxZ->SetYTitle("Entries");

  // open input file and get the TTree
  TFile inFile(aodHFFileName,"READ");
  if (!inFile.IsOpen()) return;

  TTree *aodTree = (TTree*)inFile.Get("aodTree");

  AliAODEvent *aod = new AliAODEvent();

  aod->ReadFromTree(aodTree);

  // load heavy flavour vertices
  TClonesArray *arrayVerticesHF = 
    (TClonesArray*)aod->GetList()->FindObject("VerticesHF"); 
  
  // load D0->Kpi candidates
  TClonesArray *arrayD0toKpi = 
    (TClonesArray*)aod->GetList()->FindObject("D0toKpi"); 
     

  Double_t cutsD0[9]=
  // cutsD0[0] = inv. mass half width [GeV]
  // cutsD0[1] = dca [cm]
  // cutsD0[2] = cosThetaStar
  // cutsD0[3] = pTK [GeV/c]
  // cutsD0[4] = pTPi [GeV/c]
  // cutsD0[5] = d0K [cm]   upper limit!
  // cutsD0[6] = d0Pi [cm]  upper limit!
  // cutsD0[7] = d0d0 [cm^2]
  // cutsD0[8] = cosThetaPoint
                     {1000.,
		      100000.,
		      1.1,
		      0.,
		      0.,
		      100000.,
		      100000.,
		      100000000.,
		      -1.1}; 

  Int_t nTotHF=0,nTotD0toKpi=0;
  AliAODVertex *vtx1=0;

  // loop over events
  Int_t nEvents = aodTree->GetEntries();
  for (Int_t nEv = 0; nEv < nEvents; nEv++) {
    cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;

    // read event
    aodTree->GetEvent(nEv);
    //aodTree->BranchRef();

    
    // loop over D0->Kpi candidates
    Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
    nTotD0toKpi += nD0toKpi;
    cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
    
    for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
      AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
      Int_t okD0=0,okD0bar=0; 
      if(d->SelectD0(cutsD0,okD0,okD0bar)) {
	//cout<<1e8*d->Prodd0d0()<<endl;
	hMass->Fill(d->InvMassD0(),0.5);
	hMass->Fill(d->InvMassD0bar(),0.5);
	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
	hSecVtxZ->Fill(d->GetSecVtxZ());
	//cout<<d->GetSecVtxZ() <<endl;

      }
    }
    
  }
  
  printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);

  TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
  c->Divide(2,2);
  c->cd(1);
  hCPtaVSd0d0->Draw("colz");
  c->cd(2);
  hMass->SetFillColor(4);
  hMass->Draw();
  c->cd(3);
  hSecVtxZ->SetFillColor(2);
  hSecVtxZ->Draw();

  return;
}
Commit	Line	Data
95e5b6b5	1	void ReadAODVertexingHF(const char *aodFileName="AliAOD.root",
	2	const char *aodHFFileName="AliAOD.VertexingHF.root")
	3	{
	4	//
	5	// Example macro to read D0->Kpi candidates from AOD (having the
	6	// standard AOD + a friend heavy-flavour AOD) and apply cuts
	7	// Origin: A.Dainese
	8	//
	9
	10	Bool_t useParFiles=kFALSE;
60c2149d	11	gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
95e5b6b5	12	LoadLibraries(useParFiles);
	13
	14	// create a test histogram
	15	TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
	16	hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
	17	hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
	18	TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
	19	hMass->SetXTitle("Invariant mass [GeV]");
	20	hMass->SetYTitle("Entries");
	21	TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
	22	hSecVtxZ->SetXTitle("z of decay vertex [cm]");
	23	hSecVtxZ->SetYTitle("Entries");
	24	TH1F hDeltaMassDstar = new TH1F("hDeltaMassDstar","D delta mass plot",100,0,0.3);
	25	hDeltaMassDstar->SetXTitle("M(Kpipi)-M(Kpi) [GeV]");
	26	hDeltaMassDstar->SetYTitle("Entries");
	27
	28
	29	// open input file and get the TTree
	30	TFile inFile(aodFileName,"READ");
	31	if (!inFile.IsOpen()) return;
	32
	33	TTree aodTree = (TTree)inFile.Get("aodTree");
	34	aodTree->AddFriend("aodTree",aodHFFileName);
	35
	36	AliAODEvent *aod = new AliAODEvent();
	37
	38	aod->ReadFromTree(aodTree);
	39
	40	// load heavy flavour vertices
	41	TClonesArray *arrayVerticesHF =
	42	(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
	43
	44	// load D0->Kpi candidates
	45	TClonesArray *arrayD0toKpi =
	46	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
	47
	48	// load 3prong candidates
	49	TClonesArray *array3Prong =
	50	(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
	51
	52	// load D* candidates
	53	TClonesArray *arrayDstar =
	54	(TClonesArray*)aod->GetList()->FindObject("Dstar");
	55
	56	// load cascade (V0+track) candidates
	57	TClonesArray *arrayCascades =
	58	(TClonesArray*)aod->GetList()->FindObject("CascadesHF");
	59
	60
	61	Double_t cutsD0[9]=
	62	// cutsD0[0] = inv. mass half width [GeV]
	63	// cutsD0[1] = dca [cm]
	64	// cutsD0[2] = cosThetaStar
	65	// cutsD0[3] = pTK [GeV/c]
	66	// cutsD0[4] = pTPi [GeV/c]
	67	// cutsD0[5] = d0K [cm] upper limit!
	68	// cutsD0[6] = d0Pi [cm] upper limit!
	69	// cutsD0[7] = d0d0 [cm^2]
	70	// cutsD0[8] = cosThetaPoint
	71	{1000.,
	72	100000.,
	73	1.1,
	74	0.,
	75	0.,
76	100000.,
77	100000.,
78	100000000.,
79	-1.1};
80	Double_t cutsDstar[5]=
81	// (to be passed to AliAODRecoCascadeHF::SelectDstar())
82	// 0 = inv. mass half width of D* [GeV]
83	// 1 = half width of (M_Kpipi-M_Kpi) [GeV]
84	// 2 = PtMin of pi_s [GeV/c]
85	// 3 = PtMax of pi_s [GeV/c]
86	// 4 = theta, angle between the trace of pi_s and D0 decay plane [rad]
87	{999999.,
88	999999.,
89	0.1,
90	1.0,
91	0.5};
92
93	Double_t cutsLctoV0[9]=// cuts on Lambdac candidates to V0+bachelor
94	// (to be passed to AliAODRecoDecayHF3Prong::SelectLctoV0())
95	// 0 = inv. mass half width in K0s hypothesis [GeV]
96	// 1 = inv. mass half width in Lambda hypothesis [GeV]
97	// 2 = inv. mass V0 in K0s hypothesis half width [GeV]
98	// 3 = inv. mass V0 in Lambda hypothesis half width [GeV]
99	// 4 = pT min Bachelor track [GeV/c]
100	// 5 = pT min V0-Positive track [GeV/c]
101	// 6 = pT min V0-Negative track [GeV/c]
102	// 7 = dca cut on the V0 (cm)
103	// 8 = dca cut on the cascade (cm)
104	{2.0,2.0,1.0,1.0,0.0,0.0,0.0,1000.,1000.};
105
106	Int_t nTotHF=0,nTotD0toKpi=0,nTotDstar=0,nTot3Prong=0,nTotCasc=0;
107	AliAODVertex *vtx1=0;
108
109	// loop over events
110	Int_t nEvents = aodTree->GetEntries();
111	for (Int_t nEv = 0; nEv < nEvents; nEv++) {
112	cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
113
114	// read event
115	aodTree->GetEvent(nEv);
116	//aodTree->BranchRef();
117
118	//print event info
119	aod->GetHeader()->Print();
120
121	// primary vertex
122	vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
123	vtx1->Print();
124
125	// make trkIDtoEntry register (temporary)
126	Int_t trkIDtoEntry[100000];
127	for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
128	AliAODTrack *track = aod->GetTrack(it);
129	trkIDtoEntry[track->GetID()]=it;
130	}
131
132
133	// loop over D0->Kpi candidates
134	Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
135	nTotD0toKpi += nD0toKpi;
136	cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
137
138	for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
139	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)arrayD0toKpi->UncheckedAt(iD0toKpi);
140
141	Bool_t unsetvtx=kFALSE;
142	if(!d->GetOwnPrimaryVtx()) {
143	d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
144	unsetvtx=kTRUE;
145	}
146	Int_t okD0=0,okD0bar=0;
147	if(d->SelectD0(cutsD0,okD0,okD0bar)) {
148	//cout<<1e8*d->Prodd0d0()<<endl;
149	hMass->Fill(d->InvMassD0(),0.5);
150	hMass->Fill(d->InvMassD0bar(),0.5);
151	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
152	hSecVtxZ->Fill(d->GetSecVtxZ());
153	//cout<<d->GetSecVtxX() <<endl;
154
155	// get daughter AOD tracks
156	AliAODTrack trk0 = (AliAODTrack)d->GetDaughter(0);
157	AliAODTrack trk1 = (AliAODTrack)d->GetDaughter(1);
158	if(!trk0 \|\| !trk1) {
159	trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
160	trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
161	cout<<"references to standard AOD not available"<<endl;
162	}
163	cout<<"pt of positive track: "<<trk0->Pt()<<endl;
164
165	// make a AliNeutralTrackParam from the D0
166	// and calculate impact parameters to primary vertex
167	AliNeutralTrackParam trackD0(d);
168	cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
169	//trackD0.Print();
170	Double_t dz[2],covdz[3];
171	trackD0.PropagateToDCA(vtx1,aod->GetMagneticField(),1000.,dz,covdz);
172	cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
173	}
174	if(unsetvtx) d->UnsetOwnPrimaryVtx();
175	}
176
177
178	// loop over D* candidates
179	Int_t nDstar = arrayDstar->GetEntriesFast();
180	nTotDstar += nDstar;
181	cout<<"Number of D*->D0pi: "<<nDstar<<endl;
182
183	for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) {
184	AliAODRecoCascadeHF c = (AliAODRecoCascadeHF)arrayDstar->UncheckedAt(iDstar);
185	Bool_t unsetvtx=kFALSE;
186	if(!c->GetOwnPrimaryVtx()) {
187	c->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
188	c->Get2Prong()->SetOwnPrimaryVtx(vtx1);
189	unsetvtx=kTRUE;
190	}
191	if(c->SelectDstar(cutsDstar,cutsD0)) {
192	hDeltaMassDstar->Fill(c->DeltaInvMass());
193	// get daughters
194	AliAODTrack trk = (AliAODTrack)c->GetBachelor();
195	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)c->Get2Prong();
196	cout<<"pt of soft pi: "<<trk->Pt()<<endl;
197	cout<<"pt of D0: "<<d->Pt()<<endl;
198	}
199
200	if(unsetvtx) {
201	c->UnsetOwnPrimaryVtx();
202	c->Get2Prong()->UnsetOwnPrimaryVtx();
203	}
204	}
205
206
207	// count 3prong candidates
208	Int_t n3Prong = array3Prong->GetEntriesFast();
209	nTot3Prong += n3Prong;
210	cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
211
212	// loop over HF vertices
213	Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
214	nTotHF += nVtxsHF;
215	cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
216	for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
217	AliAODVertex vtxHF = (AliAODVertex)arrayVerticesHF->UncheckedAt(iVtx);
218	// print info
219	//cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
220	}
221
222
223	// count cascade candidates
224	if (arrayCascades){
225	Int_t nCasc = arrayCascades->GetEntriesFast();
226	nTotCasc+=nCasc;
227	cout << "Number of Cascades: "<<nCasc<<endl;
228	}
229
230	}
231
232	printf("\n Total HF vertices: %d\n",nTotHF);
233	printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
234	printf("\n Total D*->D0pi: %d\n",nTotDstar);
235	printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
236	if (arrayCascades) printf("\n Total Cascades: %d\n",nTotCasc);
237
238	TCanvas *c1 = new TCanvas("c1","c1",0,0,800,700);
239	c1->Divide(2,2);
240	c1->cd(1);
241	hCPtaVSd0d0->Draw("colz");
242	c1->cd(2);
243	hMass->SetFillColor(4);
244	hMass->Draw();
245	c1->cd(3);
246	hSecVtxZ->SetFillColor(2);
247	hSecVtxZ->Draw();
248	c1->cd(4);
249	hDeltaMassDstar->SetFillColor(3);
250	hDeltaMassDstar->Draw();
251
252	return;
253	}
254	//------------------------------------------------------------------------
255	void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
256	{
257	//
258	// Example macro to read D0->Kpi candidates from a stand-alone
259	// heavy-flavour AOD (i.e. without standard AOD) and apply cuts
260	// Origin: A.Dainese
261	//
262	Bool_t useParFiles=kFALSE;
263	gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/LoadLibraries.C");
264	LoadLibraries(useParFiles);
265
266	// create a test histogram
267	TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
268	hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
269	hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
270	TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
271	hMass->SetXTitle("Invariant mass [GeV]");
272	hMass->SetYTitle("Entries");
273	TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
274	hSecVtxZ->SetXTitle("z of decay vertex [cm]");
275	hSecVtxZ->SetYTitle("Entries");
276
277	// open input file and get the TTree
278	TFile inFile(aodHFFileName,"READ");
279	if (!inFile.IsOpen()) return;
280
281	TTree aodTree = (TTree)inFile.Get("aodTree");
282
283	AliAODEvent *aod = new AliAODEvent();
284
285	aod->ReadFromTree(aodTree);
286
287	// load heavy flavour vertices
288	TClonesArray *arrayVerticesHF =
289	(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
290
291	// load D0->Kpi candidates
292	TClonesArray *arrayD0toKpi =
293	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
294
295
296	Double_t cutsD0[9]=
297	// cutsD0[0] = inv. mass half width [GeV]
298	// cutsD0[1] = dca [cm]
299	// cutsD0[2] = cosThetaStar
300	// cutsD0[3] = pTK [GeV/c]
301	// cutsD0[4] = pTPi [GeV/c]
302	// cutsD0[5] = d0K [cm] upper limit!
303	// cutsD0[6] = d0Pi [cm] upper limit!
304	// cutsD0[7] = d0d0 [cm^2]
305	// cutsD0[8] = cosThetaPoint
306	{1000.,
307	100000.,
308	1.1,
309	0.,
310	0.,
311	100000.,
312	100000.,
313	100000000.,
314	-1.1};
315
316	Int_t nTotHF=0,nTotD0toKpi=0;
317	AliAODVertex *vtx1=0;
318
319	// loop over events
320	Int_t nEvents = aodTree->GetEntries();
321	for (Int_t nEv = 0; nEv < nEvents; nEv++) {
322	cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
323
324	// read event
325	aodTree->GetEvent(nEv);
326	//aodTree->BranchRef();
327
328
329	// loop over D0->Kpi candidates
330	Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
331	nTotD0toKpi += nD0toKpi;
332	cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
333
334	for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
335	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)arrayD0toKpi->UncheckedAt(iD0toKpi);
336	Int_t okD0=0,okD0bar=0;
337	if(d->SelectD0(cutsD0,okD0,okD0bar)) {
338	//cout<<1e8*d->Prodd0d0()<<endl;
339	hMass->Fill(d->InvMassD0(),0.5);
340	hMass->Fill(d->InvMassD0bar(),0.5);
341	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
342	hSecVtxZ->Fill(d->GetSecVtxZ());
343	//cout<<d->GetSecVtxZ() <<endl;
344
345	}
346	}
347
348	}
349
350	printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
351
352	TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
353	c->Divide(2,2);
354	c->cd(1);
355	hCPtaVSd0d0->Draw("colz");
356	c->cd(2);
357	hMass->SetFillColor(4);
358	hMass->Draw();
359	c->cd(3);
360	hSecVtxZ->SetFillColor(2);
361	hSecVtxZ->Draw();
362
363	return;
364	}