FairRoot/PandaRoot
Functions
emc_correction_QA.C File Reference

Go to the source code of this file.

Functions

int emc_correction_QA (TString InputFile="emc_complete_QA.root", Int_t test=1, Int_t version=1)
 

Function Documentation

int emc_correction_QA ( TString  InputFile = "emc_complete_QA.root",
Int_t  test = 1,
Int_t  version = 1 
)

Definition at line 7 of file emc_correction_QA.C.

References alpha, c1, c2, c3, c4, c5, cluster_array, cluster_energy, cluster_phi, cluster_theta, cos(), Double_t, PndEmcAbsClusterCalibrator::Energy(), PndEmcCluster::energy(), f, f1, f2, f3, f4, PndEmcCluster::GetEnergyCorrected(), PndMCTrack::GetMomentum(), i, PndEmcClusterCalibrator::MakeEmcClusterCalibrator(), max_energy, mctrack, mctrack_array, Pi, sqrt(), t, threshold, v1, v2, PndEmcAbsClusterCalibrator::Where(), and PndEmcCluster::where().

8 {
9  // test=1
10  TH1F *h_energy= new TH1F("h_energy","Cluster energy",100,0.8,1.2);
11  TH1F *h_energy_corr1= new TH1F("h_energy_corr1","Cluster energy (corrected), (hist)",100,0.8,1.2);
12  TH1F *h_energy_corr2= new TH1F("h_energy_corr2","Cluster energy (corrected), (parametrization)",100,0.8,1.2);
13  TH1F *h_energy_corr_old= new TH1F("h_energy_corr_old","Cluster energy (corrected), (old)",100,0.8,1.2);
14  TH1F *h_energy_shashlyk= new TH1F("h_energy_shashlyk","Cluster energy (shashlyk)",100,0.8,1.2);
15  TH1F *h_energy_shashlyk_corr1= new TH1F("h_energy_shashlyk_corr1","Cluster energy (shashlyk) (corrected), (hist)",100,0.8,1.2);
16  TH1F *h_energy_shashlyk_corr2= new TH1F("h_energy_shashlyk_corr2","Cluster energy (shashlyk) (corrected), (parametrization)",100,0.8,1.2);
17  TH1F *h_energy_shashlyk_corr_old= new TH1F("h_energy_shashlyk_corr_old","Cluster energy (shashlyk) (corrected)",100,0.8,1.2);
18 
19  // test=2
20  TH2F *hEnergyRatioTheta=new TH2F("hEnergyRatioTheta","Energy_reco/Energy_MC",100, 0., 180., 100, 0.8, 1.2);
21  hEnergyRatioTheta->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
22  TH2F *hEnergyRatioEnergy=new TH2F("hEnergyRatioEnergy","Energy_reco/Energy_MC",100, 0., 10., 100, 0.8, 1.2);
23  hEnergyRatioEnergy->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
24  TH2F *hEnergyRatioCorrEnergy1=new TH2F("hEnergyRatioCorrEnergy1","Energy_reco (corrected)/Energy_MC, (hist)",100, 0., 10., 100, 0.8, 1.2);
25  hEnergyRatioCorrEnergy1->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
26  TH2F *hEnergyRatioCorrTheta1=new TH2F("hEnergyRatioCorrTheta1","Energy_reco (corrected)/Energy_MC, (hist)",100, 0., 180., 100, 0.8, 1.2);
27  hEnergyRatioCorrTheta1->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
28  TH2F *hEnergyRatioCorrEnergy2=new TH2F("hEnergyRatioCorrEnergy2","Energy_reco (corrected)/Energy_MC, (parametrization)",100, 0., 10., 100, 0.8, 1.2);
29  hEnergyRatioCorrEnergy2->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
30  TH2F *hEnergyRatioCorrTheta2=new TH2F("hEnergyRatioCorrTheta2","Energy_reco (corrected)/Energy_MC, (parametrization)",100, 0., 180., 100, 0.8, 1.2);
31  hEnergyRatioCorrTheta2->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
32  TH2F *hEnergyRatioCorrEnergyOld=new TH2F("hEnergyRatioCorrEnergyOld","Energy_reco (corrected)/Energy_MC, (old)",100, 0., 10., 100, 0.8, 1.2);
33  hEnergyRatioCorrEnergyOld->GetXaxis()->SetTitle("Cluster Reconstructed Energy (GeV)");
34  TH2F *hEnergyRatioCorrThetaOld=new TH2F("hEnergyRatioCorrThetaOld","Energy_reco (corrected)/Energy_MC, (old)",100, 0., 180., 100, 0.8, 1.2);
35  hEnergyRatioCorrThetaOld->GetXaxis()->SetTitle("Cluster Reconstructed #theta Angle (#circ)");
36 
37  // test=3
38  TH1F *h_mpi0= new TH1F("h_mpi0","pi0 invariant mass",100,0.05,0.2);
39  TH1F *h_mpi0_corr1= new TH1F("h_mpi0_corr1","pi0 invariant mass (corrected) (hist)",100,0.05,0.2);
40  TH1F *h_mpi0_corr2= new TH1F("h_mpi0_corr2","pi0 invariant mass (corrected) (parametrization)",100,0.05,0.2);
41  TH1F *h_mpi0_corr_old= new TH1F("h_mpi0_corr_old","pi0 invariant mass (corrected) (old)",100,0.05,0.2);
42 
43  TFile *f=new TFile(InputFile);
44  TTree *t=(TTree *) f->Get("pndsim") ;
45 
46  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
47  t->SetBranchAddress("MCTrack",&mctrack_array);
48 
49  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
50  t->SetBranchAddress("EmcCluster",&cluster_array);
51 
52  Double_t cluster_energy, cluster_theta, cluster_phi; //position of the cluster
53  TVector3 mc_momentum;
54  Double_t thetaMC, phiMC, energyMC;
55  Double_t max_energy;
56 
57  // Here the cluster calibrator is set
58  // first parameter is method (1-histogram, 2-parametrization)
59  // second parameter is version (defined by geometry)
60  PndEmcAbsClusterCalibrator * calibrator1= PndEmcClusterCalibrator::MakeEmcClusterCalibrator(1, version);
61  PndEmcAbsClusterCalibrator * calibrator2= PndEmcClusterCalibrator::MakeEmcClusterCalibrator(2, version);
62 
63  if ((test==1)||(test==2))
64  {
65  for (Int_t j=0; j<t->GetEntriesFast(); j++)
66  {
67  t->GetEntry(j);
68  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
69  mc_momentum=mctrack->GetMomentum();
70  energyMC=mc_momentum.Mag();
71  thetaMC=mc_momentum.Theta()*180./TMath::Pi();
72  max_energy=0;
73  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
74  {
75  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
76  cluster_energy=cluster->energy();
77 
78  if (cluster_energy>max_energy)
79  {
80  max_energy=cluster_energy;
81  TVector3 cluster_pos=cluster->where();
82  cluster_theta=cluster_pos.Theta()*180./TMath::Pi();
83  Double_t energy_corr1 = calibrator1->Energy(cluster);
84  Double_t energy_corr2 = calibrator2->Energy(cluster);
85  Double_t energy_corr_old = cluster->GetEnergyCorrected();
86  }
87 
88  Double_t e_ratio=max_energy/energyMC;
89  Double_t e_ratio_corr1=energy_corr1/energyMC;
90  Double_t e_ratio_corr2=energy_corr2/energyMC;
91  Double_t e_ratio_corr_old=energy_corr_old/energyMC;
92 
93  if (cluster_theta<5.)
94  {
95  h_energy_shashlyk->Fill(max_energy);
96  h_energy_shashlyk_corr1->Fill(energy_corr1);
97  h_energy_shashlyk_corr2->Fill(energy_corr2);
98  h_energy_shashlyk_corr_old->Fill(energy_corr_old);
99  } else {
100  h_energy->Fill(max_energy);
101  h_energy_corr1->Fill(energy_corr1);
102  h_energy_corr2->Fill(energy_corr2);
103  h_energy_corr_old->Fill(energy_corr_old);
104  }
105 
106  hEnergyRatioEnergy->Fill( energyMC, e_ratio);
107  hEnergyRatioTheta->Fill( thetaMC, e_ratio);
108  hEnergyRatioCorrEnergy1->Fill(energyMC, e_ratio_corr1);
109  hEnergyRatioCorrTheta1->Fill(thetaMC, e_ratio_corr1);
110  hEnergyRatioCorrEnergy2->Fill(energyMC, e_ratio_corr2);
111  hEnergyRatioCorrTheta2->Fill(thetaMC, e_ratio_corr2);
112  hEnergyRatioCorrEnergyOld->Fill(energyMC, e_ratio_corr_old);
113  hEnergyRatioCorrThetaOld->Fill(thetaMC, e_ratio_corr_old);
114  }
115  }
116  }
117 
118  if (test==3)
119  {
120  double threshold=0.02;
121  TVector3 v1,v2;
122  TVector3 v1_corr1,v2_corr1;
123  TVector3 v1_corr2,v2_corr2;
124  Double_t energy1, energy2, energy1_corr1, energy2_corr1;
125  Double_t energy1_corr2, energy2_corr2, energy1_corr_old, energy2_corr_old;
126  // Cluster energy
127  for (Int_t j=0; j< t->GetEntriesFast(); j++)
128  {
129  t->GetEntry(j);
130  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
131  {
132  PndEmcCluster *cluster1=(PndEmcCluster*)cluster_array->At(i);
133  energy1=cluster1->energy();
134  if (energy1<threshold) continue;
135  v1=cluster1->where();
136  energy1_corr1 = calibrator1->Energy(cluster1);
137  v1_corr1 = calibrator1->Where(cluster1);
138  energy1_corr2 = calibrator2->Energy(cluster1);
139  v1_corr2 = calibrator2->Where(cluster1);
140  energy1_corr_old = cluster1->GetEnergyCorrected();
141 
142  for (Int_t k=i; k<cluster_array->GetEntriesFast(); k++)
143  {
144  PndEmcCluster *cluster2=(PndEmcCluster*)cluster_array->At(k);
145  if (cluster2->energy()<threshold) continue;
146  energy2=cluster2->energy();
147  energy2_corr1 = calibrator1->Energy(cluster2);
148  energy2_corr2 = calibrator2->Energy(cluster2);
149  energy2_corr_old = cluster2->GetEnergyCorrected();
150  v2=cluster2->where();
151  v2_corr1 = calibrator1->Where(cluster2);
152  v2_corr2 = calibrator2->Where(cluster2);
153  double alpha=v1.Angle(v2);
154  double invMass=sqrt(2*energy1*energy2*(1-cos(alpha)));
155  h_mpi0->Fill(invMass);
156  double alpha_corr1=v1_corr1.Angle(v2_corr1);
157  double invMass_corr1=sqrt(2*energy1_corr1*energy2_corr1*(1-cos(alpha_corr1)));
158  h_mpi0_corr1->Fill(invMass_corr1);
159  double alpha_corr2=v1_corr2.Angle(v2_corr2);
160  double invMass_corr2=sqrt(2*energy1_corr2*energy2_corr2*(1-cos(alpha_corr2)));
161  h_mpi0_corr2->Fill(invMass_corr2);
162  double invMass_corr_old=sqrt(2*energy1_corr_old*energy2_corr_old*(1-cos(alpha)));
163  h_mpi0_corr_old->Fill(invMass_corr_old);
164  }
165  }
166  }
167  }
168 
169 
170 if (test==1)
171 {
172  TCanvas* c1 = new TCanvas("c1", "Energy reconstruction vs Energy", 100, 100, 800, 800);
173  c1->Divide(2,2);
174  c1->cd(1);
175  hEnergyRatioEnergy->Draw("contz");
176  TLine *l=new TLine(0.,1.,10.,1.0);
177  l->SetLineColor(2);
178  l->SetLineWidth(2);
179  l->Draw();
180 
181  c1->cd(2);
182  hEnergyRatioCorrEnergy1->Draw("contz");
183  TLine *l=new TLine(0.,1.,10.,1.0);
184  l->SetLineColor(2);
185  l->SetLineWidth(2);
186  l->Draw();
187 
188  c1->cd(3);
189  hEnergyRatioCorrEnergy2->Draw("contz");
190  TLine *l=new TLine(0.,1.,10.,1.0);
191  l->SetLineColor(2);
192  l->SetLineWidth(2);
193  l->Draw();
194 
195 // c1->cd(4);
196 // hEnergyRatioCorrEnergyOld->Draw("contz");
197 // TLine *l=new TLine(0.,1.,10.,1.0);
198 // l->SetLineColor(2);
199 // l->SetLineWidth(2);
200 // l->Draw();
201 
202  TCanvas* c2 = new TCanvas("c2", "Energy reconstruction vs Theta", 100, 100, 800, 800);
203  c2->Divide(2,2);
204  c2->cd(1);
205  hEnergyRatioTheta->Draw("contz");
206  TLine *l=new TLine(0.,1.,180.,1.0);
207  l->SetLineColor(2);
208  l->SetLineWidth(2);
209  l->Draw();
210  c2->cd(2);
211  hEnergyRatioCorrTheta1->Draw("contz");
212  TLine *l=new TLine(0.,1.,180.,1.0);
213  l->SetLineColor(2);
214  l->SetLineWidth(2);
215  l->Draw();
216  c2->cd(3);
217  hEnergyRatioCorrTheta2->Draw("contz");
218  TLine *l=new TLine(0.,1.,180.,1.0);
219  l->SetLineColor(2);
220  l->SetLineWidth(2);
221  l->Draw();
222 // c2->cd(4);
223 // hEnergyRatioCorrThetaOld->Draw("contz");
224 // TLine *l=new TLine(0.,1.,180.,1.0);
225 // l->SetLineColor(2);
226 // l->SetLineWidth(2);
227 // l->Draw();
228 }
229 
230 if (test==2)
231 {
232  TCanvas* c3 = new TCanvas("c3", "Energy (TS)", 100, 100, 800, 800);
233  c3->Divide(2,2);
234  c3->cd(1);
235  h_energy->Draw();
236  c3->cd(2);
237  h_energy_corr1->Draw();
238  c3->cd(3);
239  h_energy_corr2->Draw();
240 // c3->cd(4);
241 // h_energy_corr_old->Draw();
242 
243  TCanvas* c4 = new TCanvas("c4", "Energy (shashlyk)", 100, 100, 800, 800);
244  c4->Divide(2,2);
245  c4->cd(1);
246  h_energy_shashlyk->Draw();
247  c4->cd(2);
248  h_energy_shashlyk_corr1->Draw();
249  c4->cd(3);
250  h_energy_shashlyk_corr2->Draw();
251 // c4->cd(4);
252 // h_energy_shashlyk_corr_old->Draw();
253 }
254 if (test==3)
255 {
256  gStyle->SetOptFit(kTRUE);
257 
258  double pars[5]={800,0.13,0.01,100,-500};
259  TF1 *f1 = new TF1("f1","gaus(0)+pol1(3)",0.08,0.18);
260  TF1 *f2 = new TF1("f2","gaus(0)+pol1(3)",0.08,0.18);
261  TF1 *f3 = new TF1("f3","gaus(0)+pol1(3)",0.08,0.18);
262  TF1 *f4 = new TF1("f4","gaus(0)+pol1(3)",0.08,0.18);
263  f1->SetParameters(pars);
264  f2->SetParameters(pars);
265  f3->SetParameters(pars);
266  f4->SetParameters(pars);
267 
268  TCanvas* c5 = new TCanvas("c5", "pi0 mass", 100, 100, 800, 800);
269  c5->Divide(2,2);
270  c5->cd(1);
271  h_mpi0->Draw();
272  h_mpi0->Fit(f1,"R");
273  c5->cd(2);
274  h_mpi0_corr1->Draw();
275  h_mpi0_corr1->Fit(f2,"R");
276  c5->cd(3);
277  h_mpi0_corr2->Draw();
278  h_mpi0_corr2->Fit(f3,"R");
279 // c5->cd(4);
280 // h_mpi0_corr_old->Draw();
281 // h_mpi0_corr_old->Fit(f4,"R");
282 }
283  return 0;
284 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
c5
c5
Definition: plot_dirc.C:75
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
PndEmcCluster::where
TVector3 where() const
Definition: PndEmcCluster.cxx:111
c4
c4
Definition: plot_dirc.C:71
i
Int_t i
Definition: run_full.C:25
f1
TF1 * f1
Definition: reco_analys2.C:50
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
c2
c2
Definition: plot_dirc.C:39
f4
TFile * f4
Definition: ReadMesDecCorrRes.C:9
PndEmcAbsClusterCalibrator::Energy
virtual Double_t Energy(PndEmcCluster *clust, Int_t pid=22)=0
cluster_theta
Double_t cluster_theta
Definition: anal_hit_digi_cluster_fwendcap.C:31
f3
TFile * f3
Definition: ReadMesDecCorrRes.C:8
cluster_array
TClonesArray * cluster_array
Definition: anal_hit_digi_cluster_fwendcap.C:27
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
cluster_energy
Double_t cluster_energy
Definition: anal_hit_digi_cluster_fwendcap.C:31
f
TFile * f
Definition: bump_analys.C:12
PndEmcClusterCalibrator::MakeEmcClusterCalibrator
static PndEmcAbsClusterCalibrator * MakeEmcClusterCalibrator(Int_t method, Int_t version=1)
Definition: PndEmcClusterCalibrator.cxx:58
cluster_phi
Double_t cluster_phi
Definition: anal_hit_digi_cluster_fwendcap.C:31
PndEmcCluster
a cluster (group of neighboring crystals) of hit emc crystals
Definition: PndEmcCluster.h:29
c1
c1
Definition: plot_dirc.C:35
threshold
double threshold
Definition: outdated/dc1/stt/runreco.C:105
c3
c3
Definition: plot_dirc.C:50
PndEmcAbsClusterCalibrator::Where
virtual TVector3 Where(PndEmcCluster *clust, Int_t pid=22)=0
max_energy
double max_energy
Definition: macro/detectors/emc/dedicated/reco_analys.C:27
PndEmcCluster::GetEnergyCorrected
Double_t GetEnergyCorrected() const
Definition: PndEmcCluster.cxx:379
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
v1
TVector3 v1
Definition: bump_analys.C:40
PndEmcCluster::energy
virtual Double_t energy() const
Definition: PndEmcCluster.cxx:79
v2
TVector3 v2
Definition: bump_analys.C:40
t
TTree * t
Definition: bump_analys.C:13
alpha
double alpha
Definition: f_Init.h:9
mctrack_array
TClonesArray * mctrack_array
Definition: anal_hit_digi_cluster_fwendcap.C:11
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
Generated on Wed Apr 3 2019 10:03:13 for FairRoot/PandaRoot by   doxygen 1.8.5