[u/mrichter/AliRoot.git] / PWGHF / 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
  //

  TStopwatch t;
  t.Start();

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

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

  AliRDHFCutsD0toKpi *cutsD0toKpi = new AliRDHFCutsD0toKpi("CutsD0toKpi");
  cutsD0toKpi->SetStandardCutsPP2010();
  cutsD0toKpi->SetRemoveDaughtersFromPrim(kFALSE);

  AliRDHFCutsDStartoKpipi *cutsDStar = new AliRDHFCutsDStartoKpipi("CutsDStartoKpipi");
  cutsDStar->SetStandardCutsPP2010();

  AliRDHFCutsDplustoKpipi *cutsDplustoKpipi = new AliRDHFCutsDplustoKpipi("CutsDplustoKpipi");
  cutsDplustoKpipi->SetStandardCutsPP2010();
  cutsDplustoKpipi->SetRemoveDaughtersFromPrim(kFALSE);
  
  Int_t nTot3ProngSele=0;
  Int_t nTotD0toKpiSele=0;
  Int_t nTotDStarSele=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;
    }
    */

    // Fix references to daughter tracks
    //AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF();
    //fixer->FixReferences(aod);
    //delete fixer;
    //
    //AliRDHFCutsD0toKpi *mycuts=new AliRDHFCutsD0toKpi();
    //mycuts->SetFixRefs(kTRUE);
    //mycuts->IsEventSelected(aod);

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

      d->SetOwnPrimaryVtx(vtx1);

      /*
	// get daughter AOD tracks
	AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
	AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);

	if(trk0->GetStatus()) printf("ok %d\n",iD0toKpi);
      */
      printf("D0 %d\n",iD0toKpi);
      if(cutsD0toKpi->IsSelected(d,AliRDHFCuts::kAll)) {printf("D0 %d passed\n",iD0toKpi); nTotD0toKpiSele++;}

    }


    // 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);
      printf("D* %d\n",iDstar);
      if(cutsDStar->IsSelected(c,AliRDHFCuts::kCandidate)) {printf("D* %d passed\n",iDstar); nTotDStarSele++;}

    }

  
    // count 3prong candidates
    Int_t n3Prong = array3Prong->GetEntriesFast();
    nTot3Prong += n3Prong;
    cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;

    for (Int_t i3p = 0; i3p < n3Prong; i3p++) {
      AliAODRecoDecayHF3Prong *ccc = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3p);
      printf("3p %d\n",i3p);
      //if(cutsDplustoKpipi->IsSelected(ccc,AliRDHFCuts::kCandidate)) {printf("3p %d passed\n",i3p); nTot3ProngSele++;}

    }

    /*
    // 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 <0; 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; selected %d\n",nTotD0toKpi,nTotD0toKpiSele);
  printf("\n Total D*->D0pi: %d; selected %d\n",nTotDstar,nTotDStarSele);
  printf("\n Total Charm->3Prong: %d; selected %d\n",nTot3Prong,nTot3ProngSele);
  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();
  */

  t.Stop();
  t.Print();

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