[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / 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
  //
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libAOD.so");
  gSystem->Load("libCORRFW.so");
  gSystem->Load("libPWG3base.so");
  gSystem->Load("libPWG3vertexingHF.so");

  // 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"); 
    

  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};


  Int_t nTotHF=0,nTotD0toKpi=0,nTotDstar=0,nTot3Prong=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;
    }
    
  }
  
  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);

  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
  //
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libAOD.so");
  gSystem->Load("libPWG3base.so");

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