FairRoot/PandaRoot
PndLmdPerformanceTask.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndLmdPerformance Task source file -----
3 // ----- Based on PndLmdGeaneTask -----
4 // ----- Author: Prometeusz Jasinski -----
5 // -------------------------------------------------------------------------
6 // libc includes
7 #include <iostream>
8 
9 // Root includes
10 #include "TROOT.h"
11 #include "TClonesArray.h"
12 #include "TVector3.h"
13 #include <TMatrixDSym.h>
14 #include <TCanvas.h>
15 #include <TColor.h>
16 #include <TStyle.h>
17 #include <sstream>
18 
19 // framework includes
20 #include "FairRootManager.h"
21 #include "PndLmdPerformanceTask.h"
22 #include "FairRun.h"
23 #include "FairRuntimeDb.h"
24 #include "FairHit.h"
25 #include "PndMCTrack.h"
26 #include "FairBaseParSet.h"
27 #include "TGeant3.h"
28 #include "FairTrackParH.h"
29 #include "FairTrackParP.h"
30 #include "TDatabasePDG.h"
31 
32 // PndSds includes
33 #include "PndSdsMCPoint.h"
34 #include "PndSdsHit.h"
35 //PndLmd includes
36 #include "PndLinTrack.h"
37 #include "TGaxis.h"
38 //#include "PndLmdDim.h"
39 
40 // GenFit includes
41 #include "RKTrackRep.h"
42 #include "GFFieldManager.h"
43 #include "PndGenfitField.h"
44 //#include "GFPandaField.h"
45 
46 #include <vector>
47 #include <map>
48 
49 // ----- Default constructor -------------------------------------------
50 PndLmdPerformanceTask::PndLmdPerformanceTask() :
51  FairTask("Geane Task for PANDA Lmd"), fEventNr(0), fUseMVDPoint(false) {
52  fPro = NULL;
53  fGeoH = NULL;
54  hist_output_file = NULL;
55  last_percent = -1;
56  fnevents = -1;
57  fievent = 0;
58 }
59 // -------------------------------------------------------------------------
60 
61 PndLmdPerformanceTask::PndLmdPerformanceTask(Double_t pBeam, TVector3 IP) :
62  FairTask("Geane Task for PANDA Lmd"), fEventNr(0), fUseMVDPoint(false) {
63  fPbeam = pBeam;
64  // if(fPDGid!=-2212){ //calculate momentum if particle is not antiproton
65  // }
66  hist_output_file = NULL;
67  cout << "Beam Momentum for particle with PDGid#" << fPDGid
68  << " this run is " << fPbeam << endl;
69  vtx = IP;
70  cout << "Interaction Point:" << endl;
71  vtx.Print();
72  last_percent = -1;
73  fnevents = -1;
74  fievent = 0;
75 }
76 
77 // ----- Destructor ----------------------------------------------------
78 PndLmdPerformanceTask::~PndLmdPerformanceTask() {
79  if (hist_output_file) hist_output_file->Close();
80  delete hist_output_file;
81 }
82 
83 // ----- Public method Init --------------------------------------------
84 InitStatus PndLmdPerformanceTask::Init() {
85 
86  //if (hist_output_file)
87  // hist_output_file->cd();
88 
89  // Get RootManager
90  FairRootManager* ioman = FairRootManager::Instance();
91  if (!ioman) {
92  std::cout << "-E- PndLmdPerformanceTask::Init: "
93  << "RootManager not instantiated!" << std::endl;
94  return kFATAL;
95  } else {
96  fnevents = ioman->CheckMaxEventNo();
97  }
98 
99  /*
100  fMCHits = (TClonesArray*) ioman->GetObject("LMDPoint");
101  if (!fMCHits) {
102  std::cout << "-W- PndLmdPerformanceTask::Init: " << "No LMDPoint"
103  << " array!" << std::endl;
104  return kERROR;
105  }*/
106 
107  /*fTracks = (TClonesArray*) ioman->GetObject("LMDTrack");
108  if (!fTracks) {
109  std::cout << "-W- PndLmdPerformanceTask::Init: " << "No Track"
110  << " array!" << std::endl;
111  return kERROR;
112  }
113 
114  fHits = (TClonesArray*) ioman->GetObject("LMDHitsStrip");
115  if (!fHits) {
116  std::cout << "-W- PndLmdPerformanceTask::Init: " << "No Hits"
117  << " array!" << std::endl;
118  return kERROR;
119  }*/
120 
121  //--- MC info -----------------------------------------------------------------
122  true_tracks = (TClonesArray*) ioman->GetObject("MCTrack"); //True Track to compare
123  /* if (!true_tracks) {
124  std::cout << "-W- PndLmdPerformanceTask::Init: " << "No MCTrack"
125  << " array!" << std::endl;
126  return kFATAL;
127  }*/
128  true_points = (TClonesArray*) ioman->GetObject("LMDPoint"); //True Points to compare
129  //----------------------------------------------------------------------------------
130 
131 
132  //--- Digitization info ------------------------------------------------------------
133  //fStripClusterArray
134  // = (TClonesArray*) ioman->GetObject("LMDStripClusterCand");
135  //fStripDigiArray = (TClonesArray*) ioman->GetObject("LMDStripDigis");
136  //----------------------------------------------------------------------------------
137 
138  //--- Real Hits --------------------------------------------------------------------
139  //rechit_array = (TClonesArray*) ioman->GetObject("LMDHitsStrip"); //Points for Tracks
140  //----------------------------------------------------------------------------------
141 
142 
143  //--- Track Candidate ---------------------------------------------------------------
144  //trkcand_array = (TClonesArray*) ioman->GetObject("LMDTrackCand"); //Points for Track Canidates
145  //-----------------------------------------------------------------------------------
146 
147  //--- Real tracks -------------------------------------------------------------------
148  //rec_trk = (TClonesArray*) ioman->GetObject("LMDTrack"); //Tracks
149  //----------------------------------------------------------------------------------
150 
151  fDetName = new TClonesArray("TObjString");
152  ioman->Register("DetName", "Geane", fDetName, kTRUE);
153 
154  fPro = new FairGeanePro();
155  //fPro->PropagateToPCA(1, -1);
156  fGeoH = PndGeoHandling::Instance();
157  //FairRun* fRun = FairRun::Instance();
158  //FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
159  lmddim = PndLmdDim::Instance();
160  lmddim -> Read_transformation_matrices("",true);
161  lmddim -> Read_transformation_matrices("",false);
162  //lmddim -> Read_transformation_matrices("matrices.txt", true);
163  //lmddim -> Read_transformation_matrices("matrices_perfect.txt", false);
164  // FairBaseParSet* par=(FairBaseParSet*)
165  // (rtdb->findContainer("FairBaseParSet"));
166  // fPbeam = par->GetBeamMom();
167  // cout<<"Beam Momentum for this run is "<<fPbeam<<endl;
168 
169  GFFieldManager::getInstance()->init(new PndGenfitField());
170 
171  std::cout << " Setting up histograms in PndLmdPerformanceTask ";
172  flush(std::cout);
173 
174  TPad foo; // never remove this line :-)))
175  if(1){
176  gROOT->SetStyle("Plain");
177  const Int_t NRGBs = 5;
178  const Int_t NCont = 255;
179  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
180  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
181  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
182  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
183  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
184  gStyle->SetNumberContours(NCont);
185  gStyle->SetTitleFont(10*13+2,"xyz");
186  gStyle->SetTitleSize(0.06, "xyz");
187  gStyle->SetTitleOffset(1.3,"y");
188  gStyle->SetTitleOffset(1.3,"z");
189  gStyle->SetLabelFont(10*13+2,"xyz");
190  gStyle->SetLabelSize(0.06,"xyz");
191  gStyle->SetLabelOffset(0.009,"xyz");
192  gStyle->SetPadBottomMargin(0.16);
193  gStyle->SetPadTopMargin(0.16);
194  gStyle->SetPadLeftMargin(0.16);
195  gStyle->SetPadRightMargin(0.16);
196  gStyle->SetOptTitle(1);
197  gStyle->SetOptStat(1);
198  gROOT->ForceStyle();
199  gStyle->SetFrameFillColor(0);
200  gStyle->SetFrameFillStyle(0);
201  TGaxis::SetMaxDigits(3);
202  }
203 
204  if (!hist_output_file) SetHistFilename("_hists.root");
205  hist_output_file->cd();
206 
207  hist_angular_distr_gen = new TH2F("hist_angular_distr_gen",
208  "hist_angular_distr_gen", 400, 2.e-3, 12.e-3, 400, -3.141, 3.141);
209  hist_angular_distr_acc = new TH2F("hist_angular_distr_acc",
210  "hist_angular_distr_acc", 400, 2.e-3, 12.e-3, 400, -3.141, 3.141);
211 
212  // spatial acceptance in x and y at the first lumi plane
213  hist_spatial_distr_gen = lmddim->Get_histogram_Plane(0, 0, true, true, false);// new TH2F("hist_spatial_distr_gen",
214  //"hist_spatial_distr_gen", 1000, -100., 100., 1000, -100., 100.);
215  hist_spatial_distr_acc = lmddim->Get_histogram_Plane(1, 0, true, true, false);// new TH2F("hist_spatial_distr_acc",
216  //"hist_spatial_distr_acc", 1000, -100., 100., 1000, -100., 100.);
217 
218  hist_theta_over_mom_gen = new TH2F("hist_theta_over_mom_gen",
219  "hist_theta_over_mom_gen", 2000, 1., 20., 100, 2.e-3, 12.e-3);
220  hist_theta_over_mom_acc = new TH2F("hist_theta_over_mom_acc",
221  "hist_theta_over_mom_acc", 2000, 1., 20., 100, 2.e-3, 12.e-3);
222 
223  hist_phi_over_mom_gen = new TH2F("hist_phi_over_mom_gen",
224  "hist_phi_over_mom_gen", 2000, 1., 20., 400, -3.141, 3.141);
225  hist_phi_over_mom_acc = new TH2F("hist_phi_over_mom_acc",
226  "hist_phi_over_mom_acc", 2000, 1., 20., 400, -3.141, 3.141);
227 
228  tree_results = new TTree("tree_results", "tree_results");
229  // initial momenta of the antiprotons
230  tree_results->Branch("px_init", &px_init);
231  tree_results->Branch("py_init", &py_init);
232  tree_results->Branch("pz_init", &pz_init);
233  tree_results->Branch("ptheta_init", &ptheta_init);
234  tree_results->Branch("pphi_init", &pphi_init);
235  // momenta at the entry of a detector plane
236  // in the reference system of the LUMI
237  tree_results->Branch("px_in", &px_in);
238  tree_results->Branch("py_in", &py_in);
239  tree_results->Branch("pz_in", &pz_in);
240  tree_results->Branch("ptheta_in", &ptheta_in);
241  tree_results->Branch("pphi_in", &pphi_in);
242  // position at the entry of a detector plane
243  // in the reference system of the LUMI
244  tree_results->Branch("x_in", &x_in);
245  tree_results->Branch("y_in", &y_in);
246  tree_results->Branch("z_in", &z_in);
247  // detector and sensor id of the mc hits
248  tree_results->Branch("det_id", &det_id);
249  tree_results->Branch("sens_id", &sens_id);
250  // calculated plane and sensor in plane id
251  tree_results->Branch("plane", &plane);
252  tree_results->Branch("sensor", &sensor);
253  tree_results->Branch("missed", &missed);
254  // new HVMAPS variables
255  tree_results->Branch("ihalf", &ihalf);
256  tree_results->Branch("iplane", &iplane);
257  tree_results->Branch("imodule", &imodule);
258  tree_results->Branch("iside", &iside);
259  tree_results->Branch("idie", &idie);
260  tree_results->Branch("isensor", &isensor);
261  tree_results->Branch("x_in_mod", &x_in_mod); // reference frame on the module surface
262  tree_results->Branch("y_in_mod", &y_in_mod);
263  tree_results->Branch("z_in_mod", &z_in_mod);
264  tree_results->Branch("x_in_sens", &x_in_sens); // reference frame on the sensor
265  tree_results->Branch("y_in_sens", &y_in_sens);
266  tree_results->Branch("z_in_sens", &z_in_sens);
267  tree_results->Branch("x_in_sens_al", &x_in_sens_al); // reference frame on the sensor
268  tree_results->Branch("y_in_sens_al", &y_in_sens_al);
269  tree_results->Branch("z_in_sens_al", &z_in_sens_al);
270  tree_results->Branch("x_in_aligned", &x_in_aligned); // reference frame on the module surface
271  tree_results->Branch("y_in_aligned", &y_in_aligned); // aligned coordinates in the lumi frame
272  tree_results->Branch("z_in_aligned", &z_in_aligned);
273  tree_results->Branch("theta_prop", &theta_prop);
274  tree_results->Branch("theta_prop_aligned", &theta_prop_aligned);
275 
276  // Drawing directly from a tree is elegant but slow as every Draw call
277  // loops over the whole tree of events
278  // Better it is to create the desired histograms in advance
279  // and to fill those on event by event basis simultaneously
280 
281  // histograms per plane
282  //hist_xy[nplanes];
283  //hist_theta_init[nplanes];
284  //hist_theta_in[nplanes];
285  //hist_theta_diff[nplanes];
286  //hist_theta_diff_rel[nplanes];
287 
288  // histograms per sensor
289  //hists_xy[nplanes][nsensors_per_plane];
290  //hists_theta_init[nplanes][nsensors_per_plane];
291  //hists_theta_in[nplanes][nsensors_per_plane];
292  //hists_theta_diff[nplanes][nsensors_per_plane];
293  //hists_theta_diff_rel[nplanes][nsensors_per_plane];
294 
295  TCanvas temp_canvas("temp_canvas", "canvas for initialization", 100, 100);
296  temp_canvas.cd();
297  // loop over planes
298  cout << " constructing histograms for " << nplanes << " planes with "
299  << nsensors_per_plane << " sensors per plane " << endl;
300 
301  stringstream hist_name;
302  stringstream hist_title;
303  for (unsigned int _iplane = 0; _iplane < nplanes; _iplane++) {
304 
305  hist_name .str("");
306  hist_title.str("");
307 
308  hist_name << "hist_xy_plane_" << _iplane;
309  hist_title << "xy hit distribution plane " << _iplane;
310  hist_xy[_iplane] = new TH2F(hist_name.str().c_str(),
311  hist_title.str().c_str(), 100, -10, 10, 100, -10, 10);
312  hist_xy[_iplane]->Draw();
313  hist_xy[_iplane]->GetXaxis()->SetTitle("X [cm]");
314  hist_xy[_iplane]->GetYaxis()->SetTitle("Y [cm]");
315 
316  hist_name .str("");
317  hist_title.str("");
318 
319  hist_name << "hist_theta_init_plane_" << _iplane;
320  hist_title << "initial #Theta distribution plane " << _iplane;
321  hist_theta_init[_iplane] = new TH2F(hist_name.str().c_str(),
322  hist_title.str().c_str(), 150, 0, 15e-3, 360, -3.141, +3.141);
323  hist_theta_init[_iplane]->Draw();
324  hist_theta_init[_iplane]->GetXaxis()->SetTitle("#Theta [rad]");
325  hist_theta_init[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
326 
327  hist_name .str("");
328  hist_title.str("");
329 
330  hist_name << "hist_theta_in_plane_" << _iplane;
331  hist_title << "#Theta distribution at plane " << _iplane;
332  hist_theta_in[_iplane] = new TH2F(hist_name.str().c_str(),
333  hist_title.str().c_str(), 200, 0, 20e-3, 360, -3.141, +3.141);
334  hist_theta_in[_iplane]->Draw();
335  hist_theta_in[_iplane]->GetXaxis()->SetTitle("#Theta [rad]");
336  hist_theta_in[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
337 
338  hist_name .str("");
339  hist_title.str("");
340 
341  hist_name << "hist_rec_theta_with_plane_" << _iplane;
342  hist_title << "#Theta reconstructed distribution with plane " << _iplane;
343  hist_theta_rec[_iplane] = new TH2F(hist_name.str().c_str(),
344  hist_title.str().c_str(), 200, 0, 20e-3, 360, -3.141, +3.141);
345  hist_theta_rec[_iplane]->Draw();
346  hist_theta_rec[_iplane]->GetXaxis()->SetTitle("#Theta [rad]");
347  hist_theta_rec[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
348 
349  hist_name .str("");
350  hist_title.str("");
351 
352  hist_name << "hist_theta_diff_in_plane_" << _iplane;
353  hist_title << "#Delta#Theta distribution at plane " << _iplane;
354  hist_theta_diff[_iplane] = new TH2F(hist_name.str().c_str(),
355  hist_title.str().c_str(), 100, -40e-4, 40e-4, 360, -3.141, +3.141);
356  hist_theta_diff[_iplane]->Draw();
357  hist_theta_diff[_iplane]->GetXaxis()->SetTitle("#Delta#Theta [rad]");
358  hist_theta_diff[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
359 
360  hist_name .str("");
361  hist_title.str("");
362 
363  hist_name << "hist_theta_rec_diff_in_plane_" << _iplane;
364  hist_title << "#Delta#Theta reco distribution at plane " << _iplane;
365  hist_theta_rec_diff[_iplane] = new TH2F(hist_name.str().c_str(),
366  hist_title.str().c_str(), 100, -40e-4, 40e-4, 360, -3.141, +3.141);
367  hist_theta_rec_diff[_iplane]->Draw();
368  hist_theta_rec_diff[_iplane]->GetXaxis()->SetTitle("#Delta#Theta [rad]");
369  hist_theta_rec_diff[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
370 
371  hist_name .str("");
372  hist_title.str("");
373 
374  hist_name << "hist_theta_diff_rel_in_plane_" << _iplane;
375  hist_title << "#Delta#Theta/#Theta distribution at plane " << _iplane;
376  hist_theta_diff_rel[_iplane] = new TH2F(hist_name.str().c_str(),
377  hist_title.str().c_str(), 100, -1, 1, 360, -3.141, +3.141);
378  hist_theta_diff_rel[_iplane]->Draw();
379  hist_theta_diff_rel[_iplane]->GetXaxis()->SetTitle("#Delta#Theta/#Theta");
380  hist_theta_diff_rel[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
381 
382  hist_name .str("");
383  hist_title.str("");
384 
385  hist_name << "hist_theta_rec_diff_rel_in_plane_" << _iplane;
386  hist_title << "#Delta#Theta/#Theta reco distribution at plane " << _iplane;
387  hist_theta_rec_diff_rel[_iplane] = new TH2F(hist_name.str().c_str(),
388  hist_title.str().c_str(), 100, -1, 1, 360, -3.141, +3.141);
389  hist_theta_rec_diff_rel[_iplane]->Draw();
390  hist_theta_rec_diff_rel[_iplane]->GetXaxis()->SetTitle("#Delta#Theta/#Theta");
391  hist_theta_rec_diff_rel[_iplane]->GetYaxis()->SetTitle("#Phi [rad]");
392 
393  // loop over sensors per plane which are enumerated linearly
394  for (unsigned int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
395  hist_name .str("");
396  hist_title.str("");
397 
398  hist_name << "hist_xy_plane_" << _iplane << "_sensor_" << _isensor;
399  hist_title << "xy hit distribution plane " << _iplane << " sensor "
400  << _isensor;
401  hists_xy[_iplane][_isensor] = new TH2F(hist_name.str().c_str(),
402  hist_title.str().c_str(), 100, -10, 10, 100, -10, 10);
403  hists_xy[_iplane][_isensor]->Draw();
404  hists_xy[_iplane][_isensor]->GetXaxis()->SetTitle("X [cm]");
405  hists_xy[_iplane][_isensor]->GetYaxis()->SetTitle("Y [cm]");
406 
407  hist_name .str("");
408  hist_title.str("");
409 
410  hist_name << "hist_xy_local_plane_" << _iplane << "_sensor_" << _isensor;
411  hist_title << "xy local hit distribution plane " << _iplane << " sensor "
412  << _isensor;
413  hists_xy_local[_iplane][_isensor] = new TH2F(hist_name.str().c_str(),
414  hist_title.str().c_str(), 100, -2, 2, 100, -2, 2);
415  hists_xy_local[_iplane][_isensor]->Draw();
416  hists_xy_local[_iplane][_isensor]->GetXaxis()->SetTitle("X [cm]");
417  hists_xy_local[_iplane][_isensor]->GetYaxis()->SetTitle("Y [cm]");
418 
419  hist_name << "hist_theta_init_plane_" << _iplane << "_sensor_"
420  << _isensor;
421  hist_title << "initial #Theta distribution plane " << _iplane
422  << " sensor " << _isensor;
423  hists_theta_init[_iplane][_isensor] = new TH1F(
424  hist_name.str().c_str(), hist_title.str().c_str(), 150, 0,
425  15e-3);
426  hists_theta_init[_iplane][_isensor]->Draw();
427  hists_theta_init[_iplane][_isensor]->GetXaxis()->SetTitle(
428  "#Theta [rad]");
429  hists_theta_init[_iplane][_isensor]->GetYaxis()->SetTitle("entries");
430 
431  hist_name .str("");
432  hist_title.str("");
433 
434  hist_name << "hist_theta_in_plane_" << _iplane << "_sensor_"
435  << _isensor;
436  hist_title << "#Theta distribution at plane " << _iplane
437  << " sensor " << _isensor;
438  hists_theta_in[_iplane][_isensor] = new TH1F(hist_name.str().c_str(),
439  hist_title.str().c_str(), 200, 0, 20e-3);
440  hists_theta_in[_iplane][_isensor]->Draw();
441  hists_theta_in[_iplane][_isensor]->GetXaxis()->SetTitle(
442  "#Theta [rad]");
443  hists_theta_in[_iplane][_isensor]->GetYaxis()->SetTitle("entries");
444 
445  hist_name .str("");
446  hist_title.str("");
447 
448  hist_name << "hist_theta_diff_in_plane_" << _iplane << "_sensor_"
449  << _isensor;
450  hist_title << "#Delta#Theta distribution at plane " << _iplane
451  << " sensor " << _isensor;
452  hists_theta_diff[_iplane][_isensor] = new TH1F(
453  hist_name.str().c_str(), hist_title.str().c_str(), 100,
454  -20e-5, 20e-5);
455  hists_theta_diff[_iplane][_isensor]->Draw();
456  hists_theta_diff[_iplane][_isensor]->GetXaxis()->SetTitle(
457  "#Delta#Theta [rad]");
458  hists_theta_diff[_iplane][_isensor]->GetYaxis()->SetTitle("entries");
459 
460  hist_name .str("");
461  hist_title.str("");
462 
463  hist_name << "hist_theta_diff_rel_in_plane_" << _iplane
464  << "_sensor_" << _isensor;
465  hist_title << "#Delta#Theta/#Theta distribution at plane "
466  << _iplane << " sensor " << _isensor;
467  hists_theta_diff_rel[_iplane][_isensor] = new TH1F(
468  hist_name.str().c_str(), hist_title.str().c_str(), 100, -1,
469  1);
470  hists_theta_diff_rel[_iplane][_isensor]->Draw();
471  hists_theta_diff_rel[_iplane][_isensor]->GetXaxis()->SetTitle(
472  "#Delta#Theta/#Theta");
473  hists_theta_diff_rel[_iplane][_isensor]->GetYaxis()->SetTitle(
474  "entries");
475  }
476  }
477 
478  hist_name .str("");
479  hist_title.str("");
480 
481  hist_name << "hist_theta_diff_prop_true";
482  hist_title << "#Delta#Theta true distribution (Error by Geane back propagation) ";
483  hist_theta_diff_prop_true = new TH2F(hist_name.str().c_str(),
484  hist_title.str().c_str(), 300, -10e-5, 10e-5, 360, -3.141, +3.141);
485  hist_theta_diff_prop_true->Draw();
486  hist_theta_diff_prop_true->GetXaxis()->SetTitle("#Delta#Theta [rad]");
487  hist_theta_diff_prop_true->GetYaxis()->SetTitle("#Phi [rad]");
488 
489  hist_name .str("");
490  hist_title.str("");
491 
492  hist_name << "hist_theta_diff_prop_true_o_theta";
493  hist_title << "#Delta#Theta true distribution (Error by Geane back propagation) ";
494  hist_theta_diff_prop_true_o_theta = new TH2F(hist_name.str().c_str(),
495  hist_title.str().c_str(), 300, -10e-5, 10e-5, 200, 2.e-3, 10.e-3);
496  hist_theta_diff_prop_true_o_theta->Draw();
497  hist_theta_diff_prop_true_o_theta->GetXaxis()->SetTitle("#Delta#Theta [rad]");
498  hist_theta_diff_prop_true_o_theta->GetYaxis()->SetTitle("#Theta [rad]");
499 
500 
501  std::cout << " done " << std::endl;
502 
503  // rotation into the lumi position according to the beamline file
504  // matrix tri6 - tr=1 rot=1 refl=0 scl=0
505  // 0.999192 0.000000 0.040200 Tx = 22.893063
506  // 0.000000 1.000000 -0.000000 Ty = 0.000000
507  // -0.040200 0.000000 0.999192 Tz = 1049.770390
508  // The matrix is a rotation around y of around 0.04025 radian
509  // TRotation inv_lmdrotation;
510  inv_lmdrotation.RotateY(-0.04025);
511  inv_lmdrotation.RotateY(-40.068e-3);
512  inv_lmdtranslation.SetXYZ(-26.2461, 0.000000, -1130.);
513 
514  fgGeoMan = (TGeoManager*) gROOT->FindObject("FAIRGeom");
515  if (!fgGeoMan) cout << "Error: could not find the geometry manager!" << endl;
516 
517  return kSUCCESS;
518 }
519 // -------------------------------------------------------------------------
520 void PndLmdPerformanceTask::SetParContainers() {
521  // Get Base Container
523  //s FairRuntimeDb* rtdb=ana->GetRuntimeDb();
524 
525 }
526 
527 void PndLmdPerformanceTask::SetHistFilename(TString filename){
528  if (hist_output_file){
529  cout << " Warning in PndLmdPerformanceTask::SetHistFilename(): output file was set to " << hist_output_file->GetName() << endl;
530  cout << " Setting it to: " << filename << endl;
531  hist_output_file->Close();
532  }
533  hist_output_file = new TFile(filename, "RECREATE");
534 }
535 
536 /*
537 void PndLmdPerformanceTask::Propagate(TVector3 pos, TVector3 mom, int pdgcode,
538  TVector3 ip, TVector3& initpos, TVector3& initdir) {
539  //TVector3 initpos, TVector3 initdir, TVector3 IP, TVector3& respos, TVector3& resdir, ", &);mom){ // PDGCode
540  if (!fPro) {
541  cout << " no propagator found, creating one and initilizing the rest "
542  << endl;
543  Init();
544  //fPro = new FairGeanePro();
545  //fGeoH = PndGeoHandling::Instance();
546  //FairRun* fRun = FairRun::Instance();
547  //FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
548  }
549  fPbeam = initdir.Mag();
550  // if(fPDGid!=-2212){ //calculate momentum if particle is not antiproton
551  // }
552  //cout<<"Beam Momentum for particle with PDGid#"<<fPDGid<<" this run is "<<fPbeam<<endl;
553  vtx = ip;
554  //cout<<"Interaction Point:"<<endl;
555  //vtx.Print();
556  TVector3 StartPosErr, StartMomErr, ResPosErr, ResMomErr;
557  int charge = pdgcode > 0 ? 1 : pdgcode < 0 ? -1 : 0;
558  FairTrackParH *fStart = new FairTrackParH(pos, mom, StartPosErr,
559  StartMomErr, charge); // exchange fCharge by charge of PDGCode
560  //TClonesArray& clref = *fTrackParGeane;
561  //Int_t size = clref.GetEntriesFast();
562  FairTrackParH *fRes = new FairTrackParH(initpos, initdir, ResPosErr,
563  ResMomErr, charge);
564  fPro->SetPoint(ip);
565  // fPro->BackTrackToVertex();
566  Bool_t isProp = fPro->Propagate(fStart, fRes, pdgcode);
567  delete fStart;
568  delete fRes;
569 
570  if (!isProp)
571  cout
572  << "PndLmdPerformanceTask::Propagate(): back propagation not possible! "
573  << endl;
574 
575  //TVector3 gPos(fRes->GetX(),fRes->GetY(),fRes->GetZ());
576  //TVector3 gMom(fRes->GetPx(),fRes->GetPy(),fRes->GetPz());
577  //TVector3 gErrPos(fRes->GetDX(),fRes->GetDY(),fRes->GetDZ());
578  //TVector3 gErrMom(fRes->GetDPx(),fRes->GetDPy(),fRes->GetDPz());
579 
580  //initpos = fRes->GetPosition();
581  //initdir = fRes->GetMomentum();
582 }*/
583 
584 // ----- Public method Exec --------------------------------------------
585 void PndLmdPerformanceTask::Exec(Option_t*) {
586 
587  bool verbose(false);
588 
589  fTrackPlaneHitMap.clear();
590  fTrackMCVectMapIn.clear();
591  fTrackMCVectMapInit.clear();
592 
593  fievent++;
594  if (fnevents == -1) fnevents = 1000000;
595 
596  DrawProgressBar(50, (fievent) / ((double) fnevents));
597 
598  const int nParticles = true_tracks->GetEntriesFast();
599 
600  int npbar = 0;
601  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
602  // if(nParticles!=4) continue;
603  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(iParticle);
604  //Int_t mcID = mctrk->GetPdgCode();
605  if (mctrk->IsGeneratorCreated()) {//mcID == -2212) {
606  npbar++;
607  // if(nParticles!=4) cout<<"For event #"<<j<<" mcID="<<mcID<<endl;
608  TVector3 MomMC = mctrk->GetMomentum();
609  hist_angular_distr_gen->Fill(MomMC.Theta(), MomMC.Phi());
610  hist_theta_over_mom_gen->Fill(MomMC.Mag(), MomMC.Theta());
611  hist_phi_over_mom_gen->Fill(MomMC.Mag(), MomMC.Phi());
612  //hallTheta->Fill(MomMC.Theta());
613  //hallPhi->Fill(MomMC.Phi());
614  fTrackMCVectMapInit.push_back(MomMC);
615  }
616  }
617  if (verbose && npbar != 1)
618  cout << "found " << npbar << " anti-protons" << endl;
619 
620  int nHits = true_points->GetEntriesFast();
621  int nSdsHits = 0;
622  //cout << " hits " << nHits << endl;
623  for (Int_t iHit = 0; iHit < nHits; iHit++) {
624  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
625  PndMCTrack *mctrk =
626  (PndMCTrack*) true_tracks->At(mcpoint->GetTrackID());
627  if (mctrk->IsGeneratorCreated()) { // take only anti-protons from the generator
628  if (mcpoint) {
629  TVector3 _mcpoint((mcpoint->GetPosition())); // position at the sensor entrance
630  //cout << " momentum " << mctrk->GetMomentum().Mag();
631  //cout << " momentum difference " << (mctrk->GetMomentum().Mag() - fTrackMCVectMapInit[0].Mag())/fTrackMCVectMapInit[0].Mag();
632  //cout << " Energy loss " << mcpoint->GetEnergyLoss() << endl;
633  TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
634  mcpoint->GetPz()); // momentum of the track at the entrance
635  // variables used for conversion via root geometries
636  //double glob_pt[3] = {mcpoint->GetX(), mcpoint->GetY(), mcpoint->GetZ() + 1.e-4};
637  //double loc_pt[3] = {0,0,0};
638  if (0) { // perform some geane tests
639  if (mcpoint->GetSensorID() > 7) return;
640  TVector3 StartPos, StartPosErr, StartMom, StartMomErr, StartO, StartU, StartV;
641  //TVector3 DirVec = _momrec.Unit();
642  //cout << " here " << ihere++ << endl;
643  StartMom = _mctrack;//.Unit()*1.5*(1.-1.e-3);//_momrec.Unit()*itmominit.Mag();
644  TVector3 move_upstream_vect(0,0,-0.1);
645  StartPos = _mcpoint + move_upstream_vect ;
646  StartPosErr = TVector3(50.e-4, 50.e-4, 50.e-4);// TVector3(50.e-4, 50.e-4, 50.e-4); // 50 mu like pitch size
647  StartMomErr = _mctrack.Mag()*1e-3;//_momrec.Unit()*fPbeam*1e-3; //_momrec.Unit()*fPbeam*1e-3; // some estimate
648  //cout << " start mom " << StartMom.X() << '\t' << StartMom.Y() << '\t' << StartMom.Z() << '\t' << endl;
649  //cout << " start mom err " << StartMomErr.X() << '\t' << StartMomErr.Y() << '\t' << StartMomErr.Z() << '\t' << endl;
650  //cout << " start pos " << StartPos.X() << '\t' << StartPos.Y() << '\t' << StartPos.Z() << '\t' << endl;
651  //cout << " start pos err " << StartPosErr.X() << '\t' << StartPosErr.Y() << '\t' << StartPosErr.Z() << '\t' << endl;
652  int PDGCode = -2212;
653  TVector3 _vtx(0.,0.,0.); // point to find the closest approach to
654  TVector3 gPos(0.,0.,0.); // point of closest approach
655  TVector3 gMom(0.,0.,0.); // momentum at that point
656  if (0) {// geane methods implemented by turany
657  //if (fievent > 100000) return; // debugging
658  //if (_mctrack.Phi() < 0.09 || _mctrack.Phi() > 0.1) return;
659  //cout << " here " << ihere++ << endl;
660  FairTrackParH *fStart = new FairTrackParH(StartPos, StartMom, StartPosErr, StartMomErr, -1);
661  FairTrackParH *fRes = new FairTrackParH();
662  //cout << " here " << ihere++ << endl;
663  fPro->SetPoint(_vtx);
664  fPro->PropagateToPCA(1,-1);
665  // fPro->BackTrackToVertex();
666  Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
667  if (!isProp){
668  cout << " Warning in PndLmdPerformanceTask::Exec(): propagation failed! " << endl;
669  }
670  gPos.SetXYZ(fRes->GetX(),fRes->GetY(),fRes->GetZ());
671  gMom.SetXYZ(fRes->GetPx(),fRes->GetPy(),fRes->GetPz());
672  //_momrec = gMom;
673  delete fStart;
674  delete fRes;
675  } else { // genfit Runge Kutta Nystroem algorithm
676  RKTrackRep* trackrep = new RKTrackRep(StartPos, StartMom, StartPosErr, StartMomErr, PDGCode);
677  trackrep->extrapolateToPoint(_vtx, gPos, gMom);
678  }
679  if (fTrackMCVectMapInit.size() > 0){
680  double theta_lmd_true = gMom.Theta();
681  double theta_true = fTrackMCVectMapInit[0].Theta(); // reference momentum direction
682  TVector3 _mcpoint_rot = inv_lmdtranslation + _mcpoint;
683  _mcpoint_rot = inv_lmdrotation * _mcpoint_rot;
684  //if (-2.4 < _mcpoint_rot.Phi() && _mcpoint_rot.Phi() < -2.3){
685  hist_theta_diff_prop_true->Fill(theta_lmd_true - theta_true, fTrackMCVectMapInit[0].Phi());// _mcpoint_rot.Phi());// fTrackMCVectMapInit[0].Phi());
686  hist_theta_diff_prop_true_o_theta->Fill(theta_lmd_true - theta_true, theta_true);
687  //}
688  if (fTrackMCVectMapInit.size() > 1 && verbose){
689  cout << " Warning in PndLmdPerformanceTask::Exec(): more than 1 IP MC anti proton found -> using first found" << endl;
690  }
691  }
692  return;
693  }
694  if (1) {
695  //const double* current_point = fgGeoMan->GetCurrentPoint();
696  //cout << current_point[0] << " " << current_point[1] << " " << current_point[2] << endl;
697  TGeoNode* node =
698  fgGeoMan->FindNode(mcpoint->GetX(), mcpoint->GetY(), mcpoint->GetZ() + 1.e-4);
699  if (node){
700  missed = false;
701  //std::cout << node->GetName() << std::endl;
702  //std::cout << node->GetIndex() << std::endl;
703  //node->GetMatrix()->Print();
704  //std::cout << << std::endl;
705  //std::cout << node->GetMotherVolume()->GetName() << std::endl;
706  //std::cout << node->GetNumber() << "\n" << std::endl;
707  TGeoVolume* local_top = fgGeoMan->GetVolume("lmd_vol_ref_sys");
708  if (! local_top ){
709  std::cout << " Error: Could not set top volume to lmd_vol_ref_sys" << std::endl;
710  } else {
711  //fgGeoMan->SetTopVolume(local_top);
712  //cout << fgGeoMan->GetPath() << endl;
713  //fgGeoMan->MasterToLocal(glob_pt, loc_pt);
714 
715  }
716  //std::cout << glob_pt[0] << " " << glob_pt[1] << " " << glob_pt[2] << " " << std::endl;
717  //std::cout << loc_pt[0] << " " << loc_pt[1] << " " << loc_pt[2] << " " << std::endl;
718  } else {
719  missed = true;
720  std::cout << " Error: no node found for the current MC point " << std::endl;
721  //mcpoint->Print();
722  }
723  /*
724  double new_point[3] = {8.7, 1., 29.7375};
725  fgGeoMan->SetCurrentPoint(new_point);
726  const double* new_current_point = fgGeoMan->GetCurrentPoint();
727  cout << new_current_point[0] << " " << new_current_point[1] << " " << new_current_point[2] << endl;
728  cout << fgGeoMan->GetCurrentNode()->GetName() << endl;
729  gGeoMan->FindNode();
730  cout << fgGeoMan->GetCurrentNode()->GetName() << endl;
731 
732  string path = "/lmd_HV_MAPS_1/vol_lmd_vac_0/Lumi_HV-MAPS_0/vol_LumPassive_cvd_disc_plane_1_disc_1_0";
733  fgGeoMan->cd(path.c_str());///LumActivePixelRect_plane_0_disc_0_side_0_col_0_row_0_0");
734  cout << fgGeoMan->GetPath()<< endl;
735  TGeoPhysicalNode* node =
736  fgGeoMan->MakePhysicalNode((path + "/LumActivePixelRect_plane_0_disc_0_side_0_col_0_row_0_0").c_str());///Lumi_HV-MAPS/LumPassive_cvd_disc_plane_2_disc_3");
737  if (node){
738  //TGeoCombiTrans* align_trans = new TGeoCombiTrans(1.,0.,0.,georot_no);
739  //node->Align(align_trans);
740  //gGeoMan->CloseGeometry();
741  TGeoHMatrix* matrix = node -> GetMatrix();
742  matrix -> Print();
743  }*/
744  }
745 
746  if (0) { // translate it into the reference system of the lumi monitor
747  //_mcpoint = inv_lmdtranslation + _mcpoint;
748  //_mcpoint = inv_lmdrotation * _mcpoint;
749  }
750  if (1) {
751  _mcpoint = lmddim->Transform_global_to_lmd_local(_mcpoint, false, false);
752  }
753  TVector3 momMC = mctrk->GetMomentum(); // momentum in the primary vertex
754  TVector3 posMC = mctrk->GetStartVertex(); // position of the primary vertex
755  if (1) { // use back propagation for the track
756  //TVector3 initpos, initmom(0, 0, momMC.Mag()), StartPosErr,
757  // StartMomErr, ResPosErr, ResMomErr;
758  //int charge = mctrk->GetPdgCode() > 0 ? 1
759  // : mctrk->GetPdgCode() < 0 ? -1 : 0;
760  //FairTrackParH *fStart = new FairTrackParH(posMC, momMC,
761  // StartPosErr, StartMomErr, charge); // exchange fCharge by charge of PDGCode
762  //FairTrackParH *fRes = new FairTrackParH(initpos, initmom,
763  // ResPosErr, ResMomErr, charge);
764  //FairTrackParH *fRes = new FairTrackParH();
765  //fPro->SetPoint(vtx);
766  //fPro->PropagateToPCA(1, -1);
767  //Bool_t isProp = fPro->Propagate(fStart, fRes,
768  // mctrk->GetPdgCode());
769  //if (!isProp){
770  //cout << "Propagate(): back propagation not possible! "
771  // << endl;
772  //}
773  //else {
774  //cout << "Propagation worked! " << endl;
775  //}
776  //_mctrack = fRes->GetMomentum();//
777  //_mctrack = inv_lmdrotation*_mctrack;
778  //lmddim->Transform_global_to_lmd_local(_mctrack, true);
779  //delete fStart;
780  //delete fRes;
781  }
782  det_id = mcpoint->GetDetectorID(); // store the detector id
783  sens_id = mcpoint->GetSensorID(); // store the sensor id
784  // calculate the plane and sensor on this plane
785  lmddim->Get_sensor_by_id(sens_id, ihalf, iplane, imodule, iside, idie, isensor);
786  sensor = isensor;
787  plane = iplane;
788  if (idie == 1) sensor += 2;
789  sensor += lmddim->n_sensors * iside;
790  sensor += lmddim->n_sensors * 2 * imodule;
791  sensor += lmddim->n_sensors * 2 * lmddim->nmodules * ihalf; // I know, I know these lines are not optimized ;)
792  //plane = sens_id / nsensors_per_plane;
793  //sensor = sens_id - plane * nsensors_per_plane;
794  TVector3 _mctrack_aligned = lmddim->Transform_global_to_sensor(_mctrack, ihalf, iplane, imodule, iside, idie, isensor, true, true);
795  _mctrack = lmddim->Transform_global_to_sensor(_mctrack, ihalf, iplane, imodule, iside, idie, isensor, true, false);
796 
797  // write those values out
798  px_init = momMC.X();
799  py_init = momMC.Y();
800  pz_init = momMC.Z();
801  ptheta_init = momMC.Theta();
802  pphi_init = momMC.Phi();
803  x_in = _mcpoint.X();
804  y_in = _mcpoint.Y();
805  z_in = _mcpoint.Z();
806  px_in = _mctrack.X();
807  py_in = _mctrack.Y();
808  pz_in = _mctrack.Z();
809  ptheta_in = _mctrack.Theta();
810  pphi_in = _mctrack.Phi();
811 
812  const TVector3& _mcpoint_mod =
813  lmddim->Transform_lmd_local_to_module_side(_mcpoint,
814  ihalf, iplane, imodule, iside, false, false);
815  x_in_mod = _mcpoint_mod.X();
816  y_in_mod = _mcpoint_mod.Y();
817  z_in_mod = _mcpoint_mod.Z();
818  const TVector3& _mcpoint_sens =
819  lmddim->Transform_global_to_sensor(mcpoint->GetPosition(),
820  ihalf, iplane, imodule, iside, idie, isensor, false, false);
821  x_in_sens = _mcpoint_sens.X(); // reference frame on the sensor
822  y_in_sens = _mcpoint_sens.Y();
823  z_in_sens = _mcpoint_sens.Z();
824  const TVector3& _mcpoint_sens_al =
825  lmddim->Transform_global_to_sensor(mcpoint->GetPosition(),
826  ihalf, iplane, imodule, iside, idie, isensor, false, true);
827  x_in_sens_al = _mcpoint_sens_al.X(); // reference frame on the sensor
828  y_in_sens_al = _mcpoint_sens_al.Y();
829  z_in_sens_al = _mcpoint_sens_al.Z();
830  const TVector3& _mcpoint_aligned =
831  lmddim->Transform_global_to_lmd_local(mcpoint->GetPosition(), false, true);
832  x_in_aligned = _mcpoint_aligned.X(); // reference frame on the module surface
833  y_in_aligned = _mcpoint_aligned.Y(); // aligned coordinates in the lumi frame
834  z_in_aligned = _mcpoint_aligned.Z();
835 
836  theta_prop = _mctrack.Theta();
837  theta_prop_aligned = _mctrack_aligned.Theta();
838 
839  if (plane < 0 || plane >= nplanes || sensor < 0 || sensor
840  >= nsensors_per_plane) {
841  cout
842  << " warning: calculated plane or sensor is wrong: plane "
843  << plane << " sensor " << sensor << " sensorID " << sens_id << endl;
844  } else {
845  // we need the unrotated vectors, too
846  TVector3 mcpoint_org((mcpoint->GetPosition())); // position at the sensor entrance
847  TVector3 mctrack_org(mcpoint->GetPx(), mcpoint->GetPy(),
848  mcpoint->GetPz()); // momentum of the track at the entrance
849  fTrackPlaneHitMap[plane].push_back(mcpoint_org);
850  fTrackMCVectMapIn[plane].push_back(mctrack_org);
851 
852  hist_xy[plane] ->Fill(x_in, y_in);
853  hists_xy[plane][sensor]->Fill(x_in, y_in);
854  hists_xy_local[plane][sensor]->Fill(_mcpoint_sens.X(), _mcpoint_sens.Y());
855 
856  hist_theta_init[plane] ->Fill(ptheta_init, pphi_init);
857  hists_theta_init[plane][sensor]->Fill(ptheta_init);
858 
859  hist_theta_in[plane] ->Fill(ptheta_in, pphi_in);
860  hists_theta_in[plane][sensor]->Fill(ptheta_in);
861 
862  hist_theta_diff[plane] ->Fill(ptheta_in - ptheta_init, pphi_init);
863  hists_theta_diff[plane][sensor]->Fill(
864  ptheta_in - ptheta_init, pphi_init);
865 
866  if (ptheta_init > 0) {
867  hist_theta_diff_rel[plane] ->Fill(
868  (ptheta_in - ptheta_init) / ptheta_init, pphi_init);
869  hists_theta_diff_rel[plane][sensor]->Fill(
870  (ptheta_in - ptheta_init) / ptheta_init, pphi_init);
871  }
872 
873  tree_results->Fill();
874  }
875 
876  //cout << _mcpoint.X() << " " << _mcpoint.Y() << " "<< _mcpoint.Z() << endl;
877  hist_spatial_distr_acc->Fill(_mcpoint.X(), _mcpoint.Y());
878  //hist_angular_distr_acc->Fill(momMC.Theta(), momMC.Phi()); // moved to reconstruction point
879  //hist_theta_over_mom_acc->Fill(momMC.Mag(), momMC.Theta());
880  //hist_phi_over_mom_acc->Fill(momMC.Mag(), momMC.Phi());
881  nSdsHits++;
882  }
883  }
884  }
885  if (verbose && nHits > 0)
886  cout << "found " << nHits << " hit(s) with " << nSdsHits
887  << " sds hit(s)" << endl;
888 
889  // perform a simple tracking using hits between 2 planes without digitization
890  // to see only effects of smearing
891  std::map<int, std::vector<TVector3> >::iterator itplane1, itplane2, itmomin;
892  // find a hit in the first plane (simplified tracking here needs the first plane)
893  itplane1 = fTrackPlaneHitMap.find(0);
894  // as well as the momentum in the first plane
895  itmomin = fTrackMCVectMapIn.find(0);
896  bool accepted(false);
897  if (itplane1 != fTrackPlaneHitMap.end() && itmomin != fTrackMCVectMapIn.end()){
898  if (fTrackMCVectMapInit.size() > 0){
899  TVector3 itmominit = fTrackMCVectMapInit[0]; // reference momentum direction
900  if (fTrackMCVectMapInit.size() > 1 && verbose){
901  cout << " Warning in PndLmdPerformanceTask::Exec(): more than 1 IP MC anti proton found -> using first found" << endl;
902  }
903 
904  if ((*itmomin).second.size() > 0){
905  TVector3 _itmomin = (*itmomin).second[0]; // reference momentum direction at plane 0
906  if ((*itmomin).second.size() > 1 && verbose){
907  cout << " Warning in PndLmdPerformanceTask::Exec(): more than 1 IP MC anti proton momentum in the first plane found -> using first found" << endl;
908  }
909  // check the first found point
910  if ((*itplane1).second.size() > 0){
911  TVector3 point1 = (*itplane1).second[0];
912  if ((*itplane1).second.size() > 1 && verbose){
913  cout << " Warning in PndLmdPerformanceTask::Exec(): more than 1 IP MC anti proton hit in the plane 0 found -> using first found" << endl;
914  }
915  // find partners in other planes
916  for (int _iplane = 1; _iplane < nplanes; _iplane++){
917  itplane2 = fTrackPlaneHitMap.find(_iplane);
918  if (itplane2 != fTrackPlaneHitMap.end()){
919  // reconstruct the vector between the first plane and that one
920  if ((*itplane2).second.size() > 0){
921  TVector3 point2 = (*itplane2).second[0];
922  if ((*itplane2).second.size() > 1 && verbose){
923  cout << " Warning in PndLmdPerformanceTask::Exec(): more than 1 IP MC anti proton hit in the plane "<< _iplane <<" found -> using first found" << endl;
924  }
925  TVector3 _momrec = point2 - point1;
926  //int ihere(0);
927  //cout << " here " << ihere++ << endl;
928  // back propagation
929  if (0){
930  //cout << " here " << ihere++ << endl;
931  TVector3 StartPos, StartPosErr, StartMom, StartMomErr, StartO, StartU, StartV;
932  //TVector3 DirVec = _momrec.Unit();
933  //cout << " here " << ihere++ << endl;
934  StartMom = _momrec.Unit()*itmominit.Mag();
935  StartPos = point1;
936  StartPosErr = TVector3(50.e-4, 50.e-4, 50.e-4);// TVector3(50.e-4, 50.e-4, 50.e-4); // 50 mu like pitch size
937  StartMomErr = _momrec.Unit()*fPbeam*1e-3; //_momrec.Unit()*fPbeam*1e-3; // some estimate
938  TVector3 _vtx(0.,0.,0.); //
939  int PDGCode = -2212;
940  //cout << " start mom " << StartMom.X() << '\t' << StartMom.Y() << '\t' << StartMom.Z() << '\t' << endl;
941  //cout << " start mom err " << StartMomErr.X() << '\t' << StartMomErr.Y() << '\t' << StartMomErr.Z() << '\t' << endl;
942  //cout << " start pos " << StartPos.X() << '\t' << StartPos.Y() << '\t' << StartPos.Z() << '\t' << endl;
943  //cout << " start pos err " << StartPosErr.X() << '\t' << StartPosErr.Y() << '\t' << StartPosErr.Z() << '\t' << endl;
944  FairTrackParH *fStart = new FairTrackParH(StartPos, StartMom, StartPosErr, StartMomErr, -1);
945  FairTrackParH *fRes = new FairTrackParH();
946  //cout << " here " << ihere++ << endl;
947  TVector3 gPos(0.,0.,0.);
948  TVector3 gMom(0.,0.,0.);
949  if (0){ // use geane back propagation else ...
950  fPro->SetPoint(_vtx);
951  fPro->PropagateToPCA(1,-1);
952  // fPro->BackTrackToVertex();
953  Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
954  if (!isProp){
955  cout << " Warning in PndLmdPerformanceTask::Exec(): propagation failed! " << endl;
956  }
957  gPos.SetXYZ(fRes->GetX(),fRes->GetY(),fRes->GetZ());
958  gMom.SetXYZ(fRes->GetPx(),fRes->GetPy(),fRes->GetPz());
959 
960  delete fStart;
961  delete fRes;
962  } else { // ... use Runge-Kutta algorithm
963 
964  RKTrackRep* trackrep = new RKTrackRep(StartPos, StartMom, StartPosErr, StartMomErr, PDGCode);
965  trackrep->extrapolateToPoint(_vtx, gPos, gMom);
966 
967  }
968 
969  _momrec = gMom;
970  }
971  hist_theta_rec[_iplane]->Fill(_momrec.Theta(), _momrec.Phi());
972  hist_theta_rec_diff[_iplane]->Fill(_momrec.Theta() - itmominit.Theta(), itmominit.Phi());
973  hist_theta_rec_diff_rel[_iplane]->Fill((_momrec.Theta() - itmominit.Theta())/itmominit.Theta(), itmominit.Phi());
974  if (!accepted){ // fill the following histograms only once per event
975  hist_angular_distr_acc->Fill(itmominit.Theta(), itmominit.Phi());
976  hist_theta_over_mom_acc->Fill(itmominit.Mag(), itmominit.Theta());
977  hist_phi_over_mom_acc->Fill(itmominit.Mag(), itmominit.Phi());
978  accepted = true;
979  }
980  }
981  }
982  }
983  }
984  }
985  }
986  }
987 
988  /*
989 
990  // cout<<"PndLmdPerformanceTask::Exec starts!"<<endl;
991  std::map<int, std::vector<int> > mcHitMap;//Track -> MCHits
992  fTrackParGeane->Delete();
993  fTrackParIni->Delete();
994  fTrackParFinal->Delete();
995  fDetName->Delete();
996 
997  mcHitMap = AssignHitsToTracks();
998  fEventNr++;
999 
1000  //Charge & mass of particle
1001  Int_t PDGCode = -2212; //antiproton
1002  // Int_t PDGCode = fPDGid;
1003  TDatabasePDG *fdbPDG = TDatabasePDG::Instance();
1004  TParticlePDG *fParticle = fdbPDG->GetParticle(PDGCode);
1005  Double_t fCharge = fParticle->Charge();
1006 
1007  //go through all tracks
1008  int glI = fTracks->GetEntriesFast();
1009  // cout<<"glI = "<<glI<<endl;
1010  Int_t counterGeaneTrk = 0;
1011  for (Int_t i = 0; i < glI; i++) {
1012  // cout<<"PndLmdPerformanceTask::Exec for track#"<<i<<endl;
1013  TVector3 StartPos, StartPosErr, StartMom, StartMomErr, StartO, StartU,
1014  StartV;
1015  int p = 0;
1016 
1018  PndLinTrack* recTrack = (PndLinTrack*) (fTracks->At(i));
1019  // cout<<"recTrack="<<recTrack<<endl;
1020  TString DecName = recTrack->GetDetName();
1021  if (DecName != "Lumi")
1022  continue;
1023  // cout<<"DecName: "<<DecName<<endl;
1024  //Vector of particle momentum and starting point
1025  TVector3 DirVec = recTrack->GetDirectionVec();
1026  StartMom = TVector3(DirVec.X() * fPbeam, DirVec.Y() * fPbeam,
1027  DirVec.Z() * fPbeam);
1028  StartPos = recTrack->GetStartVec();
1029  StartPosErr = recTrack->GetStartErrVec();
1030  StartMomErr = (recTrack->GetDirectionErrVec()) * fPbeam;
1031  if (fVerbose > 2) {
1032  cout << "------------------------------------------" << endl;
1033  cout << "StartPos:" << endl;
1034  StartPos.Print();
1035  cout << "StartPosErr:" << endl;
1036  StartPosErr.Print();
1037  cout << "" << endl;
1038  cout << "StartMom: " << StartMom.Mag() << endl;
1039  StartMom.Print();
1040  cout << "StartMomErr: " << StartMomErr.Mag() << endl;
1041  StartMomErr.Print();
1042  }
1043  TClonesArray& clref1 = *fTrackParIni;
1044  Int_t size1 = clref1.GetEntriesFast();
1045 
1047  FairTrackParH *fStart = new (clref1[size1]) FairTrackParH(StartPos,
1048  StartMom, StartPosErr, StartMomErr, fCharge);
1049  TClonesArray& clref = *fTrackParGeane;
1050  Int_t size = clref.GetEntriesFast();
1051  FairTrackParH *fRes = new (clref[size]) FairTrackParH();
1052  fPro->SetPoint(vtx);
1053  fPro->PropagateToPCA(1, -1);
1054  // fPro->BackTrackToVertex();
1055  Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
1057 
1058  // ///Forwars propagate to the 1st plane to a space point---------------------------------
1059  // PndMCTrack* mctrk = (PndMCTrack*)(fMCTracks->At(i));
1060  // TVector3 MomMC = mctrk->GetMomentum();
1061  // TVector3 PosMC = mctrk->GetStartVertex();
1062  // FairTrackParH *fStart =
1063  // new (clref1[size1]) FairTrackParH(PosMC, MomMC, StartPosErr, StartMomErr, fCharge);
1064  // TClonesArray& clref = *fTrackParGeane;
1065  // Int_t size = clref.GetEntriesFast();
1066  // FairTrackParH *fRes = new(clref[size]) FairTrackParH();
1067  // fPro->SetPoint(StartPos);
1068  // //fPro->PropagateToPCA(1,-1);
1069  // fPro->PropagateToPCA(1,1);
1070  // Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
1071  // ///----------------------------------------------------------------------
1072 
1073 
1074  // ///Propagate to virtual plane at PCA ------------------------------------
1075  // TVector3 oc = (0,0,0);
1076  // TVector3 dj(0,1,0);
1077  // TVector3 dk(-1,0,0);
1078  // FairTrackParP *fStart =
1079  // new (clref1[size1]) FairTrackParP(StartPos, StartMom, StartPosErr, StartMomErr, fCharge, oc, dj, dk);
1080  // TClonesArray& clref = *fTrackParGeane;
1081  // Int_t size = clref.GetEntriesFast();
1082  // FairTrackParP *fRes = new(clref[size]) FairTrackParP();
1083  // fPro->SetPoint(vtx);
1084  // fPro->BackTrackToVirtualPlaneAtPCA(1);
1085  // Bool_t isProp = fPro->Propagate(fStart, fRes, PDGCode);
1086  // cout<<"================= %%%% ===================="<<endl;
1087  // ///----------------------------------------------------------------------
1088 
1089 
1090  //------- TEST of calculation errors -------
1091  TVector3 gPos(fRes->GetX(), fRes->GetY(), fRes->GetZ());
1092  TVector3 gMom(fRes->GetPx(), fRes->GetPy(), fRes->GetPz());
1093  TVector3 gErrPos(fRes->GetDX(), fRes->GetDY(), fRes->GetDZ());
1094  TVector3 gErrMom(fRes->GetDPx(), fRes->GetDPy(), fRes->GetDPz());
1095  // cout<<" "<<endl;
1096  if (fVerbose > 2) {
1097  // cout<<"================= %%%% ===================="<<endl;
1098 
1099  cout << "gPos:" << endl;
1100  gPos.Print();
1101  // cout<<"difference between final and initial point = "<<(-StartPos.Mag()+gPos.Mag())<<endl;
1102  cout << "gErrPos:" << endl;
1103  gErrPos.Print();
1104  // cout<<"difference between final and initial point error = "<<(-StartPosErr.Mag()+gErrPos.Mag())<<endl;
1105 
1106  cout << "gMom: " << gMom.Mag() << endl;
1107  gMom.Print();
1108  // cout<<"difference between initial and final momentum = "<<(StartMom.Mag()-gMom.Mag())<<endl;
1109  cout << "gErrMom: " << gErrMom.Mag() << endl;
1110  // cout<<"difference between initial and final momentum error= "<<(StartMomErr.Mag()-gErrMom.Mag())<<endl;
1111  gErrMom.Print();
1112  cout << "================= %%%% ====================" << endl;
1113  }
1114  //------------------------------------------
1115 
1116  if (isProp == kTRUE) {
1117  new ((*fTrackParFinal)[counterGeaneTrk]) FairTrackParH(*(fRes)); //save Track
1118  // new((*fTrackParFinal)[counterGeaneTrk]) FairTrackParP(*(fRes)); //save Track
1119  counterGeaneTrk++;
1120  cout << "***** isProp TRUE *****" << endl;
1121  } else {
1122  cout << "!!! Back-propagation with GEANE didn't return result !!!"
1123  << endl;
1124  cout << "StartPos:" << endl;
1125  StartPos.Print();
1126  cout << "StartPosErr:" << endl;
1127  StartPosErr.Print();
1128 
1129  cout << "StartMom: " << StartMom.Mag() << endl;
1130  StartMom.Print();
1131  cout << "StartMomErr: " << StartMomErr.Mag() << endl;
1132  StartMomErr.Print();
1133  new ((*fTrackParFinal)[counterGeaneTrk]) FairTrackParH(); //save NULL
1134  // new((*fTrackParFinal)[counterGeaneTrk]) FairTrackParP(); //save NULL
1135  counterGeaneTrk++;
1136  }
1137  }
1138  // fMCTracks->Delete();
1139  fMCHits->Delete();
1140  cout << "PndLmdPerformanceTask::Exec END!" << endl;*/
1141 }
1142 
1143 #include<TColor.h>
1144 #include<TStyle.h>
1145 
1146 void PndLmdPerformanceTask::Finish() {
1147  // is there a file to store the histograms to?
1148  if (hist_output_file)
1149  hist_output_file->cd();
1150 
1151  //gStyle->SetPaperSize(50,50);
1152 
1153 
1154  hist_angular_distr_acc->Divide(hist_angular_distr_gen);
1155 
1156  TGaxis::SetMaxDigits(3);
1157  TCanvas* canvas1 = new TCanvas("canvas1", "canvas1", 800, 800);
1158  canvas1->Divide(2, 2);
1159  canvas1->cd(1);
1160  hist_angular_distr_gen->Draw("COLZ");
1161  hist_angular_distr_gen->SetXTitle("#Theta [rad]");
1162  hist_angular_distr_gen->SetYTitle("#Phi [rad]");
1163  gPad->Update();
1164  canvas1->cd(2);
1165  hist_angular_distr_acc->Draw("COLZ");
1166  hist_angular_distr_acc->SetXTitle("#Theta [rad]");
1167  hist_angular_distr_acc->SetYTitle("#Phi [rad]");
1168  gPad->Update();
1169  canvas1->cd(3);
1170  hist_spatial_distr_gen->Draw("COLZ");
1171  hist_spatial_distr_gen->SetXTitle("X [cm]");
1172  hist_spatial_distr_gen->SetYTitle("Y [cm]");
1173  hist_spatial_distr_gen->GetXaxis()->SetRangeUser(-10, 10);
1174  hist_spatial_distr_gen->GetYaxis()->SetRangeUser(-10, 10);
1175  gPad->Update();
1176  canvas1->cd(4);
1177  hist_spatial_distr_acc->Draw("COLZ");
1178  hist_spatial_distr_acc->SetXTitle("X [cm]");
1179  hist_spatial_distr_acc->SetYTitle("Y [cm]");
1180  hist_spatial_distr_acc->GetXaxis()->SetRangeUser(-10, 10);
1181  hist_spatial_distr_acc->GetYaxis()->SetRangeUser(-10, 10);
1182  gPad->Update();
1183  canvas1->Print("Resolution_acceptance_results.ps(");
1184 
1185  hist_theta_over_mom_acc->Divide(hist_theta_over_mom_gen);
1186  hist_phi_over_mom_acc->Divide(hist_phi_over_mom_gen);
1187 
1188  // find the interesting momentum range
1189  int xbin_low(-1), xbin_high(-1);
1190  for (int xbin = 1; xbin <= hist_theta_over_mom_gen->GetNbinsX(); xbin++){
1191  int nentries(0);
1192  for (int ybin = 1; ybin <= hist_theta_over_mom_gen->GetNbinsY(); ybin++){
1193  nentries+=hist_angular_distr_gen->GetBinContent(xbin, ybin);
1194  }
1195  if (nentries > 0){
1196  if (xbin_low == -1){
1197  xbin_low = xbin;
1198  }
1199  } else {
1200  if (xbin_low > -1){
1201  xbin_high = xbin;
1202  break;
1203  }
1204  }
1205  }
1206  if (0 && xbin_low > -1 && xbin_high > -1){
1207  hist_theta_over_mom_gen->GetXaxis()->SetRange(xbin_low, xbin_high);
1208  hist_theta_over_mom_acc->GetXaxis()->SetRange(xbin_low, xbin_high);
1209  hist_phi_over_mom_gen->GetXaxis()->SetRange(xbin_low, xbin_high);
1210  hist_phi_over_mom_acc->GetXaxis()->SetRange(xbin_low, xbin_high);
1211  }
1212 
1213  canvas1->cd(1);
1214  hist_theta_over_mom_gen->Draw("COLZ");
1215  hist_theta_over_mom_gen->SetXTitle("p [GeV]");
1216  hist_theta_over_mom_gen->SetYTitle("#Theta [rad]");
1217  gPad->Update();
1218  canvas1->cd(2);
1219  hist_theta_over_mom_acc->Draw("COLZ");
1220  hist_theta_over_mom_acc->SetXTitle("p [GeV]");
1221  hist_theta_over_mom_acc->SetYTitle("#Theta [rad]");
1222  gPad->Update();
1223  canvas1->cd(3);
1224  hist_phi_over_mom_gen->Draw("COLZ");
1225  hist_phi_over_mom_gen->SetXTitle("p [GeV]");
1226  hist_phi_over_mom_gen->SetYTitle("#Phi [rad]");
1227  gPad->Update();
1228  canvas1->cd(4);
1229  hist_phi_over_mom_acc->Draw("COLZ");
1230  hist_phi_over_mom_acc->SetXTitle("p [GeV]");
1231  hist_phi_over_mom_acc->SetYTitle("#Phi [rad]");
1232  gPad->Update();
1233 
1234  canvas1->Print("Resolution_acceptance_results.ps(");
1235 
1236  // draw individual plane and sensor properties
1237  //gStyle->SetPaperSize(25,25);
1238  TCanvas canvas_properties_per_plane("canvas_properties_per_plane",
1239  "properties per plane", 800, 800);
1240  canvas_properties_per_plane.cd();
1241  // assuming to have 4 planes
1242  if (nplanes != 4)
1243  cout
1244  << " warning: attempting to draw histograms for a number of planes not 4! "
1245  << endl;
1246  canvas_properties_per_plane.Divide(2, 2);
1247 
1248  TCanvas canvas_properties_per_sensor("canvas_properties_per_sensor",
1249  "properties per sensor", 800, 800);
1250  canvas_properties_per_sensor.cd();
1251  // assuming to have 100 sensors per plane
1252  if (nsensors_per_plane != 100)
1253  cout
1254  << " warning: attempting to draw histograms for a number of sensors per plane not 100! "
1255  << endl;
1256  canvas_properties_per_sensor.Divide(10, 10);
1257 
1258  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1259  canvas_properties_per_plane.cd(_iplane + 1);
1260  hist_xy[_iplane]->Draw("COLZ");
1261  lmddim->Draw_Sensors(iplane);
1262  }
1263  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1264  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1265  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1266  canvas_properties_per_sensor.cd(_isensor + 1);
1267  hists_xy[_iplane][_isensor]->Draw("COLZ");
1268  }
1269  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1270  }
1271  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1272  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1273  canvas_properties_per_sensor.cd(_isensor + 1);
1274  hists_xy_local[_iplane][_isensor]->Draw("COLZ");
1275  }
1276  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1277  }
1278 
1279  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1280  canvas_properties_per_plane.cd(_iplane + 1);
1281  hist_theta_init[_iplane]->Draw("COLZ");
1282  }
1283  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1284  canvas_properties_per_plane.cd(_iplane + 1);
1285  hist_theta_rec[_iplane]->Draw("COLZ");
1286  }
1287  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1288  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1289  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1290  canvas_properties_per_sensor.cd(_isensor + 1);
1291  hists_theta_init[_iplane][_isensor]->Draw();
1292  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1293  if (_isensor > 0)
1294  hists_theta_init[_iplane][_isensor]->Draw("same E");
1295  else
1296  hists_theta_init[_iplane][_isensor]->Draw("E");
1297  //hists_theta_init[_iplane][_isensor]->SetLineColor(
1298  // kAzure - 9 + _isensor);
1299  }
1300  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1301  }
1302  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1303  canvas_properties_per_sensor.cd(_isensor + 1);
1304  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1305  if (_iplane > 0)
1306  hists_theta_init[_iplane][_isensor]->Draw("same");
1307  else
1308  hists_theta_init[_iplane][_isensor]->Draw();
1309  hists_theta_init[_iplane][_isensor]->SetLineColor(kGray + _iplane);
1310  }
1311  }
1312  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1313  gPad->Clear();
1314  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1315 
1316  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1317  canvas_properties_per_plane.cd(_iplane + 1);
1318  hist_theta_in[_iplane]->Draw("COLZ");
1319  }
1320  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1321  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1322  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1323  canvas_properties_per_sensor.cd(_isensor + 1);
1324  hists_theta_in[_iplane][_isensor]->Draw();
1325  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1326  if (_isensor > 0)
1327  hists_theta_in[_iplane][_isensor]->Draw("same E");
1328  else
1329  hists_theta_in[_iplane][_isensor]->Draw("E");
1330  //hists_theta_in[_iplane][_isensor]->SetLineColor(kAzure - 9 + _isensor);
1331  }
1332  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1333  }
1334  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1335  canvas_properties_per_sensor.cd(_isensor + 1);
1336  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1337  if (_iplane > 0)
1338  hists_theta_in[_iplane][_isensor]->Draw("same");
1339  else
1340  hists_theta_in[_iplane][_isensor]->Draw();
1341  hists_theta_in[_iplane][_isensor]->SetLineColor(kGray + _iplane);
1342  }
1343  }
1344  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1345  gPad->Clear();
1346  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1347 
1348  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1349  canvas_properties_per_plane.cd(_iplane + 1);
1350  hist_theta_diff[_iplane]->Draw("COLZ");
1351  }
1352  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1353  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1354  canvas_properties_per_plane.cd(_iplane + 1);
1355  hist_theta_rec_diff[_iplane]->Draw("COLZ");
1356  }
1357  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1358  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1359  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1360  canvas_properties_per_sensor.cd(_isensor + 1);
1361  hists_theta_diff[_iplane][_isensor]->Draw();
1362  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1363  if (_isensor > 0)
1364  hists_theta_diff[_iplane][_isensor]->Draw("same E");
1365  else
1366  hists_theta_diff[_iplane][_isensor]->Draw("E");
1367  //hists_theta_diff[_iplane][_isensor]->SetLineColor(
1368  // kAzure - 9 + _isensor);
1369  }
1370  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1371  }
1372  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1373  canvas_properties_per_sensor.cd(_isensor + 1);
1374  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1375  if (_iplane > 0)
1376  hists_theta_diff[_iplane][_isensor]->Draw("same");
1377  else
1378  hists_theta_diff[_iplane][_isensor]->Draw();
1379  hists_theta_diff[_iplane][_isensor]->SetLineColor(kGray + _iplane);
1380  }
1381  }
1382  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1383  gPad->Clear();
1384  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1385 
1386  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1387  canvas_properties_per_plane.cd(_iplane + 1);
1388  hist_theta_diff_rel[_iplane]->Draw("COLZ");
1389  }
1390  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1391  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1392  canvas_properties_per_plane.cd(_iplane + 1);
1393  hist_theta_rec_diff_rel[_iplane]->Draw("COLZ");
1394  }
1395  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps(");
1396  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1397  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1398  canvas_properties_per_sensor.cd(_isensor + 1);
1399  hists_theta_diff_rel[_iplane][_isensor]->Draw();
1400  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1401  if (_isensor > 0)
1402  hists_theta_diff_rel[_iplane][_isensor]->Draw("same E");
1403  else
1404  hists_theta_diff_rel[_iplane][_isensor]->Draw("E");
1405  //hists_theta_diff_rel[_iplane][_isensor]->SetLineColor(
1406  // kAzure - 9 + _isensor);
1407  }
1408  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1409  }
1410  for (int _isensor = 0; _isensor < nsensors_per_plane; _isensor++) {
1411  canvas_properties_per_sensor.cd(_isensor + 1);
1412  for (int _iplane = 0; _iplane < nplanes; _iplane++) {
1413  if (_iplane > 0)
1414  hists_theta_diff_rel[_iplane][_isensor]->Draw("same");
1415  else
1416  hists_theta_diff_rel[_iplane][_isensor]->Draw();
1417  hists_theta_diff_rel[_iplane][_isensor]->SetLineColor(kGray + _iplane);
1418  }
1419  }
1420  canvas_properties_per_sensor.cd(nsensors_per_plane+1);
1421  gPad->Clear();
1422  canvas_properties_per_sensor.Print("Resolution_acceptance_results.ps(");
1423 
1424  canvas_properties_per_plane.Clear();
1425  canvas_properties_per_plane.Print("Resolution_acceptance_results.ps)");
1426  // convert to pdf
1427  cout << "converting to pdf " << system(
1428  "ps2pdf Resolution_acceptance_results.ps") << endl;
1429  cout << "deleting ps file " << system(
1430  "rm Resolution_acceptance_results.ps") << endl;
1431 
1432  if (hist_output_file)
1433  hist_output_file->Write();
1434  //output_histfile->Write();
1435 
1436 }
1437 
1438 std::map<int, std::vector<int> > PndLmdPerformanceTask::AssignHitsToTracks() {
1439  /*
1440  std::map<int, std::vector<int> > result;
1441  for (int i = 0; i < fMCHits->GetEntriesFast(); i++) { //get all MC Hits
1442  PndSdsMCPoint* myPoint = (PndSdsMCPoint*) (fMCHits->At(i)); //sort MCHits with Tracks
1443  //PndMCTrack* myTrack = (PndMCTrack*)(fMCTracks->At(myPoint->GetTrackID()));
1444  result[myPoint->GetTrackID()].push_back(i);
1445 
1446  }
1447  return result;*/
1448 }
1449 
1450 void PndLmdPerformanceTask::DrawProgressBar(int len, double percent) {
1451  if ((int) (last_percent * 100) == (int) (percent * 100))
1452  return;
1453  //cout << " drawing " << endl;
1454  cout << "\x1B[2K"; // Erase the entire current line.
1455  cout << "\x1B[0E"; // Move to the beginning of the current line.
1456  string progress;
1457  for (int i = 0; i < len; ++i) {
1458  if (i < static_cast<int> (len * percent)) {
1459  progress += "=";
1460  } else {
1461  progress += " ";
1462  }
1463  }
1464  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
1465  << "%";
1466  flush( cout); // Required.
1467  last_percent = percent;
1468 }
1469 
1470 ClassImp(PndLmdPerformanceTask)
1471 ;
nplanes
const int nplanes(120)
PndLmdPerformanceTask::DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: PndLmdPerformanceTask.cxx:1450
i
Int_t i
Definition: run_full.C:25
PndLmdPerformanceTask::hists_xy_local
TH2 * hists_xy_local[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:213
PndLmdPerformanceTask::hists_theta_diff_rel
TH1 * hists_theta_diff_rel[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:218
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndLmdDim::Instance
static PndLmdDim * Instance()
Definition: PndLmdDim.cxx:249
ClassImp
ClassImp(PndLmdPerformanceTask)
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndLmdPerformanceTask::Exec
virtual void Exec(Option_t *opt)
Definition: PndLmdPerformanceTask.cxx:585
PndLmdPerformanceTask::fTrackMCVectMapIn
std::map< int, std::vector< TVector3 > > fTrackMCVectMapIn
Definition: PndLmdPerformanceTask.h:101
PndLmdPerformanceTask::hists_theta_in
TH1 * hists_theta_in[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:216
PndSdsMCPoint::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
PndMCTrack::IsGeneratorCreated
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
PndGenfitField
Magnetic field.
Definition: PndGenfitField.h:12
RKTrackRep
Track Representation module based on a Runge-Kutta algorithm including a full material model...
Definition: tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h:49
nHits
int nHits
Definition: RiemannTest.C:16
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndLmdPerformanceTask::SetHistFilename
void SetHistFilename(TString filename)
Definition: PndLmdPerformanceTask.cxx:527
PndLmdPerformanceTask::fTrackMCVectMapInit
std::vector< TVector3 > fTrackMCVectMapInit
Definition: PndLmdPerformanceTask.h:102
nsensors_per_plane
const unsigned int nsensors_per_plane(100)
PndLmdPerformanceTask::hist_theta_rec_diff_rel
TH2 * hist_theta_rec_diff_rel[nplanes]
Definition: PndLmdPerformanceTask.h:209
PndGeoHandling::Instance
static PndGeoHandling * Instance()
Definition: PndGeoHandling.cxx:36
RKTrackRep::extrapolateToPoint
void extrapolateToPoint(const TVector3 &pos, TVector3 &poca, TVector3 &dirInPoca)
This method is to extrapolate the track to point of closest approach to a point in space...
Definition: RKTrackRep.cxx:329
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
GFFieldManager::init
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
Definition: GFFieldManager.h:65
PndLmdPerformanceTask::Init
virtual InitStatus Init()
Definition: PndLmdPerformanceTask.cxx:84
PndLmdPerformanceTask::hists_xy
TH2 * hists_xy[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:212
PndLmdPerformanceTask::fTrackPlaneHitMap
std::map< int, std::vector< TVector3 > > fTrackPlaneHitMap
Definition: PndLmdPerformanceTask.h:99
PndLmdPerformanceTask::hists_theta_diff
TH1 * hists_theta_diff[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:217
GFFieldManager::getInstance
static GFFieldManager * getInstance()
Definition: GFFieldManager.h:69
PndLmdPerformanceTask::hist_theta_diff_rel
TH2 * hist_theta_diff_rel[nplanes]
Definition: PndLmdPerformanceTask.h:207
PndLmdPerformanceTask::hists_theta_init
TH1 * hists_theta_init[nplanes][nsensors_per_plane]
Definition: PndLmdPerformanceTask.h:215
PndLmdPerformanceTask::AssignHitsToTracks
std::map< int, std::vector< int > > AssignHitsToTracks()
Definition: PndLmdPerformanceTask.cxx:1438
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndLmdPerformanceTask::hist_theta_rec_diff
TH2 * hist_theta_rec_diff[nplanes]
Definition: PndLmdPerformanceTask.h:208
filename
const string filename
Definition: Modify_root_geometry.C:23
Generated on Sun Aug 18 2019 03:01:36 for FairRoot/PandaRoot by   doxygen 1.8.8