All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
CreatePDFProcessor.cc
Go to the documentation of this file.
1 #include <cmath>
2 #include <vector>
3 
4 #include <marlin/Global.h>
5 #include "DD4hep/Detector.h"
6 #include "DDRec/Vector3D.h"
7 #include "DD4hep/DD4hepUnits.h"
8 
9 #include "EVENT/ReconstructedParticle.h"
10 #include "EVENT/LCRelation.h"
11 #include "EVENT/MCParticle.h"
12 #include "IMPL/ClusterImpl.h"
13 #include <lcio.h>
14 
15 #include "TVector3.h"
16 
17 #include "CreatePDFProcessor.hh"
18 
19 CreatePDFProcessor aCreatePDFProcessor ;
20 
21 CreatePDFProcessor::CreatePDFProcessor()
22  : Processor("CreatePDFProcessor") {
23 
24  // Processor description
25  _description = "Processor for PDF templetes creation";
26 
27 
28  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
29  "OutputPDFFilename" ,
30  "Input collection of Reconstrcuted Particle",
31  _filename,
32  std::string("pdf_standard_v01.root"));
33 
34  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
35  "RecoParticleCollection" ,
36  "Input collection of Reconstrcuted Particle",
37  _PfoCollection,
38  std::string("PandoraPFOs"));
39 
40  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
41  "MCTruthLinkCollection" ,
42  "Input collection of Reconstrcuted Particle",
43  _LinkCollection,
44  std::string("RecoMCTruthLink"));
45 
46  std::vector< float > parselectron,parsmuon,parspion,parskaon,parsproton;
47  parselectron.push_back(-2.40638e-03);
48  parselectron.push_back(7.10337e-01);
49  parselectron.push_back(2.87718e-01);
50  parselectron.push_back(-1.71591e+00);
51  parselectron.push_back(0.0);
52  registerProcessorParameter( "dEdxParameter_electron" ,
53  "dEdx Parameters for electron",
54  _dEdxParamsElectron,
55  parselectron);
56 
57  parsmuon.push_back(8.11408e-02);
58  parsmuon.push_back(9.92207e-01);
59  parsmuon.push_back(7.58509e+05);
60  parsmuon.push_back(-1.70167e-01);
61  parsmuon.push_back(4.63670e-04);
62  registerProcessorParameter( "dEdxParameter_muon" ,
63  "dEdx Parameters for muon",
64  _dEdxParamsMuon,
65  parsmuon);
66 
67  parspion.push_back(8.10756e-02);
68  parspion.push_back(-1.45051e+06);
69  parspion.push_back(-3.09843e+04);
70  parspion.push_back(2.84056e-01);
71  parspion.push_back(3.38131e-04);
72  registerProcessorParameter( "dEdxParameter_pion" ,
73  "dEdx Parameters for pion",
74  _dEdxParamsPion,
75  parspion);
76 
77  parskaon.push_back(7.96117e-02);
78  parskaon.push_back(4.13335e+03);
79  parskaon.push_back(1.13577e+06);
80  parskaon.push_back(1.80555e-01);
81  parskaon.push_back(-3.15083e-04);
82  registerProcessorParameter( "dEdxParameter_kaon" ,
83  "dEdx Parameters for Kaon",
84  _dEdxParamsKaon,
85  parskaon);
86 
87  parsproton.push_back(7.78772e-02);
88  parsproton.push_back(4.49300e+04);
89  parsproton.push_back(9.13778e+04);
90  parsproton.push_back(1.50088e-01);
91  parsproton.push_back(-6.64184e-04);
92  registerProcessorParameter( "dEdxParameter_proton" ,
93  "dEdx Parameters for proton",
94  _dEdxParamsProton,
95  parsproton);
96 
97  registerProcessorParameter( "dEdxNormalization" ,
98  "dEdx Normalization",
99  _dEdxNormalization,
100  float(1.350e-7));
101 
102  registerProcessorParameter( "dEdxErrorFactor" ,
103  "dEdx Error factor",
104  _dEdxErrorFactor,
105  float(7.55));
106 }
107 
108 void CreatePDFProcessor::init() {
109  //streamlog_out(DEBUG) << " init called " << std::endl ;
110  //dEdx parameters
111  double pars[28];
112  for(int i=0;i<5;i++) pars[i]=_dEdxParamsElectron[i];
113  for(int i=0;i<5;i++) pars[i+5]=_dEdxParamsMuon[i];
114  for(int i=0;i<5;i++) pars[i+10]=_dEdxParamsPion[i];
115  for(int i=0;i<5;i++) pars[i+15]=_dEdxParamsKaon[i];
116  for(int i=0;i<5;i++) pars[i+20]=_dEdxParamsProton[i];
117  pars[25] =true; //(double) _UseBayes;
118  pars[26] =(double) _dEdxNormalization;
119  pars[27] =(double) _dEdxErrorFactor;
120 
121  //create histograms
122  std::string histtxt,labeltxt;
123  for(int i=0;i<6;i++){
124  for(int j=0;j<21;j++){
125  switch(j){
126  case 0:
127  //momentum
128  histtxt="hmomentum" + itos(i+1);
129  labeltxt="p(track)";
130  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",50,0.0,250.0);
131  pidvariable[i][j]->SetMarkerStyle(20);
132  pidvariable[i][j]->SetFillColor(0);
133  pidvariable[i][j]->SetLineColor(1);
134  pidvariable[i][j]->SetLineWidth(2.0);
135  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
136  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
137  pidvariable[i][j]->SetMinimum(0.0);
138  break;
139  case 1:
140  //ep
141  histtxt="hep" + itos(i+1);
142  labeltxt="E/p";
143  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",40,0.0,2.0);
144  pidvariable[i][j]->SetMarkerStyle(20);
145  pidvariable[i][j]->SetFillColor(0);
146  pidvariable[i][j]->SetLineColor(1);
147  pidvariable[i][j]->SetLineWidth(2.0);
148  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
149  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
150  pidvariable[i][j]->SetMinimum(0.0);
151  break;
152  case 2:
153  //ehad
154  histtxt="hehad" + itos(i+1);
155  labeltxt="ECAL/(ECAL+HCAL)";
156  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",50,0.0,1.0);
157  pidvariable[i][j]->SetMarkerStyle(20);
158  pidvariable[i][j]->SetFillColor(0);
159  pidvariable[i][j]->SetLineColor(1);
160  pidvariable[i][j]->SetLineWidth(2.0);
161  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
162  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
163  pidvariable[i][j]->SetMinimum(0.0);
164  break;
165  case 3:
166  //photonEnergy(reserved)
167  histtxt="hphotonEnergy" + itos(i+1);
168  labeltxt="E(photon) (GeV)";
169  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",50,0.0,50.0);
170  pidvariable[i][j]->SetMarkerStyle(20);
171  pidvariable[i][j]->SetFillColor(0);
172  pidvariable[i][j]->SetLineColor(1);
173  pidvariable[i][j]->SetLineWidth(2.0);
174  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
175  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
176  pidvariable[i][j]->SetMinimum(0.0);
177  break;
178  case 4:
179  //mucal
180  histtxt="hmucal" + itos(i+1);
181  labeltxt="MuCAL (GeV)";
182  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,20.0);
183  pidvariable[i][j]->SetMarkerStyle(20);
184  pidvariable[i][j]->SetFillColor(0);
185  pidvariable[i][j]->SetLineColor(1);
186  pidvariable[i][j]->SetLineWidth(2.0);
187  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
188  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
189  pidvariable[i][j]->SetMinimum(0.0);
190  break;
191  case 5:
192  //chi2
193  histtxt="hchi2" + itos(i+1);
194  labeltxt="#chi^{2}";
195  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",104,-2.0,50.0);
196  pidvariable[i][j]->SetMarkerStyle(20);
197  pidvariable[i][j]->SetFillColor(0);
198  pidvariable[i][j]->SetLineColor(1);
199  pidvariable[i][j]->SetLineWidth(2.0);
200  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
201  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
202  pidvariable[i][j]->SetMinimum(0.0);
203  break;
204  case 6:
205  //shower max/expected shower max
206  histtxt="hldiscrepancy" + itos(i+1);
207  labeltxt="ShowerMax/EM ShowerMax";
208  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",65,-30.0,100.0);
209  pidvariable[i][j]->SetMarkerStyle(20);
210  pidvariable[i][j]->SetFillColor(0);
211  pidvariable[i][j]->SetLineColor(1);
212  pidvariable[i][j]->SetLineWidth(2.0);
213  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
214  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
215  pidvariable[i][j]->SetMinimum(0.0);
216  break;
217  case 7:
218  //absorption length/expected shower max
219  histtxt="htdiscrepancy" + itos(i+1);
220  labeltxt="absorption Length/Rm";
221  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,2.0);
222  pidvariable[i][j]->SetMarkerStyle(20);
223  pidvariable[i][j]->SetFillColor(0);
224  pidvariable[i][j]->SetLineColor(1);
225  pidvariable[i][j]->SetLineWidth(2.0);
226  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
227  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
228  pidvariable[i][j]->SetMinimum(0.0);
229  break;
230  case 8:
231  //shower max using hits
232  histtxt="hxl20" + itos(i+1);
233  labeltxt="xl20";
234  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,100.0);
235  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
236  pidvariable[i][j]->SetMarkerStyle(20);
237  pidvariable[i][j]->SetFillColor(0);
238  pidvariable[i][j]->SetLineColor(1);
239  pidvariable[i][j]->SetLineWidth(2.0);
240  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
241  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
242  pidvariable[i][j]->SetMinimum(0.0);
243  break;
244  case 9:
245  //likelihood using electron hypothesis
246  histtxt="hlikeliele" + itos(i+1);
247  labeltxt="Log(likelihood) (electron)";
248  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-100.0,0.0);
249  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
250  pidvariable[i][j]->SetMarkerStyle(20);
251  pidvariable[i][j]->SetFillColor(0);
252  pidvariable[i][j]->SetLineColor(1);
253  pidvariable[i][j]->SetLineWidth(2.0);
254  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
255  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
256  pidvariable[i][j]->SetMinimum(0.0);
257  break;
258  case 10:
259  //likelihood using muon hypothesis
260  histtxt="hlikelimuo" + itos(i+1);
261  labeltxt="Log(likelihood) (#mu)";
262  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-100.0,0.0);
263  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
264  pidvariable[i][j]->SetMarkerStyle(20);
265  pidvariable[i][j]->SetFillColor(0);
266  pidvariable[i][j]->SetLineColor(1);
267  pidvariable[i][j]->SetLineWidth(2.0);
268  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
269  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
270  pidvariable[i][j]->SetMinimum(0.0);
271  break;
272  case 11:
273  //likelihood using pion hypothesis
274  histtxt="hlikelipi" + itos(i+1);
275  labeltxt="Log(likelihood) (#pi)";
276  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-100.0,0.0);
277  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
278  pidvariable[i][j]->SetMarkerStyle(20);
279  pidvariable[i][j]->SetFillColor(0);
280  pidvariable[i][j]->SetLineColor(1);
281  pidvariable[i][j]->SetLineWidth(2.0);
282  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
283  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
284  pidvariable[i][j]->SetMinimum(0.0);
285  break;
286  case 12:
287  //likelihood using kaon hypothesis
288  histtxt="hlikelik" + itos(i+1);
289  labeltxt="Log(likelihood) (K)";
290  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-100.0,0.0);
291  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
292  pidvariable[i][j]->SetMarkerStyle(20);
293  pidvariable[i][j]->SetFillColor(0);
294  pidvariable[i][j]->SetLineColor(1);
295  pidvariable[i][j]->SetLineWidth(2.0);
296  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
297  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
298  pidvariable[i][j]->SetMinimum(0.0);
299  break;
300  case 13:
301  //likelihood using proton hypothesis
302  histtxt="hlikelip" + itos(i+1);
303  labeltxt="Log(likelihood) (p)";
304  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-100.0,0.0);
305  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
306  pidvariable[i][j]->SetMarkerStyle(20);
307  pidvariable[i][j]->SetFillColor(0);
308  pidvariable[i][j]->SetLineColor(1);
309  pidvariable[i][j]->SetLineWidth(2.0);
310  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
311  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
312  pidvariable[i][j]->SetMinimum(0.0);
313  break;
314  case 14:
315  //likelihood ratio
316  histtxt="hlikelikpi" + itos(i+1);
317  labeltxt="Log(likelihood(K)/likelihood(#pi))";
318  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",160,-50.0,30.0);
319  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
320  pidvariable[i][j]->SetMarkerStyle(20);
321  pidvariable[i][j]->SetFillColor(0);
322  pidvariable[i][j]->SetLineColor(1);
323  pidvariable[i][j]->SetLineWidth(2.0);
324  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
325  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
326  pidvariable[i][j]->SetMinimum(0.0);
327  break;
328  case 15:
329  //likelihood ratio
330  histtxt="hlikelipk" + itos(i+1);
331  labeltxt="Log(likelihood(p)/likelihood(K))";
332  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",160,-50.0,30.0);
333  //pidvariable[i][j]=new TH1F(histtxt.c_str(),"",80,0.0,40.0);
334  pidvariable[i][j]->SetMarkerStyle(20);
335  pidvariable[i][j]->SetFillColor(0);
336  pidvariable[i][j]->SetLineColor(1);
337  pidvariable[i][j]->SetLineWidth(2.0);
338  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
339  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
340  pidvariable[i][j]->SetMinimum(0.0);
341  break;
342  case 16:
343  //delta x
344  histtxt="hdeltax" + itos(i+1);
345  labeltxt="#Deltax#timesQ";
346  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-5.0,5.0);
347  pidvariable[i][j]->SetMarkerStyle(20);
348  pidvariable[i][j]->SetFillColor(0);
349  pidvariable[i][j]->SetLineColor(1);
350  pidvariable[i][j]->SetLineWidth(2.0);
351  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
352  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
353  pidvariable[i][j]->SetMinimum(0.0);
354  break;
355  case 17:
356  //delta z
357  histtxt="hdeltaz" + itos(i+1);
358  labeltxt="#Deltaz";
359  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,-5.0,5.0);
360  pidvariable[i][j]->SetMarkerStyle(20);
361  pidvariable[i][j]->SetFillColor(0);
362  pidvariable[i][j]->SetLineColor(1);
363  pidvariable[i][j]->SetLineWidth(2.0);
364  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
365  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
366  pidvariable[i][j]->SetMinimum(0.0);
367  break;
368  case 18:
369  //cluster depth
370  histtxt="hcludepth" + itos(i+1);
371  labeltxt="cluster depth";
372  pidvariable[i][j]=new TH1F(histtxt.c_str(),"",200,1500.0,3500.0);
373  pidvariable[i][j]->SetMarkerStyle(20);
374  pidvariable[i][j]->SetFillColor(0);
375  pidvariable[i][j]->SetLineColor(1);
376  pidvariable[i][j]->SetLineWidth(2.0);
377  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
378  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
379  pidvariable[i][j]->SetMinimum(0.0);
380  break;
381  case 19:
382  //hit mean
383  histtxt="hhitmean" + itos(i+1);
384  labeltxt="hit mean";
385  if(i==0) pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,2500.0);
386  if(i!=0) pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,1000.0);
387  pidvariable[i][j]->SetMarkerStyle(20);
388  pidvariable[i][j]->SetFillColor(0);
389  pidvariable[i][j]->SetLineColor(1);
390  pidvariable[i][j]->SetLineWidth(2.0);
391  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
392  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
393  pidvariable[i][j]->SetMinimum(0.0);
394  break;
395  case 20:
396  //hitrms
397  histtxt="hhitrms" + itos(i+1);
398  labeltxt="hitrms";
399  if(i==0) pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,2500.0);
400  if(i!=0) pidvariable[i][j]=new TH1F(histtxt.c_str(),"",100,0.0,1000.0);
401  pidvariable[i][j]->SetMarkerStyle(20);
402  pidvariable[i][j]->SetFillColor(0);
403  pidvariable[i][j]->SetLineColor(1);
404  pidvariable[i][j]->SetLineWidth(2.0);
405  pidvariable[i][j]->GetYaxis()->SetTitle("Normalized Events");
406  pidvariable[i][j]->GetXaxis()->SetTitle(labeltxt.c_str());
407  pidvariable[i][j]->SetMinimum(0.0);
408  break;
409  }
410  if(i==0) pidvariable[i][j]->SetLineColor(5);
411  if(i==1) pidvariable[i][j]->SetLineColor(6);
412  if(i==2) pidvariable[i][j]->SetLineColor(2);
413  if(i==3) pidvariable[i][j]->SetLineColor(3);
414  if(i==4) pidvariable[i][j]->SetLineColor(4);
415  if(i==5) pidvariable[i][j]->SetLineColor(7);
416  }
417  }
418 
419  _myPID = new LikelihoodPID(pars);
420 
421  //get B field
422  double bfield[3]={};
423  dd4hep::Detector& lcdd = dd4hep::Detector::getInstance();
424  lcdd.field().magneticField({0.0,0.0,0.0}, bfield);
425  _bfield = (float)bfield[2]/dd4hep::tesla;
426 
427  printParameters();
428 }
429 
430 void CreatePDFProcessor::processRunHeader( LCRunHeader* /*run*/) {
431 }
432 
433 void CreatePDFProcessor::processEvent( LCEvent * evt ) {
434  _PFOCol = evt->getCollection( _PfoCollection ) ;
435  _LinkCol = evt->getCollection( _LinkCollection ) ;
436  int npfo = _PFOCol->getNumberOfElements();
437  int nln = _LinkCol->getNumberOfElements();
438 
439  std::map<int, int> pid;
440  pid.insert(std::map<int, int>::value_type(11, 0));
441  pid.insert(std::map<int, int>::value_type(13, 1));
442  pid.insert(std::map<int, int>::value_type(211, 2));
443  pid.insert(std::map<int, int>::value_type(321, 3));
444  pid.insert(std::map<int, int>::value_type(2212, 4));
445 
446  for(int i=0;i<npfo;i++){
447  ReconstructedParticle* part = (ReconstructedParticle*) _PFOCol->getElementAt(i);
448  if(part->getCharge()==0) continue; //neutral
449 
450  EVENT::Track* trk = part->getTracks()[0];
451  EVENT::ClusterVec cluvec = part->getClusters();
452 
453  TVector3 pp = part->getMomentum();
454 
455  //get deposit energy
456  double ecal=0.0,hcal=0.0,mucal=0.0;
457  if(cluvec.size()!=0){
458  for(unsigned int i=0;i<cluvec.size();i++){
459  ecal+=cluvec[i]->getSubdetectorEnergies()[0];
460  hcal+=cluvec[i]->getSubdetectorEnergies()[1];
461  mucal+=cluvec[i]->getSubdetectorEnergies()[2];
462  }
463  }
464 
465  //check MC particle
466  int pdfid=-1;
467  for(int j=0;j< nln;j++){
468  LCRelation* plc = (LCRelation*) _LinkCol->getElementAt( j );
469  ReconstructedParticle* fp = (ReconstructedParticle*) plc->getFrom();
470 
471  if(part == fp){
472  MCParticle* partmc = (MCParticle*) plc->getTo();
473  int mcpdg = fabs(partmc->getPDG());
474  if(mcpdg!=11 && mcpdg!=13 && mcpdg!=211 && mcpdg!=321 && mcpdg!=2212) continue;
475  pdfid = pid[mcpdg];
476  if(pdfid==1 && mucal==0.0) pdfid=5;
477  break;
478  }
479  } //done.
480 
481  if(pdfid<0) continue; //do not match to MC particle
482 
483  //start to fill historams
484  //momentum
485  pidvariable[pdfid][0]->Fill(pp.Mag());
486 
487  if(cluvec.size()!=0){
488  pidvariable[pdfid][1]->Fill((ecal+hcal)/pp.Mag());
489  pidvariable[pdfid][2]->Fill(ecal/(ecal+hcal));
490  pidvariable[pdfid][4]->Fill(mucal);
491 
492  //get shower profiles
493  FloatVec shapes=cluvec[0]->getShape();
494  pidvariable[pdfid][5]->Fill(shapes[0+4]);
495  pidvariable[pdfid][6]->Fill(shapes[5+4]);
496  pidvariable[pdfid][7]->Fill(fabs(shapes[3+4])/shapes[6+4]);
497  pidvariable[pdfid][8]->Fill(fabs(shapes[15+4])/7.00);
498 
499  //for mu/pi<5.0GeV?
500  if(pp.Mag()>0.3 && pp.Mag()<6.0){
501  pidvariable[pdfid][18]->Fill(shapes[17+4]);
502  pidvariable[pdfid][19]->Fill(shapes[18+4]);
503  pidvariable[pdfid][20]->Fill(shapes[19+4]);
504  }
505  }
506 
507  //get dedx info
508  double dedx = trk->getdEdx();
509  double dedxerr = trk->getdEdxError();
510  double ehypo = -0.5*fabs(_myPID->get_dEdxChi2(0, pp, dedxerr, dedx))+_myPID->get_dEdxFactor(0,pp, dedxerr, dedx);
511  double muhypo = -0.5*fabs(_myPID->get_dEdxChi2(1, pp, dedxerr, dedx))+_myPID->get_dEdxFactor(1, pp, dedxerr, dedx);
512  double pihypo = -0.5*fabs(_myPID->get_dEdxChi2(2, pp, dedxerr, dedx))+_myPID->get_dEdxFactor(2, pp, dedxerr, dedx);
513  double khypo = -0.5*fabs(_myPID->get_dEdxChi2(3, pp, dedxerr, dedx))+_myPID->get_dEdxFactor(3, pp, dedxerr, dedx);
514  double phypo = -0.5*fabs(_myPID->get_dEdxChi2(4, pp, dedxerr, dedx))+_myPID->get_dEdxFactor(4, pp, dedxerr, dedx);
515 
516  pidvariable[pdfid][9]->Fill(ehypo);
517  pidvariable[pdfid][10]->Fill(muhypo);
518  pidvariable[pdfid][11]->Fill(pihypo);
519  pidvariable[pdfid][12]->Fill(khypo);
520  pidvariable[pdfid][13]->Fill(phypo);
521 
522  pidvariable[pdfid][14]->Fill(std::log(exp(khypo)/exp(pihypo)));
523  pidvariable[pdfid][15]->Fill(std::log(exp(phypo)/exp(khypo)));
524 
525  //deltax and deltaz
526  float delpos[3]={0.0,0.0,0.0};
527  if(cluvec.size()!=0){ // && pp.Mag()<6.0){
528  CalculateDeltaPosition(part->getCharge(), pp, cluvec[0]->getPosition(), delpos);
529  pidvariable[pdfid][16]->Fill(delpos[0]);
530  pidvariable[pdfid][17]->Fill(delpos[2]);
531  }
532  }
533 }
534 
535 void CreatePDFProcessor::check( LCEvent * /*evt*/ ) {
536 }
537 
538 void CreatePDFProcessor::end() {
539  //delete _myShowerShapes;
540  //open tfile
541  _fpdf = new TFile(_filename.c_str(), "RECREATE");
542 
543  for(int i=0;i<6;i++){
544  for(int j=0;j<21;j++){
545  pidvariable[i][j]->Write();
546  }
547  }
548 
549  _fpdf->Close();
550 
551  delete _fpdf;
552  delete _myPID;
553 }
554 
555 void CreatePDFProcessor::CalculateDeltaPosition(float charge, TVector3& p, const float* calpos, float* delpos){
556  //calculate some variables for preparation
557  //transverse momentum
558  double prphi=sqrt(pow(p[0],2)+pow(p[1],2));
559  //momentum
560  double pp=p.Mag();
561  //radius in r-phi plane
562  double radius=prphi/(0.3*_bfield);
563  //tangent lambda
564  double tanlam=p[2]/prphi;
565  //change to unit vector
566  p[0]=p[0]/prphi;
567  p[1]=p[1]/prphi;
568  p[2]=p[2]/pp;
569 
570  //chenge position vector to meter
571  float calcalpos[3];
572  calcalpos[0]=calpos[0]/1000.0;
573  calcalpos[1]=calpos[1]/1000.0;
574  calcalpos[2]=calpos[2]/1000.0;
575 
576 
577  //radius at the tracker end
578  double rradius=sqrt(pow(calcalpos[0],2)+pow(calcalpos[1],2));
579 
580  //cout << "check val. " << prphi << " " << radius << " " << tanlam
581  // << " " << rradius << endl;
582 
583  //cal. the position of the center of track circle
584  TVector3 cc;
585  cc[0]=-charge*p[1]*prphi/(0.3*_bfield);
586  cc[1]=charge*p[0]*prphi/(0.3*_bfield);
587  cc[2]=radius*tanlam;
588 
589  //cal. sign and cosine 2theta
590  double sintheta=charge*rradius/(2*radius);
591  double costheta=sqrt(1-sintheta*sintheta);
592  double sin2theta=2*sintheta*costheta;
593  double cos2theta=costheta*costheta-sintheta*sintheta;
594 
595  TVector3 calpos2;
596  calpos2[2]=rradius*tanlam;
597  calpos2[0]=cc[0]*(1-cos2theta)+cc[1]*sin2theta;
598  calpos2[1]=-cc[0]*sin2theta+cc[1]*(1-cos2theta);
599 
600  //cal. difference of the position
601  //float delpos[3];
602  delpos[0]=charge*fabs(calcalpos[0]-calpos2[0]);
603  delpos[1]=charge*fabs(calcalpos[1]-calpos2[1]);
604  delpos[2]=fabs(calcalpos[2]-calpos2[2]);
605 
606 
607  //cout << "calpos: " << calcalpos[0] << " " << calcalpos[1] << " " << calcalpos[2] << endl;
608  //cout << "calpos2: " << calpos2[0] << " " << calpos2[1] << " " << calpos2[2] << endl;
609  //cout << "result: " << delpos[0] << " " << delpos[1] << " " << delpos[2] << endl;
610  return;
611 }
612 
CreatePDFProcessor::_dEdxParamsMuon
std::vector< float > _dEdxParamsMuon
Definition: CreatePDFProcessor.hh:43
CreatePDFProcessor::_myPID
LikelihoodPID * _myPID
Definition: CreatePDFProcessor.hh:50
CreatePDFProcessor::_dEdxParamsProton
std::vector< float > _dEdxParamsProton
Definition: CreatePDFProcessor.hh:46
LikelihoodPID::get_dEdxChi2
double get_dEdxChi2(int parttype, TVector3 p, float hit, double dEdx)
Definition: LikelihoodPID.cc:369
CreatePDFProcessor::_dEdxParamsElectron
std::vector< float > _dEdxParamsElectron
Definition: CreatePDFProcessor.hh:42
LikelihoodPID::get_dEdxFactor
double get_dEdxFactor(int parttype, TVector3 p, float hit, double dEdx)
Definition: LikelihoodPID.cc:411
CreatePDFProcessor::CalculateDeltaPosition
void CalculateDeltaPosition(float charge, TVector3 &p, const float *caylpos, float *delpos)
Definition: CreatePDFProcessor.cc:555
CreatePDFProcessor::pidvariable
TH1F * pidvariable[6][21]
Definition: CreatePDFProcessor.hh:52
CreatePDFProcessor::itos
std::string itos(int i)
Definition: CreatePDFProcessor.hh:55
CreatePDFProcessor::_dEdxParamsKaon
std::vector< float > _dEdxParamsKaon
Definition: CreatePDFProcessor.hh:45
CreatePDFProcessor::_LinkCol
LCCollection * _LinkCol
Definition: CreatePDFProcessor.hh:40
CreatePDFProcessor::_dEdxParamsPion
std::vector< float > _dEdxParamsPion
Definition: CreatePDFProcessor.hh:44
CreatePDFProcessor::check
virtual void check(LCEvent *evt)
Definition: CreatePDFProcessor.cc:535
CreatePDFProcessor::processEvent
virtual void processEvent(LCEvent *evt)
Definition: CreatePDFProcessor.cc:433
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
aCreatePDFProcessor
CreatePDFProcessor aCreatePDFProcessor
Definition: CreatePDFProcessor.cc:19
CreatePDFProcessor::_PFOCol
LCCollection * _PFOCol
Definition: CreatePDFProcessor.hh:39
CreatePDFProcessor::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Definition: CreatePDFProcessor.cc:430