[u/mrichter/AliRoot.git] / PWG / muon / ReadAOD_MCBranch.C

void ReadAOD_MCBranch(char *filein="AliAODs.root"){

//=============================================================================
// Macro to read AOD+MC Branch: 
// in this example the macro reads generated and reconstructed JPsi 
// Plots refer to all MC generated particles, to MC particles which have
// been reconstructed and reconstructed tracks
//
// Each aod event has 3 MC particles (J/Psi, mu1, mu2)
// To access MC muons:   mctrackall->IsPhysicalPrimary()
// To access JPsi muons: mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()

// AliStack::IsPhysicalPrimary() 
// Test if a particle is a physical primary according to the following definition:
// Particles produced in the collision including products of strong and
// electromagnetic decay and excluding feed-down from weak decays of strange
// particles.
//=============================================================================

  gStyle->SetFrameFillColor(0);
  gStyle->SetCanvasColor(0);
  gStyle->SetPalette(1); 

// Open AOD file    
  TFile *file = new TFile(filein);
  TTree *aodTree = file->Get("aodTree");
  
// Read AOD event    
  AliAODEvent *event = new AliAODEvent();
  event->ReadFromTree(aodTree);
  
// Read MC info   
  TClonesArray *mcarray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
  if(!mcarray) {
    printf("No MC info in AOD!\n");
    return;
  }  
  
// Book Histos
  TH2D *hpt_mc_rec = new TH2D("hpt_mc_rec","Pt correlations MC-REC",50,0.,10.,50,0.,10.);
  TH1D *hpt_rec = new TH1D("hpt_rec","Pt REC tracks",20,0.,10.);
  TH1D *hpt_mc = new TH1D("hpt_mc","PT MC particles corresponding to rec tracks",20,0.,10.);
  TH1D *hpt_allmc = new TH1D("hpt_allmc","PT all MC particles",20,0.,10.);
  TH1D *hpt_primarymc = new TH1D("hpt_primarymc","PT MC primaries (J/Psi and muons)",20,0.,10.);
  TH1D *hpt_physprimarymc = new TH1D("hpt_physprimarymc","PT Muons MC",20,0.,10.);
  TH1D *hpt_notphysprimarymc = new TH1D("hpt_notphysprimarymc","Pt JPsi MC",20,0.,10.);
  TH1D *hpt_diff = new TH1D("hpt_diff","Pt: Rec - MC",100,-5.,5.);

  TH1D *hCharge_primarymc = new TH1D("hCharge_primarymc","Charge MC primary particle",10,-5.,5.);
  TH1D *hCharge_physprimmc = new TH1D("hCharge_physprimmc","Charge Muons MC",10,-5.,5.);
  TH1D *hCharge_notphysprimmc = new TH1D("hCharge_notphysprimmc","Charge JPsi MC",10,-5.,5.);
  TH1D *hCharge_rec = new TH1D("hCharge_rec","hCharge_rec",10,-5.,5.);
  TH1D *hCharge_mc = new TH1D("hCharge_mc","hCharge_mc",10,-5.,5.);

  TH1D *hMass_mc = new TH1D("hMass_mc","Mass MC particles",100,0.,5..);
  TH1D *hMass_rec = new TH1D("hMass_rec","Mass rec tracks",100,0.,5..);
  TH1D *hMass_allmc = new TH1D("hMass_allmc","Mass all MC particles",100,0.,5..);
  TH1D *hMass_primarymc = new TH1D("hMass_primarymc","Mass primary particles MC",100,0.,5..);
  TH1D *hMass_physprimarymc = new TH1D("hMass_physprimarymc","Mass muons MC",100,0.,5..);
  TH1D *hMass_notphysprimarymc = new TH1D("hMass_notphysprimarymc","Mass JPsi MC",100,0.,5..);
  TH1D *hMassJPsi_mc = new TH1D("hMassJPsi_mc","Mass JPsi MC",100,0.,5..);
  
// Loop on AOD events
  Int_t nev=0;

  while(aodTree->GetEvent(nev++)){  

// Loop on all MC tracks    
    for(int ii=0;ii<mcarray->GetEntries();ii++){
      AliAODMCParticle *mctrackall = (AliAODMCParticle*) mcarray->At(ii);     
      // all particles
      hpt_allmc->Fill(mctrackall->Pt());
      hMass_allmc->Fill(mctrackall->M());
      
      // Muons
      if(mctrackall->IsPhysicalPrimary()) {
        hCharge_physprimmc->Fill(mctrackall->Charge());
        hpt_physprimarymc->Fill(mctrackall->Pt());
        hMass_physprimarymc->Fill(mctrackall->M());
      }
      	
      if(mctrackall->IsPrimary()) {
        hpt_primarymc->Fill(mctrackall->Pt());
        hMass_primarymc->Fill(mctrackall->M());
	hCharge_primarymc->Fill(mctrackall->Charge());
      }
        
      // JPsi
      if(mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()) {
        hCharge_notphysprimmc->Fill(mctrackall->Charge());
        hMass_notphysprimarymc->Fill(mctrackall->M());
        hpt_notphysprimarymc->Fill(mctrackall->Pt());
      }	
    }
    
// Loop on reconstructed tracks    
    for(int i=0;i<event->GetNumberOfTracks();i++){
      AliAODTrack *track = dynamic_cast<AliAODTrack*>(event->GetTrack(i));
      if(!track) {
        printf("Error: no track! \n");
        continue;
      } 
      hpt_rec->Fill(track->Pt());
      hMass_rec->Fill(track->M());
      hCharge_rec->Fill(track->Charge());
      
// Read MC particle corresponding to reconstructed tracks   
      AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(track->GetLabel());     
      if(!mctrack){
        printf("Error: no MC track! \n");
	continue;
      }
      //printf("nev=%d i=%d mc=%f rec=%f\n",nev-1,i,mctrack->Pt(),track->Pt());
      hpt_mc->Fill(mctrack->Pt());
      hpt_mc_rec->Fill(track->Pt(),mctrack->Pt());
      hpt_diff->Fill(track->Pt()-mctrack->Pt());
      hMass_mc->Fill(mctrack->M());
      hCharge_mc->Fill(mctrack->Charge());
   }
 }
 
// Plots 
 TCanvas *c = new TCanvas("c","PT spectra",20,20,600,600);
 c->Divide(3,3);
 c->cd(1);
 hpt_primarymc->Draw();
 hpt_primarymc->GetXaxis()->SetTitle("Pt MC primary particles");
 c->cd(2);  
 hpt_mc->GetXaxis()->SetTitle("Pt MC particles");
 hpt_mc->Draw();
 hpt_mc->SetLineColor(4);
 c->cd(3);
 hpt_rec->GetXaxis()->SetTitle("Rec Pt track");
 hpt_rec->Draw();
 hpt_rec->SetLineColor(2);
 c->cd(4);
 hpt_physprimarymc->Draw();
 hpt_physprimarymc->GetXaxis()->SetTitle("Pt MC muon track");
 c->cd(5);
 hpt_notphysprimarymc->Draw();
 hpt_notphysprimarymc->GetXaxis()->SetTitle("Pt MC JPsi track");
 c->cd(7);
 hpt_mc_rec->Draw("colz");
 hpt_mc_rec->GetXaxis()->SetTitle("Pt rec track");
 hpt_mc_rec->GetYaxis()->SetTitle("Pt MC track");
 c->cd(8);
 hpt_diff->Draw();
 hpt_diff->GetXaxis()->SetTitle("Rec-MC Pt");

 TCanvas *c1 = new TCanvas("c1","Mass distributions",40,40,600,600);
 c1->Divide(2,2);
 c1->cd(1);
 hMass_primarymc->Draw();
 hMass_primarymc->GetXaxis()->SetTitle("Mass MC primary particles");
 c1->cd(2);
 hMass_physprimarymc->Draw();
 hMass_physprimarymc->GetXaxis()->SetTitle("Mass MC muon particles");
 c1->cd(3);
 hMass_notphysprimarymc->Draw();
 hMass_notphysprimarymc->GetXaxis()->SetTitle("Mass MC J/Psi");
 c1->cd(4);
 hMass_rec->Draw();
 hMass_rec->GetXaxis()->SetTitle("Mass Rec tracks");

 TCanvas *c2= new TCanvas("c2","Charge distributions",60,60,600,600);
 c2->Divide(2,2);
 c2->cd(1);
 hCharge_primarymc->Draw();
 hCharge_primarymc->GetXaxis()->SetTitle("Charge MC primary particles");
 c2->cd(2);
 hCharge_physprimmc->Draw();
 hCharge_physprimmc->GetXaxis()->SetTitle("Charge muon MC");;
 c2->cd(3);
 hCharge_notphysprimmc->Draw();
 hCharge_notphysprimmc->GetXaxis()->SetTitle("Charge J/Psi MC");;
 c2->cd(4);
 hCharge_rec->Draw();
 hCharge_rec->GetXaxis()->SetTitle("Charge rec tracks");
}
Commit	Line	Data
af2b48fb	1	void ReadAOD_MCBranch(char *filein="AliAODs.root"){
	2
	3	//=============================================================================
	4	// Macro to read AOD+MC Branch:
	5	// in this example the macro reads generated and reconstructed JPsi
	6	// Plots refer to all MC generated particles, to MC particles which have
	7	// been reconstructed and reconstructed tracks
	8	//
	9	// Each aod event has 3 MC particles (J/Psi, mu1, mu2)
	10	// To access MC muons: mctrackall->IsPhysicalPrimary()
	11	// To access JPsi muons: mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()
	12
	13	// AliStack::IsPhysicalPrimary()
	14	// Test if a particle is a physical primary according to the following definition:
	15	// Particles produced in the collision including products of strong and
	16	// electromagnetic decay and excluding feed-down from weak decays of strange
	17	// particles.
	18	//=============================================================================
	19
	20	gStyle->SetFrameFillColor(0);
	21	gStyle->SetCanvasColor(0);
	22	gStyle->SetPalette(1);
	23
	24	// Open AOD file
	25	TFile *file = new TFile(filein);
	26	TTree *aodTree = file->Get("aodTree");
	27
	28	// Read AOD event
	29	AliAODEvent *event = new AliAODEvent();
	30	event->ReadFromTree(aodTree);
	31
	32	// Read MC info
	33	TClonesArray mcarray = dynamic_cast<TClonesArray>(event->FindListObject(AliAODMCParticle::StdBranchName()));
	34	if(!mcarray) {
	35	printf("No MC info in AOD!\n");
	36	return;
	37	}
	38
	39	// Book Histos
	40	TH2D *hpt_mc_rec = new TH2D("hpt_mc_rec","Pt correlations MC-REC",50,0.,10.,50,0.,10.);
	41	TH1D *hpt_rec = new TH1D("hpt_rec","Pt REC tracks",20,0.,10.);
	42	TH1D *hpt_mc = new TH1D("hpt_mc","PT MC particles corresponding to rec tracks",20,0.,10.);
	43	TH1D *hpt_allmc = new TH1D("hpt_allmc","PT all MC particles",20,0.,10.);
	44	TH1D *hpt_primarymc = new TH1D("hpt_primarymc","PT MC primaries (J/Psi and muons)",20,0.,10.);
	45	TH1D *hpt_physprimarymc = new TH1D("hpt_physprimarymc","PT Muons MC",20,0.,10.);
	46	TH1D *hpt_notphysprimarymc = new TH1D("hpt_notphysprimarymc","Pt JPsi MC",20,0.,10.);
	47	TH1D *hpt_diff = new TH1D("hpt_diff","Pt: Rec - MC",100,-5.,5.);
	48
	49	TH1D *hCharge_primarymc = new TH1D("hCharge_primarymc","Charge MC primary particle",10,-5.,5.);
	50	TH1D *hCharge_physprimmc = new TH1D("hCharge_physprimmc","Charge Muons MC",10,-5.,5.);
	51	TH1D *hCharge_notphysprimmc = new TH1D("hCharge_notphysprimmc","Charge JPsi MC",10,-5.,5.);
	52	TH1D *hCharge_rec = new TH1D("hCharge_rec","hCharge_rec",10,-5.,5.);
	53	TH1D *hCharge_mc = new TH1D("hCharge_mc","hCharge_mc",10,-5.,5.);
	54
	55	TH1D *hMass_mc = new TH1D("hMass_mc","Mass MC particles",100,0.,5..);
	56	TH1D *hMass_rec = new TH1D("hMass_rec","Mass rec tracks",100,0.,5..);
	57	TH1D *hMass_allmc = new TH1D("hMass_allmc","Mass all MC particles",100,0.,5..);
	58	TH1D *hMass_primarymc = new TH1D("hMass_primarymc","Mass primary particles MC",100,0.,5..);
	59	TH1D *hMass_physprimarymc = new TH1D("hMass_physprimarymc","Mass muons MC",100,0.,5..);
	60	TH1D *hMass_notphysprimarymc = new TH1D("hMass_notphysprimarymc","Mass JPsi MC",100,0.,5..);
	61	TH1D *hMassJPsi_mc = new TH1D("hMassJPsi_mc","Mass JPsi MC",100,0.,5..);
	62
	63	// Loop on AOD events
	64	Int_t nev=0;
65
66	while(aodTree->GetEvent(nev++)){
67
68	// Loop on all MC tracks
69	for(int ii=0;ii<mcarray->GetEntries();ii++){
70	AliAODMCParticle mctrackall = (AliAODMCParticle) mcarray->At(ii);
71	// all particles
72	hpt_allmc->Fill(mctrackall->Pt());
73	hMass_allmc->Fill(mctrackall->M());
74
75	// Muons
76	if(mctrackall->IsPhysicalPrimary()) {
77	hCharge_physprimmc->Fill(mctrackall->Charge());
78	hpt_physprimarymc->Fill(mctrackall->Pt());
79	hMass_physprimarymc->Fill(mctrackall->M());
80	}
81
82	if(mctrackall->IsPrimary()) {
83	hpt_primarymc->Fill(mctrackall->Pt());
84	hMass_primarymc->Fill(mctrackall->M());
85	hCharge_primarymc->Fill(mctrackall->Charge());
86	}
87
88	// JPsi
89	if(mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()) {
90	hCharge_notphysprimmc->Fill(mctrackall->Charge());
91	hMass_notphysprimarymc->Fill(mctrackall->M());
92	hpt_notphysprimarymc->Fill(mctrackall->Pt());
93	}
94	}
95
96	// Loop on reconstructed tracks
97	for(int i=0;i<event->GetNumberOfTracks();i++){
98	AliAODTrack track = dynamic_cast<AliAODTrack>(event->GetTrack(i));
99	if(!track) {
100	printf("Error: no track! \n");
101	continue;
102	}
103	hpt_rec->Fill(track->Pt());
104	hMass_rec->Fill(track->M());
105	hCharge_rec->Fill(track->Charge());
106
107	// Read MC particle corresponding to reconstructed tracks
108	AliAODMCParticle mctrack = (AliAODMCParticle) mcarray->At(track->GetLabel());
109	if(!mctrack){
110	printf("Error: no MC track! \n");
111	continue;
112	}
113	//printf("nev=%d i=%d mc=%f rec=%f\n",nev-1,i,mctrack->Pt(),track->Pt());
114	hpt_mc->Fill(mctrack->Pt());
115	hpt_mc_rec->Fill(track->Pt(),mctrack->Pt());
116	hpt_diff->Fill(track->Pt()-mctrack->Pt());
117	hMass_mc->Fill(mctrack->M());
118	hCharge_mc->Fill(mctrack->Charge());
119	}
120	}
121
122	// Plots
123	TCanvas *c = new TCanvas("c","PT spectra",20,20,600,600);
124	c->Divide(3,3);
125	c->cd(1);
126	hpt_primarymc->Draw();
127	hpt_primarymc->GetXaxis()->SetTitle("Pt MC primary particles");
128	c->cd(2);
129	hpt_mc->GetXaxis()->SetTitle("Pt MC particles");
130	hpt_mc->Draw();
131	hpt_mc->SetLineColor(4);
132	c->cd(3);
133	hpt_rec->GetXaxis()->SetTitle("Rec Pt track");
134	hpt_rec->Draw();
135	hpt_rec->SetLineColor(2);
136	c->cd(4);
137	hpt_physprimarymc->Draw();
138	hpt_physprimarymc->GetXaxis()->SetTitle("Pt MC muon track");
139	c->cd(5);
140	hpt_notphysprimarymc->Draw();
141	hpt_notphysprimarymc->GetXaxis()->SetTitle("Pt MC JPsi track");
142	c->cd(7);
143	hpt_mc_rec->Draw("colz");
144	hpt_mc_rec->GetXaxis()->SetTitle("Pt rec track");
145	hpt_mc_rec->GetYaxis()->SetTitle("Pt MC track");
146	c->cd(8);
147	hpt_diff->Draw();
148	hpt_diff->GetXaxis()->SetTitle("Rec-MC Pt");
149
150	TCanvas *c1 = new TCanvas("c1","Mass distributions",40,40,600,600);
151	c1->Divide(2,2);
152	c1->cd(1);
153	hMass_primarymc->Draw();
154	hMass_primarymc->GetXaxis()->SetTitle("Mass MC primary particles");
155	c1->cd(2);
156	hMass_physprimarymc->Draw();
157	hMass_physprimarymc->GetXaxis()->SetTitle("Mass MC muon particles");
158	c1->cd(3);
159	hMass_notphysprimarymc->Draw();
160	hMass_notphysprimarymc->GetXaxis()->SetTitle("Mass MC J/Psi");
161	c1->cd(4);
162	hMass_rec->Draw();
163	hMass_rec->GetXaxis()->SetTitle("Mass Rec tracks");
164
165	TCanvas *c2= new TCanvas("c2","Charge distributions",60,60,600,600);
166	c2->Divide(2,2);
167	c2->cd(1);
168	hCharge_primarymc->Draw();
169	hCharge_primarymc->GetXaxis()->SetTitle("Charge MC primary particles");
170	c2->cd(2);
171	hCharge_physprimmc->Draw();
172	hCharge_physprimmc->GetXaxis()->SetTitle("Charge muon MC");;
173	c2->cd(3);
174	hCharge_notphysprimmc->Draw();
175	hCharge_notphysprimmc->GetXaxis()->SetTitle("Charge J/Psi MC");;
176	c2->cd(4);
177	hCharge_rec->Draw();
178	hCharge_rec->GetXaxis()->SetTitle("Charge rec tracks");
179	}