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

  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,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->GetSecVtxZ() <<endl;

	// get daughter AOD tracks
	AliAODTrack *trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
	AliAODTrack *trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
	cout<<"pt of positive track: "<<trk0->Pt()<<endl;
	// this will be replaced by 
	//AliAODTrack *trk0 = (AliAODTrack*)(d->GetSecondaryVtx()->GetDaughter(0));
      }
      if(unsetvtx) d->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 Charm->3Prong: %d\n",nTot3Prong);

  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;
}
//------------------------------------------------------------------------
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");
	12	gSystem->Load("libPWG3base.so");
	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
	25	// open input file and get the TTree
	26	TFile inFile(aodFileName,"READ");
	27	if (!inFile.IsOpen()) return;
	28
	29	TTree aodTree = (TTree)inFile.Get("aodTree");
	30	aodTree->AddFriend("aodTree",aodHFFileName);
	31
	32	AliAODEvent *aod = new AliAODEvent();
	33
	34	aod->ReadFromTree(aodTree);
	35
	36	// load heavy flavour vertices
	37	TClonesArray *arrayVerticesHF =
	38	(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
	39
	40	// load D0->Kpi candidates
	41	TClonesArray *arrayD0toKpi =
	42	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
	43
	44	// load 3prong candidates
	45	TClonesArray *array3Prong =
	46	(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
	47
	48
	49	Double_t cutsD0[9]=
	50	// cutsD0[0] = inv. mass half width [GeV]
	51	// cutsD0[1] = dca [cm]
	52	// cutsD0[2] = cosThetaStar
	53	// cutsD0[3] = pTK [GeV/c]
	54	// cutsD0[4] = pTPi [GeV/c]
	55	// cutsD0[5] = d0K [cm] upper limit!
	56	// cutsD0[6] = d0Pi [cm] upper limit!
	57	// cutsD0[7] = d0d0 [cm^2]
	58	// cutsD0[8] = cosThetaPoint
	59	{1000.,
	60	100000.,
	61	1.1,
	62	0.,
	63	0.,
	64	100000.,
	65	100000.,
	66	100000000.,
	67	-1.1};
	68
	69	Int_t nTotHF=0,nTotD0toKpi=0,nTot3Prong=0;
	70	AliAODVertex *vtx1=0;
71
72	// loop over events
73	Int_t nEvents = aodTree->GetEntries();
74	for (Int_t nEv = 0; nEv < nEvents; nEv++) {
75	cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
76
77	// read event
78	aodTree->GetEvent(nEv);
79	//aodTree->BranchRef();
80
81	//print event info
82	aod->GetHeader()->Print();
83
84	// primary vertex
85	vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
86	vtx1->Print();
87
88	// make trkIDtoEntry register (temporary)
89	Int_t trkIDtoEntry[100000];
90	for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
91	AliAODTrack *track = aod->GetTrack(it);
92	trkIDtoEntry[track->GetID()]=it;
93	}
94
95	// loop over D0->Kpi candidates
96	Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
97	nTotD0toKpi += nD0toKpi;
98	cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
99
100	for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
101	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)arrayD0toKpi->UncheckedAt(iD0toKpi);
102	Bool_t unsetvtx=kFALSE;
103	if(!d->GetOwnPrimaryVtx()) {
104	d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
105	unsetvtx=kTRUE;
106	}
107	Int_t okD0=0,okD0bar=0;
108	if(d->SelectD0(cutsD0,okD0,okD0bar)) {
109	//cout<<1e8*d->Prodd0d0()<<endl;
110	hMass->Fill(d->InvMassD0(),0.5);
111	hMass->Fill(d->InvMassD0bar(),0.5);
112	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
113	hSecVtxZ->Fill(d->GetSecVtxZ());
114	//cout<<d->GetSecVtxZ() <<endl;
115
116	// get daughter AOD tracks
117	AliAODTrack *trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
118	AliAODTrack *trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
119	cout<<"pt of positive track: "<<trk0->Pt()<<endl;
120	// this will be replaced by
121	//AliAODTrack trk0 = (AliAODTrack)(d->GetSecondaryVtx()->GetDaughter(0));
122	}
123	if(unsetvtx) d->UnsetOwnPrimaryVtx();
124	}
fcca9870	125
	126	// count 3prong candidates
	127	Int_t n3Prong = array3Prong->GetEntriesFast();
	128	nTot3Prong += n3Prong;
	129	cout<<"Number of Charm->3Prong: "<<n3Prong<<endl;
	130
	131	// loop over HF vertices
	132	Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
	133	nTotHF += nVtxsHF;
	134	cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
	135	for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
	136	AliAODVertex vtxHF = (AliAODVertex)arrayVerticesHF->UncheckedAt(iVtx);
	137	// print info
	138	//cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
	139	}
	140
	141	}
	142
	143	printf("\n Total HF vertices: %d\n",nTotHF);
	144	printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
	145	printf("\n Total Charm->3Prong: %d\n",nTot3Prong);
	146
	147	TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
	148	c->Divide(2,2);
	149	c->cd(1);
	150	hCPtaVSd0d0->Draw("colz");
	151	c->cd(2);
	152	hMass->SetFillColor(4);
	153	hMass->Draw();
	154	c->cd(3);
	155	hSecVtxZ->SetFillColor(2);
	156	hSecVtxZ->Draw();
	157
	158	return;
	159	}
0bd2832b	160	//------------------------------------------------------------------------
	161	void ReadAODVertexingHFsa(const char *aodHFFileName="AliAOD.VertexingHF.sa.root")
	162	{
	163	//
	164	// Example macro to read D0->Kpi candidates from a stand-alone
	165	// heavy-flavour AOD (i.e. without standard AOD) and apply cuts
	166	// Origin: A.Dainese
	167	//
	168	gSystem->Load("libANALYSIS.so");
	169	gSystem->Load("libANALYSISalice.so");
	170	gSystem->Load("libAOD.so");
	171	gSystem->Load("libPWG3base.so");
	172
	173	// create a test histogram
	174	TH2F *hCPtaVSd0d0 = new TH2F("hCPtaVSd0d0","D^{0} correlation plot",1000,-50000,50000,1000,-1,1);
	175	hCPtaVSd0d0->SetXTitle("Product of impact parameters [#mu m^{2}]");
	176	hCPtaVSd0d0->SetYTitle("Cosine of pointing angle");
	177	TH1F *hMass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
	178	hMass->SetXTitle("Invariant mass [GeV]");
	179	hMass->SetYTitle("Entries");
	180	TH1F *hSecVtxZ = new TH1F("hSecVtxZ","D^{0} decay vertex z",1000,-10,10);
	181	hSecVtxZ->SetXTitle("z of decay vertex [cm]");
	182	hSecVtxZ->SetYTitle("Entries");
	183
	184	// open input file and get the TTree
	185	TFile inFile(aodHFFileName,"READ");
	186	if (!inFile.IsOpen()) return;
	187
	188	TTree aodTree = (TTree)inFile.Get("aodTree");
	189
	190	AliAODEvent *aod = new AliAODEvent();
	191
	192	aod->ReadFromTree(aodTree);
	193
	194	// load heavy flavour vertices
	195	TClonesArray *arrayVerticesHF =
	196	(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
	197
	198	// load D0->Kpi candidates
	199	TClonesArray *arrayD0toKpi =
	200	(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
	201
	202
	203	Double_t cutsD0[9]=
	204	// cutsD0[0] = inv. mass half width [GeV]
	205	// cutsD0[1] = dca [cm]
	206	// cutsD0[2] = cosThetaStar
	207	// cutsD0[3] = pTK [GeV/c]
	208	// cutsD0[4] = pTPi [GeV/c]
	209	// cutsD0[5] = d0K [cm] upper limit!
	210	// cutsD0[6] = d0Pi [cm] upper limit!
	211	// cutsD0[7] = d0d0 [cm^2]
	212	// cutsD0[8] = cosThetaPoint
	213	{1000.,
	214	100000.,
	215	1.1,
	216	0.,
	217	0.,
	218	100000.,
	219	100000.,
	220	100000000.,
	221	-1.1};
	222
	223	Int_t nTotHF=0,nTotD0toKpi=0;
224	AliAODVertex *vtx1=0;
225
226	// loop over events
227	Int_t nEvents = aodTree->GetEntries();
228	for (Int_t nEv = 0; nEv < nEvents; nEv++) {
229	cout<<"\n------------ Event: "<<nEv<<" ------------------"<<endl;
230
231	// read event
232	aodTree->GetEvent(nEv);
233	//aodTree->BranchRef();
234
235
236	// loop over D0->Kpi candidates
237	Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
238	nTotD0toKpi += nD0toKpi;
239	cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
240
241	for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
242	AliAODRecoDecayHF2Prong d = (AliAODRecoDecayHF2Prong)arrayD0toKpi->UncheckedAt(iD0toKpi);
243	Int_t okD0=0,okD0bar=0;
244	if(d->SelectD0(cutsD0,okD0,okD0bar)) {
245	//cout<<1e8*d->Prodd0d0()<<endl;
246	hMass->Fill(d->InvMassD0(),0.5);
247	hMass->Fill(d->InvMassD0bar(),0.5);
248	hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
249	hSecVtxZ->Fill(d->GetSecVtxZ());
250	//cout<<d->GetSecVtxZ() <<endl;
251
252	}
253	}
254
255	}
256
257	printf("\n Total D0->Kpi: %d\n",nTotD0toKpi);
258
259	TCanvas *c = new TCanvas("c","c",0,0,1000,1000);
260	c->Divide(2,2);
261	c->cd(1);
262	hCPtaVSd0d0->Draw("colz");
263	c->cd(2);
264	hMass->SetFillColor(4);
265	hMass->Draw();
266	c->cd(3);
267	hSecVtxZ->SetFillColor(2);
268	hSecVtxZ->Draw();
269
270	return;
271	}