FairRoot/PandaRoot
Functions | Variables
emc_correction_parametrization.C File Reference

Go to the source code of this file.

Functions

int emc_correction_parametrization (Int_t version, TString InputFile1, TString particle="gamma", Bool_t debug=false, Bool_t useStoredHistos=true)
 
Double_t FitFunction1 (Double_t *x, Double_t *par)
 
Double_t FitFunction2 (Double_t *x, Double_t *par)
 
Int_t GetThetaBin (Double_t val, Int_t range_set)
 
Int_t GetEnergyBin (Double_t val, Int_t range_set)
 

Variables

Double_t energyIntervalsBarrelLow [] = {0,0.03, 0.05, 0.07, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}
 
Double_t thetaIntervalsBarrelLow [] = {22.,23.,24.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100.,105.,110.,115.,120.,125.,130.,135.,137., 139., 140., 141.}
 
Double_t energyIntervalsBarrelHigh [] = {1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}
 
Double_t thetaIntervalsBarrelHigh [] = {22.,23.,24.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100.,105.,110.,115.,120.,125.,130.,135.,137., 139., 140., 141.}
 
Double_t energyIntervalsFwd [] = {0,0.03, 0.05, 0.07, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}
 
Double_t thetaIntervalsFwd [] = {5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,20.,21.,22.}
 
Double_t energyIntervalsBwd [] = {0,0.03,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0}
 
Double_t thetaIntervalsBwd [] = {147.,148.,149.,150.,151.,152.,153.,154.,155.,156.,157.,158.,159.,160.,161.,162.,163.,164.,166.,168.,170.,172.}
 
Double_t energyIntervalsShashlyk [] = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}
 
Int_t nrEnergyIntervals [5]
 
Int_t nrThetaIntervals [4]
 

Function Documentation

int emc_correction_parametrization ( Int_t  version,
TString  InputFile1,
TString  particle = "gamma",
Bool_t  debug = false,
Bool_t  useStoredHistos = true 
)

Definition at line 24 of file emc_correction_parametrization.C.

References Bool_t, c, c1, c2, c3, c4, c5, cluster_array, cluster_energy, cluster_phi, cluster_theta, ctime, Double_t, PndEmcCluster::energy(), energyIntervalsBarrelHigh, energyIntervalsBarrelLow, energyIntervalsBwd, energyIntervalsFwd, energyIntervalsShashlyk, f1, f2, fabs(), file_name, FitFunction1(), FitFunction2(), PndMCTrack::Get4Momentum(), GetEnergyBin(), GetEntries(), PndEmcCluster::GetModule(), PndMCTrack::GetMomentum(), GetThetaBin(), i, infile(), nrEnergyIntervals, nrThetaIntervals, p1, p2, p4mom, particle, phi_diff, Pi, printf(), rtime, PndEmcClusterCalibrationParObject::SetCalibrationPar(), thetaIntervalsBarrelHigh, thetaIntervalsBarrelLow, thetaIntervalsBwd, thetaIntervalsFwd, timer, track, track_array, TString, and PndEmcCluster::where().

25 {
26 TStopwatch timer;
27 timer.Start();
28 
29  // If exists intermidiate file with histograms and useStoredHistos==true
30  // perform the fit only
31  // If useStoredHistos==true but file "correction_histos.root" does not exist
32  // then it is created
33  Bool_t storeHistos=false;
34  TFile *f_histos;
35  if (useStoredHistos)
36  {
37  f_histos=new TFile("correction_histos_fit.root");
38  if (f_histos->IsZombie())
39  {
40  storeHistos=true;
41  }
42  f_histos->Close();
43  }
44 
45  TH2F *hisEnergyRatioBarrelLow, *hisEnergyRatioBarrelHigh;
46  TH2F *hisEnergyRatioFwd, *hisEnergyRatioBwd;
47  TH1F *hisEnergyRatioShashlyk;
48 
49  TString OutputFile="emc_correction_par_";
50  OutputFile = OutputFile +particle+"_";
51  OutputFile +=version;
52  OutputFile += ".root";
53 
54 
55  if (!useStoredHistos||storeHistos)
56  {
57  TChain *c=new TChain("pndsim");
58  if (InputFile1.Contains(".root"))
59  {
60  TFile *f1=new TFile(InputFile1);
61  if (f1->IsZombie())
62  {
63  std::cout<<"Input file "<<InputFile1<<" does not exists"<<std::endl;
64  abort();
65  }
66  c->Add(InputFile1);
67  }
68  else if (InputFile1.Contains(".txt"))
69  {
70  string file_name;
71  ifstream infile(InputFile1.Data(), ios::in);
72 
73  while (getline(infile,file_name, '\n'))
74  {
75  if (std::string::npos != file_name.find(".root"))
76  {
77  TFile *f2=new TFile(file_name.c_str());
78  if (f2->IsZombie())
79  {
80  std::cout<<"Input file "<<file_name<<" does not exists"<<std::endl;
81  abort();
82  }
83  c->Add(file_name.c_str());
84  }
85  }
86  }
87  else
88  {
89  std::cout<<"Wrong input file: "<<InputFile1<<std::endl;
90  abort();
91  }
92 
93  nrEnergyIntervals[0]=sizeof(energyIntervalsBarrelLow)/sizeof(Double_t);
94  nrEnergyIntervals[1]=sizeof(energyIntervalsBarrelHigh)/sizeof(Double_t);
95  nrEnergyIntervals[2]=sizeof(energyIntervalsFwd)/sizeof(Double_t);
96  nrEnergyIntervals[3]=sizeof(energyIntervalsBwd)/sizeof(Double_t);
97  nrEnergyIntervals[4]=sizeof(energyIntervalsShashlyk)/sizeof(Double_t);
98 
99  nrThetaIntervals[0]=sizeof(thetaIntervalsBarrelLow)/sizeof(Double_t);
100  nrThetaIntervals[1]=sizeof(thetaIntervalsBarrelHigh)/sizeof(Double_t);
101  nrThetaIntervals[2]=sizeof(thetaIntervalsFwd)/sizeof(Double_t);
102  nrThetaIntervals[3]=sizeof(thetaIntervalsBwd)/sizeof(Double_t);
103 
104  hisEnergyRatioBarrelLow=new TH2F("hisEnergyRatioBarrelLow","Ene_MC/Ene_reco: GetMean()",(nrEnergyIntervals[0]-1),energyIntervalsBarrelLow,(nrThetaIntervals[0]-1),thetaIntervalsBarrelLow);
105  hisEnergyRatioBarrelLow->GetXaxis()->SetTitle("Cluster Reconstructed Photon Energy (GeV)");
106  hisEnergyRatioBarrelLow->GetYaxis()->SetTitle("Cluster Reconstructed #theta Photon Angle (#circ)");
107 
108  hisEnergyRatioBarrelHigh=new TH2F("hisEnergyRatioBarrelHigh","Ene_MC/Ene_reco: GetMean()",(nrEnergyIntervals[1]-1),energyIntervalsBarrelHigh,(nrThetaIntervals[1]-1),thetaIntervalsBarrelHigh);
109  hisEnergyRatioBarrelHigh->GetXaxis()->SetTitle("Cluster Reconstructed Photon Energy (GeV)");
110  hisEnergyRatioBarrelHigh->GetYaxis()->SetTitle("Cluster Reconstructed #theta Photon Angle (#circ)");
111 
112  hisEnergyRatioFwd=new TH2F("hisEnergyRatioFwd","Ene_MC/Ene_reco: GetMean()",(nrEnergyIntervals[2]-1),energyIntervalsFwd,(nrThetaIntervals[2]-1),thetaIntervalsFwd);
113  hisEnergyRatioFwd->GetXaxis()->SetTitle("Cluster Reconstructed Photon Energy (GeV)");
114  hisEnergyRatioFwd->GetYaxis()->SetTitle("Cluster Reconstructed #theta Photon Angle (#circ)");
115 
116  hisEnergyRatioBwd=new TH2F("hisEnergyRatioBwd","Ene_MC/Ene_reco: GetMean()",(nrEnergyIntervals[3]-1),energyIntervalsBwd,(nrThetaIntervals[3]-1),thetaIntervalsBwd);
117  hisEnergyRatioBwd->GetXaxis()->SetTitle("Cluster Reconstructed Photon Energy (GeV)");
118  hisEnergyRatioBwd->GetYaxis()->SetTitle("Cluster Reconstructed #theta Photon Angle (#circ)");
119 
120  hisEnergyRatioShashlyk=new TH1F("hisEnergyRatioShashlyk","Ene_MC/Ene_reco _Shashlyk_: GetMean()",(nrEnergyIntervals[4]-1),energyIntervalsShashlyk);
121  hisEnergyRatioShashlyk->GetXaxis()->SetTitle("Cluster Reconstructed Photon Energy (GeV)");
122 
123  // energy: "E_reco / E_MC"
124  TH1F hisEnergyRatio1[50][50]; // barrel low
125  TH1F hisEnergyRatio2[50][50]; // barrel high
126  TH1F hisEnergyRatio3[50][50]; // forward endcap
127  TH1F hisEnergyRatio4[50][50]; // backward endcap
128  TH1F hisEnergyRatio5[50]; // shashlyk
129 
130  for (Int_t i=0; i<(nrEnergyIntervals[0]-1); i++)
131  {
132  for (Int_t j=0; j<(nrThetaIntervals[0]-1); j++)
133  {
134  hisEnergyRatio1[i][j].SetBins(100,0.,2.);
135  hisEnergyRatioBarrelLow->SetBinContent(i+1,j+1,1.);
136  }
137  }
138 
139  for (Int_t i=0; i<(nrEnergyIntervals[1]-1); i++)
140  {
141  for (Int_t j=0; j<(nrThetaIntervals[1]-1); j++)
142  {
143  hisEnergyRatio2[i][j].SetBins(100,0.,2.);
144  hisEnergyRatioBarrelHigh->SetBinContent(i+1,j+1,1.);
145  }
146  }
147 
148  for (Int_t i=0; i<(nrEnergyIntervals[2]-1); i++)
149  {
150  for (Int_t j=0; j<(nrThetaIntervals[2]-1); j++)
151  {
152  hisEnergyRatio3[i][j].SetBins(100,0.,2.);
153  hisEnergyRatioFwd->SetBinContent(i+1,j+1,1.);
154  }
155  }
156 
157  for (Int_t i=0; i<(nrEnergyIntervals[3]-1); i++)
158  {
159  for (Int_t j=0; j<(nrThetaIntervals[3]-1); j++)
160  {
161  hisEnergyRatio4[i][j].SetBins(100,0.,2.);
162  hisEnergyRatioBwd->SetBinContent(i+1,j+1,1.);
163  }
164  }
165 
166  for (Int_t i=0; i<(nrEnergyIntervals[4]-1); i++)
167  {
168  hisEnergyRatio5[i].SetBins(100,0.,2.);
169  hisEnergyRatioShashlyk->SetBinContent(i+1,j+1,1.);
170  }
171 
172  TClonesArray* track_array=new TClonesArray("PndMCTrack");
173  c->SetBranchAddress("MCTrack",&track_array);
174  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
175  c->SetBranchAddress("EmcCluster",&cluster_array);
176 
177  Double_t cluster_energy, cluster_theta, cluster_phi; //position of the cluster
178  Double_t thMC, phiMC, enMC;
179 
180 timer.Stop();
181 Double_t rtime = timer.RealTime();
182 Double_t ctime = timer.CpuTime();
183 timer.Continue();
184 printf("Intitialization done RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
185 
186 cout<<"c->GetEntries()="<<c->GetEntries()<<endl;
187  for (Int_t j=0; j< c->GetEntries(); j++)
188  {
189 
190  if (j%1000==0)
191  {
192  cout<<"Event "<<j<<endl;
193  }
194  c->GetEntry(j);
195 
196  // Select cluster of highest energy
197  if (cluster_array->GetEntriesFast()>0)
198  {
199  PndMCTrack *track=(PndMCTrack*)track_array->At(0);
200  TLorentzVector p4mom=track->Get4Momentum();
201  TVector3 mc_momentum=track->GetMomentum();
202 
203  thMC = mc_momentum.Theta()*(180./TMath::Pi());
204  phiMC = mc_momentum.Phi()*(180./TMath::Pi());
205 
206  enMC = p4mom.E();
207 
208  Int_t idWithHighestEnergy = 0;
209  Double_t highestEnergy = -1.;
210 
211  for (Int_t i=0; i<cluster_array->GetEntriesFast(); i++)
212  {
213  PndEmcCluster *cluster;
214  cluster=(PndEmcCluster*)cluster_array->At(i);
215  cluster_energy=cluster->energy();
216 
217  if (cluster_energy>highestEnergy)
218  {
219  idWithHighestEnergy = i;
220  highestEnergy = cluster_energy;
221  }
222  }
223 
224  // Lets analyze that cluster!
225  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(idWithHighestEnergy);
226  TVector3 cluster_pos=cluster->where();
227  cluster_theta=cluster_pos.Theta()*180./TMath::Pi();
228  cluster_phi=cluster_pos.Phi()*180./TMath::Pi();
229  cluster_energy=highestEnergy;
230  Int_t module=cluster->GetModule();
231 
232  Int_t range_set;
233  if (((module==1)||(module==2))&&(cluster_energy<=1.0))
234  {
235  range_set=1;
236  }
237  else if (((module==1)||(module==2))&&(cluster_energy>1.0))
238  {
239  range_set=2;
240  }
241  else
242  {
243  range_set=module;
244  }
245 
246  Int_t thetaBin=GetThetaBin(cluster_theta, range_set); // bin number for cluster_theta
247  Int_t energyBin=GetEnergyBin(cluster_energy, range_set);// bin number for cluster_energy
248 
249  if (cluster_theta >= 141. && cluster_theta < 147.) continue; // Avoid the egges between barrel and bwendcap...
250 
251  Double_t en_div, phi_diff;
252  if ((thetaBin<0&&range_set!=5) || energyBin<0)
253  {
254  cout<<"!!!!!!!!!!!!!!!!!!!!!!"<<endl;
255  cout<<"Module="<<module<<" range_set="<<range_set<<endl;
256  cout << "theta_reco =" << cluster_theta << "/ energy_reco =" << cluster_energy << endl;
257  cout << "thetaBin= "<<thetaBin<<" & energyBin= "<<energyBin<<endl;
258  }
259  else
260  {
261  if (cluster_theta > (thMC-5) && cluster_theta <=(thMC+5) ){
262 
263  en_div = enMC/cluster_energy;
264  phi_diff = phiMC-cluster_phi;
265 
266  if (fabs(phi_diff)<2.5 || fabs(phi_diff+360)<2.5 || fabs(phi_diff-360)<2.5)
267  {
268  switch (range_set)
269  {
270  case 1:
271  {
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  }
276  break;
277  }
278  case 2:
279  {
280  if ((fabs(phiMC) <5.)||fabs(phiMC-180.)<5.||fabs(phiMC-360.)<5.) break;
281  if (en_div > 0.7 && en_div <1.3){
282  hisEnergyRatio2[energyBin][thetaBin].Fill(en_div);
283  }
284  break;
285  }
286  case 3:
287  {
288  if (cluster_theta <11.) break; // Avoid the egges for FwEndCap...
289  if (en_div > 0.7 && en_div <1.3){
290  hisEnergyRatio3[energyBin][thetaBin].Fill(en_div);
291  }
292  break;
293  }
294  case 4:
295  {
296  if (en_div > 0.7 && en_div <1.3){
297  hisEnergyRatio4[energyBin][thetaBin].Fill(en_div);
298  }
299  break;
300  }
301  case 5:
302  {
303  if ((en_div > 0.7 && en_div <1.3)&&(thMC>1.)){ // avoid edge around 0 degree
304  hisEnergyRatio5[energyBin].Fill(en_div);
305  }
306  break;
307  }
308  default:
309  {
310  abort();
311  }
312  }
313  } // Phi
314  } // Theta
315  } // ThetaBin & EnergyBin > 0
316  }
317  }
318 
319  std::cout<<"Array of histogram is filled"<<std::endl;
320 timer.Stop();
321 Double_t rtime = timer.RealTime();
322 Double_t ctime = timer.CpuTime();
323 timer.Continue();
324 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
325  // Loop over energy and theta bins
326  Double_t mean_en_ratio=1.;
327 
328  for (Int_t i=0; i<(nrEnergyIntervals[0]-1); i++)
329  {
330  for (Int_t j=0; j<(nrThetaIntervals[0]-1); j++)
331  {
332  if ( hisEnergyRatio1[i][j].GetEntries() == 0) {
333  mean_en_ratio = 1;
334  }else{
335  mean_en_ratio = hisEnergyRatio1[i][j].GetMean();
336  }
337  hisEnergyRatioBarrelLow->SetBinContent(i+1,j+1,mean_en_ratio);
338  }
339  }
340 
341  for (Int_t i=0; i<(nrEnergyIntervals[1]-1); i++)
342  {
343  for (Int_t j=0; j<(nrThetaIntervals[1]-1); j++)
344  {
345  if ( hisEnergyRatio2[i][j].GetEntries() == 0) {
346  mean_en_ratio = 1;
347  }else{
348  mean_en_ratio = hisEnergyRatio2[i][j].GetMean();
349  }
350  hisEnergyRatioBarrelHigh->SetBinContent(i+1,j+1,mean_en_ratio);
351  }
352  }
353 
354  for (Int_t i=0; i<(nrEnergyIntervals[2]-1); i++)
355  {
356  for (Int_t j=0; j<(nrThetaIntervals[2]-1); j++)
357  {
358  if ( hisEnergyRatio3[i][j].GetEntries() == 0) {
359  mean_en_ratio = 1;
360  }else{
361  mean_en_ratio = hisEnergyRatio3[i][j].GetMean();
362  }
363  hisEnergyRatioFwd->SetBinContent(i+1,j+1,mean_en_ratio);
364  }
365  }
366 
367  for (Int_t i=0; i<(nrEnergyIntervals[3]-1); i++)
368  {
369  for (Int_t j=0; j<(nrThetaIntervals[3]-1); j++)
370  {
371  if ( hisEnergyRatio4[i][j].GetEntries() == 0) {
372  mean_en_ratio = 1;
373  }else{
374  mean_en_ratio = hisEnergyRatio4[i][j].GetMean();
375  }
376  hisEnergyRatioBwd->SetBinContent(i+1,j+1,mean_en_ratio);
377  }
378  }
379 
380  for (Int_t i=0; i<(nrEnergyIntervals[4]-1); i++)
381  {
382  if ( hisEnergyRatio5[i].GetEntries() == 0) {
383  mean_en_ratio = 1;
384  }else{
385  mean_en_ratio = hisEnergyRatio5[i].GetMean();
386  }
387  hisEnergyRatioShashlyk->SetBinContent(i+1,mean_en_ratio);
388  }
389 
390  if (storeHistos)
391  {
392  f_histos=new TFile("correction_histos_fit.root","RECREATE");
393  hisEnergyRatioBarrelLow->Write();
394  hisEnergyRatioBarrelHigh->Write();
395  hisEnergyRatioFwd->Write();
396  hisEnergyRatioBwd->Write();
397  hisEnergyRatioShashlyk->Write();
398  f_histos->Close();
399  }
400  std::cout<<"Histograms are filled with mean values"<<std::endl;
401  }
402 
403  if (useStoredHistos&&(!storeHistos))
404  {
405  f_histos=new TFile("correction_histos_fit.root");
406  hisEnergyRatioBarrelLow=(TH2F*)f_histos->Get("hisEnergyRatioBarrelLow");
407  hisEnergyRatioBarrelHigh=(TH2F*)f_histos->Get("hisEnergyRatioBarrelHigh");
408  hisEnergyRatioFwd=(TH2F*)f_histos->Get("hisEnergyRatioFwd");
409  hisEnergyRatioBwd=(TH2F*)f_histos->Get("hisEnergyRatioBwd");
410  hisEnergyRatioShashlyk=(TH1F*)f_histos->Get("hisEnergyRatioShashlyk");
411  }
412 
413  TF2 * func1 = new TF2("func1",FitFunction1,0.03,1.0,25.,135., 10);
414  double iniParams1[10] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
415  func1->SetParameters(iniParams1);
416  TCanvas *c1=new TCanvas("c1");
417  hisEnergyRatioBarrelLow->Fit(func1,"RN"); //"SURF1"
418  func1->Draw("SURF1");
419 
420  TF2 * func2 = new TF2("func2",FitFunction1,1.0,8.0,25.,135., 10);
421  double iniParams2[10] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
422  func2->SetParameters(iniParams2);
423  TCanvas *c2=new TCanvas("c2");
424  hisEnergyRatioBarrelHigh->Fit(func2,"RN");
425  func2->Draw("SURF1");//"cont1"
426 
427  TF2 * func3 = new TF2("func3",FitFunction1,0.03,8.0,11.,22., 10);
428  double iniParams3[10] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
429  func3->SetParameters(iniParams3);
430  TCanvas *c3=new TCanvas("c3");
431  hisEnergyRatioFwd->Fit(func3,"RN");
432  func3->Draw("SURF1");
433 
434  TF2 * func4 = new TF2("func4",FitFunction1,0.03,2.0,150.,163., 10);
435  double iniParams4[10] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
436  func4->SetParameters(iniParams4);
437  TCanvas *c4=new TCanvas("c4");
438  hisEnergyRatioBwd->Fit(func4,"RN");
439  func4->Draw("SURF1");
440 
441  TF1 * func5 = new TF1("func5",FitFunction2,0.0,10.0, 5);
442  double iniParams5[10] = { 0.01, 0.01, 0.01, 0.01, 0.01};
443  func5->SetParameters(iniParams5);
444  TCanvas *c5=new TCanvas("c5");
445  hisEnergyRatioShashlyk->Fit(func5,"R");
446 
447 
448  if (debug)
449  {
450  c1->SaveAs("c1.root");
451  c2->SaveAs("c2.root");
452  c3->SaveAs("c3.root");
453  c4->SaveAs("c4.root");
454  c5->SaveAs("c5.root");
455  }
456  TFile fout(OutputFile,"recreate"); // output file with corrections
457 
458  PndEmcClusterCalibrationParObject *parObject=new PndEmcClusterCalibrationParObject();
459  enum {barrel_low=1, barrel_high=2, fwcap=3, bwcap=4, fsc=5};
460  Double_t pars1[10], pars2[10], pars3[10], pars4[10], pars5[5];
461 
462  Double_t *p1=func1->GetParameters();
463  for (int i=0;i<10;i++)
464  pars1[i]=p1[i];
465  Double_t *p2=func2->GetParameters();
466  for (int i=0;i<10;i++)
467  pars2[i]=p2[i];
468  Double_t *p3=func3->GetParameters();
469  for (int i=0;i<10;i++)
470  pars3[i]=p3[i];
471  Double_t *p4=func4->GetParameters();
472  for (int i=0;i<10;i++)
473  pars4[i]=p4[i];
474  Double_t *p5=func5->GetParameters();
475  for (int i=0;i<5;i++)
476  pars5[i]=p5[i];
477 
478  parObject->SetCalibrationPar(barrel_low, pars1);
479  parObject->SetCalibrationPar(barrel_high, pars2);
480  parObject->SetCalibrationPar(fwcap, pars3);
481  parObject->SetCalibrationPar(bwcap, pars4);
482  parObject->SetCalibrationPar(fsc, pars5);
483 
484  parObject->Write();
485 
486  fout.Close();
487 timer.Stop();
488 Double_t rtime = timer.RealTime();
489 Double_t ctime = timer.CpuTime();
490 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
491 
492  return 0;
493 }
energyIntervalsFwd
Double_t energyIntervalsFwd[]
Definition: emc_correction_parametrization.C:11
thetaIntervalsBwd
Double_t thetaIntervalsBwd[]
Definition: emc_correction_hist.C:15
phi_diff
double phi_diff
Definition: macro/detectors/emc/dedicated/reco_analys.C:25
c5
c5
Definition: plot_dirc.C:75
FitFunction1
Double_t FitFunction1(Double_t *x, Double_t *par)
Definition: emc_correction_parametrization.C:497
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
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
PndEmcCluster::GetModule
Short_t GetModule() const
Definition: PndEmcCluster.cxx:282
f1
TF1 * f1
Definition: reco_analys2.C:50
energyIntervalsBarrelLow
Double_t energyIntervalsBarrelLow[]
Definition: emc_correction_parametrization.C:3
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
c2
c2
Definition: plot_dirc.C:39
FitFunction2
Double_t FitFunction2(Double_t *x, Double_t *par)
Definition: emc_correction_parametrization.C:519
PndMCTrack::Get4Momentum
TLorentzVector Get4Momentum() const
Definition: PndMCTrack.cxx:102
particle
const int particle
Definition: runGeaneTestSimBox.cxx:51
thetaIntervalsBarrelLow
Double_t thetaIntervalsBarrelLow[]
Definition: emc_correction_parametrization.C:4
cluster_theta
Double_t cluster_theta
Definition: anal_hit_digi_cluster_fwendcap.C:31
cluster_array
TClonesArray * cluster_array
Definition: anal_hit_digi_cluster_fwendcap.C:27
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
energyIntervalsBarrelHigh
Double_t energyIntervalsBarrelHigh[]
Definition: emc_correction_parametrization.C:7
timer
TStopwatch timer
Definition: hit_dirc.C:51
energyIntervalsBwd
Double_t energyIntervalsBwd[]
Definition: emc_correction_parametrization.C:15
track
PndMCTrack * track
Definition: anaLmdCluster.C:89
cluster_energy
Double_t cluster_energy
Definition: anal_hit_digi_cluster_fwendcap.C:31
nrEnergyIntervals
Int_t nrEnergyIntervals
Definition: emc_correction_hist.C:18
cluster_phi
Double_t cluster_phi
Definition: anal_hit_digi_cluster_fwendcap.C:31
energyIntervalsShashlyk
Double_t energyIntervalsShashlyk[]
Definition: emc_correction_parametrization.C:19
nrThetaIntervals
Int_t nrThetaIntervals[4]
Definition: emc_correction_hist.C:19
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
PndEmcCluster
a cluster (group of neighboring crystals) of hit emc crystals
Definition: PndEmcCluster.h:29
c1
c1
Definition: plot_dirc.C:35
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
c3
c3
Definition: plot_dirc.C:50
ctime
Double_t ctime
Definition: hit_dirc.C:114
p2
TPad * p2
Definition: hist-t7.C:117
PndEmcClusterCalibrationParObject::SetCalibrationPar
void SetCalibrationPar(Int_t iParSet, Double_t *pars)
Definition: PndEmcClusterCalibrationPar.cxx:38
f2
TFile * f2
Definition: ReadMesDecCorrRes.C:5
GetEnergyBin
Int_t GetEnergyBin(Double_t val, Int_t module)
Definition: emc_correction_hist.C:503
GetEntries
cout<<"will loop over "<< t-> GetEntries()
Definition: root2ascii.C:17
PndEmcCluster::energy
virtual Double_t energy() const
Definition: PndEmcCluster.cxx:79
p1
TPad * p1
Definition: hist-t7.C:116
thetaIntervalsBarrelHigh
Double_t thetaIntervalsBarrelHigh[]
Definition: emc_correction_parametrization.C:8
track_array
TClonesArray * track_array
Definition: plot_rk.C:34
p4mom
TLorentzVector p4mom
Definition: anal_hit_digi_cluster_fwendcap.C:46
rtime
Double_t rtime
Definition: hit_dirc.C:113
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
infile
TFile infile("dedx_out.root","READ")
thetaIntervalsFwd
Double_t thetaIntervalsFwd[]
Definition: emc_correction_hist.C:14
file_name
string file_name
Definition: structure_analysis_emc.C:23
GetThetaBin
Int_t GetThetaBin(Double_t val, Int_t module)
Definition: emc_correction_hist.C:466
Double_t FitFunction1 ( Double_t *  x,
Double_t *  par 
)

Definition at line 497 of file emc_correction_parametrization.C.

References cos(), exp(), log(), res, and theta.

Referenced by emc_correction_parametrization().

498 {
499  double e=x[0];
500  double theta=x[1]*TMath::DegToRad();
501 
502  double factor1= par[0]
503  +par[1]*log(e)
504  +par[2]*log(e)*log(e)
505  +par[3]*log(e)*log(e)*log(e)
506  +par[4]*cos(theta)
507  +par[5]*cos(theta)*cos(theta)
508  +par[6]*cos(theta)*cos(theta)*cos(theta)
509  +par[7]*cos(theta)*cos(theta)*cos(theta)*cos(theta)
510  +par[8]*cos(theta)*cos(theta)*cos(theta)*cos(theta)*cos(theta)
511  +par[9]*log(e)*cos(theta);
512 
513  double res=exp(factor1);
514 
515  return res;
516 }
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: P4_F32vec4.h:109
res
Int_t res
Definition: anadigi.C:166
par
Double_t par[3]
Definition: anaLmdCluster.C:146
log
friend F32vec4 log(const F32vec4 &a)
Definition: P4_F32vec4.h:110
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Double_t FitFunction2 ( Double_t *  x,
Double_t *  par 
)

Definition at line 519 of file emc_correction_parametrization.C.

References res, and sqrt().

Referenced by emc_correction_parametrization().

520 {
521  double e=x[0];
522  double res=(par[0]-par[1]/sqrt(e)+par[2]/e+par[3]/(e*e)-par[4]/(e*sqrt(e)));
523  return res;
524 }
res
Int_t res
Definition: anadigi.C:166
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
par
Double_t par[3]
Definition: anaLmdCluster.C:146
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Int_t GetEnergyBin ( Double_t  val,
Int_t  range_set 
)

Definition at line 565 of file emc_correction_parametrization.C.

References Double_t, energyIntervals, energyIntervalsBarrelHigh, energyIntervalsBarrelLow, energyIntervalsBwd, energyIntervalsFwd, energyIntervalsShashlyk, i, and nrEnergyIntervals.

566 {
567  Double_t *energyIntervals=0;
568  Int_t nrEnergyInt=0;
569  switch (range_set)
570  {
571  case 1:
572  nrEnergyInt=nrEnergyIntervals[0];
573  energyIntervals=energyIntervalsBarrelLow;
574  break;
575  case 2:
576  nrEnergyInt=nrEnergyIntervals[1];
577  energyIntervals=energyIntervalsBarrelHigh;
578  break;
579  case 3:
580  nrEnergyInt=nrEnergyIntervals[2];
581  energyIntervals=energyIntervalsFwd;
582  break;
583  case 4:
584  nrEnergyInt=nrEnergyIntervals[3];
585  energyIntervals=energyIntervalsBwd;
586  break;
587  case 5:
588  nrEnergyInt=nrEnergyIntervals[4];
589  energyIntervals=energyIntervalsShashlyk;
590  break;
591  default:
592  std::cout<<"Wrong EMC module: "<<range_set<<std::endl;
593  abort();
594  }
595 
596  for (Int_t i=0; i<(nrEnergyInt-1); i++)
597  {
598  if (val>=energyIntervals[i] && val<energyIntervals[i+1])
599  {
600  return (i);
601  }
602  }
603  return -1;
604 
605 }
energyIntervalsFwd
Double_t energyIntervalsFwd[]
Definition: emc_correction_parametrization.C:11
i
Int_t i
Definition: run_full.C:25
val
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
energyIntervalsBarrelLow
Double_t energyIntervalsBarrelLow[]
Definition: emc_correction_parametrization.C:3
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
energyIntervalsBarrelHigh
Double_t energyIntervalsBarrelHigh[]
Definition: emc_correction_parametrization.C:7
energyIntervalsBwd
Double_t energyIntervalsBwd[]
Definition: emc_correction_parametrization.C:15
nrEnergyIntervals
Int_t nrEnergyIntervals
Definition: emc_correction_hist.C:18
energyIntervalsShashlyk
Double_t energyIntervalsShashlyk[]
Definition: emc_correction_parametrization.C:19
energyIntervals
Double_t energyIntervals[]
Definition: emc_correction_hist.C:12
Int_t GetThetaBin ( Double_t  val,
Int_t  range_set 
)

Definition at line 526 of file emc_correction_parametrization.C.

References Double_t, i, nrThetaIntervals, thetaIntervalsBarrelHigh, thetaIntervalsBarrelLow, thetaIntervalsBwd, and thetaIntervalsFwd.

527 {
528  Double_t *thetaIntervals=0;
529  Int_t nrThetaInt=0;
530  switch (range_set)
531  {
532  case 1:
533  nrThetaInt=nrThetaIntervals[0];
534  thetaIntervals=thetaIntervalsBarrelLow;
535  break;
536  case 2:
537  nrThetaInt=nrThetaIntervals[1];
538  thetaIntervals=thetaIntervalsBarrelHigh;
539  break;
540  case 3:
541  nrThetaInt=nrThetaIntervals[2];
542  thetaIntervals=thetaIntervalsFwd;
543  break;
544  case 4:
545  nrThetaInt=nrThetaIntervals[3];
546  thetaIntervals=thetaIntervalsBwd;
547  break;
548  case 5:
549  return -1;
550  default:
551  std::cout<<"Wrong EMC module: "<<range_set<<std::endl;
552  abort();
553  }
554 
555  for (Int_t i=0; i<(nrThetaInt-1); i++)
556  {
557  if (val>=thetaIntervals[i] && val<thetaIntervals[i+1])
558  {
559  return (i);
560  }
561  }
562  return -1;
563 }
thetaIntervalsBwd
Double_t thetaIntervalsBwd[]
Definition: emc_correction_hist.C:15
i
Int_t i
Definition: run_full.C:25
val
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
thetaIntervalsBarrelLow
Double_t thetaIntervalsBarrelLow[]
Definition: emc_correction_parametrization.C:4
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
nrThetaIntervals
Int_t nrThetaIntervals[4]
Definition: emc_correction_hist.C:19
thetaIntervalsBarrelHigh
Double_t thetaIntervalsBarrelHigh[]
Definition: emc_correction_parametrization.C:8
thetaIntervalsFwd
Double_t thetaIntervalsFwd[]
Definition: emc_correction_hist.C:14

Variable Documentation

Double_t energyIntervalsBarrelHigh[] = {1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}

Definition at line 7 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetEnergyBin().

Double_t energyIntervalsBarrelLow[] = {0,0.03, 0.05, 0.07, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}

Definition at line 3 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetEnergyBin().

Double_t energyIntervalsBwd[] = {0,0.03,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0}

Definition at line 15 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetEnergyBin().

Double_t energyIntervalsFwd[] = {0,0.03, 0.05, 0.07, 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}

Definition at line 11 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetEnergyBin().

Double_t energyIntervalsShashlyk[] = {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,9.0,10.0}

Definition at line 19 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetEnergyBin().

Int_t nrEnergyIntervals[5]

Definition at line 21 of file emc_correction_parametrization.C.

Int_t nrThetaIntervals[4]

Definition at line 22 of file emc_correction_parametrization.C.

Double_t thetaIntervalsBarrelHigh[] = {22.,23.,24.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100.,105.,110.,115.,120.,125.,130.,135.,137., 139., 140., 141.}

Definition at line 8 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetThetaBin().

Double_t thetaIntervalsBarrelLow[] = {22.,23.,24.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,95.,100.,105.,110.,115.,120.,125.,130.,135.,137., 139., 140., 141.}

Definition at line 4 of file emc_correction_parametrization.C.

Referenced by emc_correction_parametrization(), and GetThetaBin().

Double_t thetaIntervalsBwd[] = {147.,148.,149.,150.,151.,152.,153.,154.,155.,156.,157.,158.,159.,160.,161.,162.,163.,164.,166.,168.,170.,172.}

Definition at line 16 of file emc_correction_parametrization.C.

Double_t thetaIntervalsFwd[] = {5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,20.,21.,22.}

Definition at line 12 of file emc_correction_parametrization.C.

Generated on Wed Apr 3 2019 10:03:13 for FairRoot/PandaRoot by   doxygen 1.8.5