27 TChain *
c=
new TChain(
"pndsim");
29 if (InputFile1.Contains(
".root"))
31 TObjArray* tmp=InputFile1.Tokenize(
"_");
32 TObjString *os=(TObjString *)tmp->At(2);
33 particle=os->GetString();
34 TFile *
f1=
new TFile(InputFile1);
37 std::cout<<
"Input file "<<InputFile1<<
" does not exists"<<std::endl;
42 else if (InputFile1.Contains(
".txt"))
45 ifstream
infile(InputFile1.Data(), ios::in);
47 while (getline(
infile,file_name,
'\n'))
49 if (std::string::npos != file_name.find(
".root"))
51 TFile *
f2=
new TFile(file_name.c_str());
54 std::cout<<
"Input file "<<file_name<<
" does not exists"<<std::endl;
57 c->Add(file_name.c_str());
60 TObjArray* tmp=
TString(file_name).Tokenize(
"_");
61 TObjString *os=(TObjString *)tmp->At(2);
62 particle=os->GetString();
69 std::cout<<
"Wrong input file: "<<InputFile1<<std::endl;
73 TString OutputFile=
"emc_correction_hist_";
74 OutputFile = OutputFile +particle+
"_";
76 OutputFile +=
".root";
84 TH2F hisEnergyRatioBarrel, hisEnergyRatioFwd, hisEnergyRatioBwd, hisEnergyRatioShashlyk;
85 TH2F hisThetaDiffBarrel, hisThetaDiffFwd, hisThetaDiffBwd, hisThetaDiffShashlyk;
88 hisEnergyRatioBarrel.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
89 hisEnergyRatioBarrel.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
92 hisThetaDiffBarrel.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
93 hisThetaDiffBarrel.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
96 hisEnergyRatioFwd.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
97 hisEnergyRatioFwd.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
100 hisThetaDiffFwd.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
101 hisThetaDiffFwd.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
104 hisEnergyRatioBwd.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
105 hisEnergyRatioBwd.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
108 hisThetaDiffBwd.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
109 hisThetaDiffBwd.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
112 hisEnergyRatioShashlyk.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
113 hisEnergyRatioShashlyk.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
116 hisThetaDiffShashlyk.GetXaxis()->SetTitle(
"Cluster Reconstructed Photon Energy (GeV)");
117 hisThetaDiffShashlyk.GetYaxis()->SetTitle(
"Cluster Reconstructed #theta Photon Angle (#circ)");
121 TH1F hisEnergyRatio1[100][100];
122 TH1F hisEnergyRatio2[100][100];
123 TH1F hisEnergyRatio3[100][100];
124 TH1F hisEnergyRatio4[100][100];
127 TH1F hisThetaDiff1[100][100];
128 TH1F hisThetaDiff2[100][100];
129 TH1F hisThetaDiff3[100][100];
130 TH1F hisThetaDiff4[100][100];
138 sprintf(cmd,
"hisEnergyRatio1_%i_%i",
i,j);
139 hisEnergyRatio1[
i][j].SetName(cmd);
140 hisEnergyRatio1[
i][j].SetTitle(cmd);
141 hisEnergyRatio1[
i][j].SetBins(100,0.,2.);
142 sprintf(cmd,
"hisThetaDiff1_%i_%i",
i,j);
143 hisThetaDiff1[
i][j].SetName(cmd);
144 hisThetaDiff1[
i][j].SetTitle(cmd);
145 hisThetaDiff1[
i][j].SetBins(100,-5.,5);
146 hisEnergyRatioBarrel.SetBinContent(
i+1,j+1,1.);
147 hisThetaDiffBarrel.SetBinContent(
i+1,j+1,0.);
151 sprintf(cmd,
"hisEnergyRatio2_%i_%i",
i,j);
152 hisEnergyRatio2[
i][j].SetName(cmd);
153 hisEnergyRatio2[
i][j].SetTitle(cmd);
154 hisEnergyRatio2[
i][j].SetBins(100,0.,2.);
155 sprintf(cmd,
"hisThetaDiff2_%i_%i",
i,j);
156 hisThetaDiff2[
i][j].SetName(cmd);
157 hisThetaDiff2[
i][j].SetTitle(cmd);
158 hisThetaDiff2[
i][j].SetBins(100,-5.,5);
159 hisEnergyRatioFwd.SetBinContent(
i+1,j+1,1.);
160 hisThetaDiffFwd.SetBinContent(
i+1,j+1,0.);
164 sprintf(cmd,
"hisEnergyRatio3_%i_%i",
i,j);
165 hisEnergyRatio3[
i][j].SetName(cmd);
166 hisEnergyRatio3[
i][j].SetTitle(cmd);
167 hisEnergyRatio3[
i][j].SetBins(100,0.,2.);
168 sprintf(cmd,
"hisThetaDiff3_%i_%i",
i,j);
169 hisThetaDiff3[
i][j].SetName(cmd);
170 hisThetaDiff3[
i][j].SetTitle(cmd);
171 hisThetaDiff3[
i][j].SetBins(100,-5.,5);
172 hisEnergyRatioBwd.SetBinContent(
i+1,j+1,1.);
173 hisThetaDiffBwd.SetBinContent(
i+1,j+1,0.);
177 sprintf(cmd,
"hisEnergyRatio4_%i_%i",
i,j);
178 hisEnergyRatio4[
i][j].SetName(cmd);
179 hisEnergyRatio4[
i][j].SetTitle(cmd);
180 hisEnergyRatio4[
i][j].SetBins(100,0.,2.);
181 sprintf(cmd,
"hisThetaDiff4_%i_%i",
i,j);
182 hisThetaDiff4[
i][j].SetName(cmd);
183 hisThetaDiff4[
i][j].SetTitle(cmd);
184 hisThetaDiff4[
i][j].SetBins(100,-5.,5);
185 hisEnergyRatioShashlyk.SetBinContent(
i+1,j+1,1.);
186 hisThetaDiffShashlyk.SetBinContent(
i+1,j+1,0.);
193 TClonesArray*
track_array=
new TClonesArray(
"PndMCTrack");
194 c->SetBranchAddress(
"MCTrack",&track_array);
195 TClonesArray*
cluster_array=
new TClonesArray(
"PndEmcCluster");
196 c->SetBranchAddress(
"EmcCluster",&cluster_array);
201 for (Int_t j=0; j< c->GetEntries(); j++)
204 cout<<
"Event "<<j<<endl;
212 thMC = mc_momentum.Theta()*(180./
TMath::Pi());
213 phiMC = mc_momentum.Phi()*(180./
TMath::Pi());
218 if (cluster_array->GetEntriesFast()>0)
220 Int_t idWithHighestEnergy = 0;
223 for (Int_t
i=0;
i<cluster_array->GetEntriesFast();
i++)
227 cluster_energy=cluster->
energy();
229 if (cluster_energy>highestEnergy)
231 idWithHighestEnergy =
i;
238 TVector3 cluster_pos=cluster->
where();
239 cluster_theta=cluster_pos.Theta()*180./
TMath::Pi();
240 cluster_phi=cluster_pos.Phi()*180./
TMath::Pi();
241 cluster_energy=cluster->
energy();
247 if (cluster_theta >= 141. && cluster_theta < 147.)
continue;
250 if (thetaBin<0 || energyBin<0)
252 cout<<
"!!!!!!!!!!!!!!!!!!!!!!"<<endl;
253 cout<<
"Module="<<module<<endl;
254 cout <<
"theta_reco =" << cluster_theta <<
"/ energy_reco =" << cluster_energy << endl;
255 cout <<
"thetaBin= "<<thetaBin<<
" & energyBin= "<<energyBin<<endl;
259 if (cluster_theta > (thMC-5) && cluster_theta <=(thMC+5) ){
261 en_div = cluster_energy/enMC;
265 if (
fabs(phi_diff)<2.5 ||
fabs(phi_diff+360)<2.5 ||
fabs(phi_diff-360)<2.5){
272 if ((
fabs(phiMC) <5.)||
fabs(phiMC-180.)<5.||
fabs(phiMC-360.)<5.)
break;
273 if (en_div > 0.7 && en_div <1.3){
274 hisEnergyRatio1[energyBin][thetaBin].Fill(en_div);
275 hisThetaDiff1[energyBin][thetaBin].Fill(theta_diff);
281 if (cluster_theta <6.)
break;
282 if (en_div > 0.7 && en_div <1.3){
283 hisEnergyRatio2[energyBin][thetaBin].Fill(en_div);
284 hisThetaDiff2[energyBin][thetaBin].Fill(theta_diff);
290 if (en_div > 0.7 && en_div <1.3){
291 hisEnergyRatio3[energyBin][thetaBin].Fill(en_div);
292 hisThetaDiff3[energyBin][thetaBin].Fill(theta_diff);
298 if (en_div > 0.7 && en_div <1.3){
299 hisEnergyRatio4[energyBin][thetaBin].Fill(en_div);
300 hisThetaDiff4[energyBin][thetaBin].Fill(theta_diff);
320 TFile fout(OutputFile,
"recreate");
321 Double_t mean_th_dif, mean_en_ratio;
333 mean_th_dif = hisThetaDiff1[
i][j].GetMean();
336 hisThetaDiffBarrel.SetBinContent(
i+1,j+1,mean_th_dif);
341 mean_en_ratio = hisEnergyRatio1[
i][j].GetMean();
344 hisEnergyRatioBarrel.SetBinContent(
i+1,j+1,mean_en_ratio);
354 mean_th_dif = hisThetaDiff2[
i][j].GetMean();
357 hisThetaDiffFwd.SetBinContent(
i+1,j+1,mean_th_dif);
362 mean_en_ratio = hisEnergyRatio2[
i][j].GetMean();
365 hisEnergyRatioFwd.SetBinContent(
i+1,j+1,mean_en_ratio);
375 mean_th_dif = hisThetaDiff3[
i][j].GetMean();
378 hisThetaDiffBwd.SetBinContent(
i+1,j+1,mean_th_dif);
383 mean_en_ratio = hisEnergyRatio3[
i][j].GetMean();
386 hisEnergyRatioBwd.SetBinContent(
i+1,j+1,mean_en_ratio);
396 mean_th_dif = hisThetaDiff4[
i][j].GetMean();
399 hisThetaDiffShashlyk.SetBinContent(
i+1,j+1,mean_th_dif);
405 mean_en_ratio = hisEnergyRatio4[
i][j].GetMean();
408 hisEnergyRatioShashlyk.SetBinContent(
i+1,j+1,mean_en_ratio);
412 hisEnergyRatioBarrel.Write();
413 hisThetaDiffBarrel.Write();
414 hisEnergyRatioFwd.Write();
415 hisThetaDiffFwd.Write();
416 hisEnergyRatioBwd.Write();
417 hisThetaDiffBwd.Write();
418 hisEnergyRatioShashlyk.Write();
419 hisThetaDiffShashlyk.Write();
427 TH1F *
h1=(TH1F *)hisEnergyRatio1[
i][j].Clone();
429 TH1F *h1=(TH1F *)hisThetaDiff1[
i][j].Clone();
434 TH1F *h1=(TH1F *)hisEnergyRatio2[
i][j].Clone();
436 TH1F *h1=(TH1F *)hisThetaDiff2[
i][j].Clone();
441 TH1F *h1=(TH1F *)hisEnergyRatio3[
i][j].Clone();
443 TH1F *h1=(TH1F *)hisThetaDiff3[
i][j].Clone();
448 TH1F *h1=(TH1F *)hisEnergyRatio4[
i][j].Clone();
450 TH1F *h1=(TH1F *)hisThetaDiff4[
i][j].Clone();
461 printf(
"RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
Double_t thetaIntervalsBwd[]
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
Short_t GetModule() const
Double_t thetaIntervalsBarrel[]
TVector3 GetMomentum() const
TLorentzVector Get4Momentum() const
TClonesArray * cluster_array
Int_t nrThetaIntervals[4]
friend F32vec4 fabs(const F32vec4 &a)
a cluster (group of neighboring crystals) of hit emc crystals
Double_t energyIntervals[]
Double_t thetaIntervalsShashlyk[]
Int_t GetEnergyBin(Double_t val, Int_t module)
cout<<"will loop over "<< t-> GetEntries()
virtual Double_t energy() const
TClonesArray * track_array
TFile infile("dedx_out.root","READ")
Double_t thetaIntervalsFwd[]
Int_t GetThetaBin(Double_t val, Int_t module)