StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
digparticle.cxx
1 // //
3 // //
4 // DIGParticle //
5 // //
6 // particle class //
7 // contains: -entry and exit point of the particle into the plane //
8 // -segment list (Charge, position) of the track //
9 // -pixel list (number, charge) where charge has been collected //
10 // //
11 // //
13 #include <digparticle.h>
14 
15 #include <TROOT.h> // for gROOT object
16 #include <TMath.h>
17 #include <TMatrixD.h>
18 #include <TCanvas.h>
19 #include <TGraph.h>
20 #include <TAxis.h>
21 #include <TRandom3.h>
22 #include <TFile.h>
23 #include <TTree.h>
24 #include <TBranch.h>
25 #include <TClonesArray.h>
26 #include <TF2.h>
27 #include <TProfile2D.h>
28 #include <TH2.h>
29 
30 
31 #include "digplane.h"
32 #include "digadc.h"
33 #include "digtransport.h"
34 
35 using namespace std;
36 
37 extern Int_t GlobalSeed;
38 
39 //_______________________________________________________________________________________
40 //
41  Double_t Lorentz2D(Double_t *x, Double_t *par){
42  //x[0] = x
43  //x[1] = y
44  // par[0] = Gamma
45  // par[1] = x0
46  // par[2] = y0
47  // par[3] = norm
48  if(par[0]>0){
49  Double_t Pi = 3.141592653;
50  return par[3]*par[0]/Pi/((x[0]-par[1])*(x[0]-par[1])+(x[1]-par[2])*(x[1]-par[2])+par[0]*par[0]) ;
51  }else{
52  return 0;
53  }
54 }
55 //_______________________________________________________________________________________
56 //
57 Double_t SumGaus2D(Double_t *x, Double_t *par){
58  //par[0] Norm_1
59  //par[1] x0_1
60  //par[2] sigma_x_1
61  //par[3] y0_1
62  //par[4] sigma_y_1
63  //par[5] Norm_2
64  //par[6] x0_2
65  //par[7] sigma_x_2
66  //par[8] y0_2
67  //par[9] sigma_y_2
68  if((par[2]!=0.0) && (par[4]!=0.0) && (par[7]!=0.0) && (par[9]!=0.0)){
69  double rx = (x[0]-par[1])/par[2];
70  double ry = (x[1]-par[3])/par[4];
71  double rx2 = (x[0]-par[6])/par[7];
72  double ry2 = (x[1]-par[8])/par[9];
73  return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)+par[5]*TMath::Exp(-(rx2*rx2+ry2*ry2)/2.0) ) ;
74  }else{
75  return par[0]+par[1]+par[2]+par[3]+par[4]+par[5]+par[6]+par[7]+par[8]+par[9];
76  //return 0;
77  }
78 }
79 //_______________________________________________________________________________________
80 //
81 Double_t SumGausLorentz2D(Double_t *x, Double_t *par){
82  //par[0] Norm_1
83  //par[1] x0_1
84  //par[2] sigma_x_1
85  //par[3] y0_1
86  //par[4] sigma_y_1
87  // par[5] = Gamma
88  // par[6] = x0
89  // par[7] = y0
90  // par[8] = norm
91  //0 0.146571+-0 0+-0 7.55129+-0 0+-0 5.94133+-0 3.2357+-0 0+-0 0+-0 27.3458+-0 chi2 / NDF = 3126.16 / 10009 = 0.312335
92  //1 0.1459+-0 0+-0 14.7913+-0 0+-0 11.7393+-0 6.86661+-0 0+-0 0+-0 58.7235+-0 chi2 / NDF = 7458.34 / 10009 = 0.745163
93  //2 0.149045+-0 0+-0 22.2327+-0 0+-0 17.6162+-0 9.20632+-0 0+-0 0+-0 90.0768+-0 chi2 / NDF = 12700.5 / 10009 = 1.26891
94  //3 0.168693+-0 0+-0 27.9107+-0 0+-0 21.8062+-0 12.9946+-0 0+-0 0+-0 102.019+-0 chi2 / NDF = 8275.98 / 10009 = 0.826854
95  Double_t Pi = 3.141592653;
96  if((par[2]!=0.0) && (par[4]!=0.0) ){
97  double rx = (x[0]-par[1])/par[2];
98  double ry = (x[1]-par[3])/par[4];
99  return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)
100  +par[8]*par[5]/Pi/((x[0]-par[6])*(x[0]-par[6])+(x[1]-par[7])*(x[1]-par[7])+par[5]*par[5])
101  );
102  }else{
103  return 0;
104  }
105 }
106 //_______________________________________________________________________________________
107 //
108 Double_t SumGausLorentz2Dnew(Double_t *x, Double_t *par){
109  //par[0] Norm_1
110  //par[1] x0_1
111  //par[2] sigma_x_1
112  //par[3] y0_1
113  //par[4] sigma_y_1
114  // par[5] = Gamma
115  // par[6] = x0
116  // par[7] = y0
117  // par[8] = norm
118  // par[9] = normgaus2
119  // par[10] = sigma_2
120 
121  /* for (int j=0;j<10; j++) {
122  cout<<"par " <<j<<" "<<par[j]<<endl;
123  }
124  */
125  Double_t Pi = 3.141592653;
126  if((par[2]!=0.0) && (par[4]!=0.0) ){
127  double rx = (x[0]-par[1])/par[2];
128  double ry = (x[1]-par[3])/par[4];
129 
130  double rx2 = (x[0])/par[10];
131  double ry2 = (x[1])/par[10];
132 
133 
134  if( TMath::Sqrt(x[0]*x[0]+x[1]*x[1])<200.0 ){
135  return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)
136  +par[8]*par[5]/Pi/((x[0]-par[6])*(x[0]-par[6])+(x[1]-par[7])*(x[1]-par[7])+par[5]*par[5])
137  +par[9]*TMath::Exp(-(rx2*rx2+ry2*ry2)/2.0)
138  );
139  }else{
140  return par[0]*( TMath::Exp(-(rx*rx+ry*ry)/2.0)
141  +par[8]*par[5]/Pi/((x[0]-par[6])*(x[0]-par[6])+(x[1]-par[7])*(x[1]-par[7])+par[5]*par[5])
142 
143  );
144  }
145  }else{
146  return 0;
147 
148  }
149 }
150 //_______________________________________________________________________________________
151 //
152 Double_t SumGausLorentz1D(Double_t *x, Double_t *par){
153 
154  //par[0] Norm_g
155  //par[1] x0_g
156  //par[2] sigma_g
157  // par[3] = Gamma_lor
158  // par[4] = x0_lor
159  // par[5] = norm_lor
160 
161 
162 
163  /* for (int j=0;j<10; j++) {
164  cout<<"par " <<j<<" "<<par[j]<<endl;
165  }
166  */
167  Double_t Pi = 3.141592653;
168  if((par[2]!=0.0) ){
169  Double_t rx = (x[0]-par[1])/par[2];
170  Double_t tempoutput;
171  tempoutput= par[0]*( TMath::Exp(-(rx*rx)/2.0)
172  +par[5]*par[3]/Pi/ ((x[0]-par[4])*(x[0]-par[4]) +par[3]*par[3])
173  );
174  // cout<<"SumGausLorentz1D " <<tempoutput<<endl;
175  return tempoutput;
176 
177  }else{
178  return 0;
179  }
180 }
181 //==============================================================================
182 ClassImp(DIGParticle)
183 DIGParticle::DIGParticle()
184 {
185  //
186  // default constructor
187  //
188 
189 }
190 //______________________________________________________________________________
191 //
192 DIGParticle::DIGParticle(Float_t EntryX, Float_t EntryY, Float_t EntryZ,
193  Float_t ExitX, Float_t ExitY, Float_t ExitZ, Float_t Energy_deposited)
194 {
195  fEntryX = EntryX;
196  fEntryY = EntryY;
197  fEntryZ = EntryZ;
198  fExitX = ExitX ;
199  fExitY = ExitY ;
200  fExitZ = ExitZ ;
201  fEnergy_deposited = Energy_deposited;
202  fNSegment =0;
203  fSegmentX.clear();
204  fSegmentY.clear();
205  fSegmentZ.clear();
206  fSegmentCharge.clear();
207 
208  fNpixels=0;
209  fPixelMap.clear();
210  fAnalogChargeMap.clear();
211  fDigitalChargeMap.clear();//not used ?
212 
213 }
214 //______________________________________________________________________________
215 //
216 DIGParticle::DIGParticle(DIGParticle & adigparticle) : TObject()
217 {
218  fEntryX = adigparticle.GetEntryX();
219  fEntryY = adigparticle.GetEntryY();
220  fEntryZ = adigparticle.GetEntryZ();
221  fExitX = adigparticle.GetExitX();
222  fExitY = adigparticle.GetExitY();
223  fExitZ = adigparticle.GetExitZ();
224  fEnergy_deposited = adigparticle.GetEnergy_deposited();
225  fNSegment =adigparticle.GetNSegment();
226  fSegmentX.resize(fNSegment);
227  fSegmentY.resize(fNSegment);
228  fSegmentZ.resize(fNSegment);
229  fSegmentCharge.resize(fNSegment);
230  // a modifier:
231  for (Int_t i=0 ; i<fNSegment ; i++){
232  fSegmentX[i] = adigparticle.GetSegmentX()[i];
233  fSegmentY[i] = adigparticle.GetSegmentY()[i];
234  fSegmentZ[i] = adigparticle.GetSegmentZ()[i];
235  fSegmentCharge[i] = adigparticle.GetSegmentCharge()[i];
236  }
237  fNpixels = adigparticle.GetNpixels();
238  fPixelMap.resize(fNpixels);
239  fAnalogChargeMap.resize(fNpixels);
240  fDigitalChargeMap.resize(fNpixels);
241  for (Int_t i=0 ; i<fNpixels ; i++){
242  fPixelMap[i] = adigparticle.GetPixelMap()[i];
243  fAnalogChargeMap[i] = adigparticle.GetAnalogCharge()[i];
244  fDigitalChargeMap[i] = adigparticle.GetDigitalCharge()[i];
245  }
246 }
247 //______________________________________________________________________________
248 //
249 
250 
251 DIGParticle::~DIGParticle() { //
252  // virtual destructor
253  //
254  // delete fLayers;
255 }
256 //______________________________________________________________________________
257 //
258 void DIGParticle::Clear(const Option_t *)
259 {
260 
261 }
262 //______________________________________________________________________________
263 //
264 void DIGParticle::ComputeChargeDeposition(Float_t StartingSegmentSize, Float_t MaximumSegmentSize,
265  Float_t MaximumChargePerSegment)
266 {
267  Float_t SegmentSize = StartingSegmentSize;
268  Float_t TotalLength = TMath::Sqrt((fExitX-fEntryX)*(fExitX-fEntryX)
269  +(fExitY-fEntryY)*(fExitY-fEntryY)
270  +(fExitZ-fEntryZ)*(fExitZ-fEntryZ));
271 
272  Float_t ChargePerSegment = 0.0;
273  if(SegmentSize<0.0){
274  SegmentSize=0.001;
275  }
276  fNSegment = int(TotalLength*1.000001/SegmentSize) ;
277  if(fNSegment<1){
278  fNSegment=1;
279  }
280  SegmentSize = TotalLength/float(fNSegment);
281  ChargePerSegment = fEnergy_deposited / float(fNSegment);
282 
283  while((SegmentSize>MaximumSegmentSize)&&(ChargePerSegment > MaximumChargePerSegment)){
284  Int_t newNSegment = int(fNSegment *1.1);
285  if(newNSegment==fNSegment){fNSegment++;}
286  SegmentSize = TotalLength/float(fNSegment);
287  ChargePerSegment = fEnergy_deposited / float(fNSegment);
288  }
289  Float_t xstep = fExitX-fEntryX;
290  Float_t ystep = fExitY-fEntryY;
291  Float_t zstep = fExitZ-fEntryZ;
292 
293 
294  for (Int_t i=0 ; i<fNSegment ; i++){
295  fSegmentX.push_back(fEntryX + (float(i+0.5)* xstep/float(fNSegment)) );
296  fSegmentY.push_back(fEntryY + (float(i+0.5)* ystep/float(fNSegment)) );
297  fSegmentZ.push_back(fEntryZ + (float(i+0.5)* zstep/float(fNSegment)) );
298  fSegmentCharge.push_back(ChargePerSegment );
299  }
300 
301 
302 }
303 //______________________________________________________________________________
304 //
305 
306 void DIGParticle::PrintInfo() {
307  std::cout<<"---------Particle properties------------- "<<endl;
308  std::cout<<"fEntryX fEntryY fEntryZ"<<endl;
309  std::cout<<fEntryX<<" "<< fEntryY<<" "<<fEntryZ<<endl;
310  std::cout<<"fExitX fExitY fExitZ fEnergy_deposited"<<endl;
311  std::cout<<fExitX <<" "<<fExitY <<" "<<fExitZ <<" "<<fEnergy_deposited <<endl;
312  std::cout<<fNSegment<<" Segments X Y Z Charge"<<endl;
313  /* for (Int_t i=0 ; i<fNSegment ; i++){
314  std::cout<<i<<" "<< fSegmentX[i]<<" "<<fSegmentY[i] <<" "<<fSegmentZ[i] <<" "<< fSegmentCharge[i]<<endl;
315  }*/
316  cout<<" size vectors "<< fPixelMap.size()<<" "<<fAnalogChargeMap.size()<<" "<<fDigitalChargeMap.size()<<endl;
317  std::cout<<fNpixels<<" fNpixels map analog digital "<<endl;
318  for (Int_t i=0 ; i<fNpixels ; i++){
319  std::cout<<i<<" "<<fPixelMap[i]<<" "<<fAnalogChargeMap[i]<<" "<<fDigitalChargeMap[i]<<endl;
320  }
321 }
322 //______________________________________________________________________________
323 //
324 void DIGParticle::ComputeChargeTransport(DIGPlane *aDIGPlane,DIGTransport *aDIGTransport){
325 
326  Float_t pitchX = aDIGPlane->GetPitchX();
327  Float_t pitchY = aDIGPlane->GetPitchY();
328 
329  Float_t lorentz2Dmodel_Cp0 = aDIGTransport->GetLorentz2DModel_Cp0();
330  Float_t lorentz2Dmodel_Cp1 = aDIGTransport->GetLorentz2DModel_Cp1();
331  Float_t rangelimit_inpitchunit = aDIGTransport->GetRangeLimit_InPitchUnit();
332 
333  Int_t ChargeModel = aDIGTransport->GetChargeModel();
334 
335  Float_t Gauss2DModel_sigma1_Cp0 = aDIGTransport->GetGauss2DModel_sigma1_Cp0();
336  Float_t Gauss2DModel_sigma1_Cp1 = aDIGTransport->GetGauss2DModel_sigma1_Cp1();
337  Float_t Gauss2DModel_sigma2_Cp0 = aDIGTransport->GetGauss2DModel_sigma2_Cp0();
338  Float_t Gauss2DModel_sigma2_Cp1 = aDIGTransport->GetGauss2DModel_sigma2_Cp1();
339  Float_t Gauss2DModel_weight = aDIGTransport->GetGauss2DModel_weight();
340 
341  Float_t LorGaussModel_Norm1_Cp0 = aDIGTransport->GetLorGaussModel_Norm1_Cp0();
342  Float_t LorGaussModel_Norm1_Cp1 = aDIGTransport->GetLorGaussModel_Norm1_Cp1();
343  Float_t LorGaussModel_Norm1_Cp2 = aDIGTransport->GetLorGaussModel_Norm1_Cp2();
344  Float_t LorGaussModel_sigma_Cp0 = aDIGTransport->GetLorGaussModel_sigma_Cp0();
345  Float_t LorGaussModel_sigma_Cp1 = aDIGTransport->GetLorGaussModel_sigma_Cp1();
346  Float_t LorGaussModel_C_Cp0 = aDIGTransport->GetLorGaussModel_C_Cp0();
347  Float_t LorGaussModel_C_Cp1 = aDIGTransport->GetLorGaussModel_C_Cp1();
348  Float_t LorGaussModel_Norm_Cp0 = aDIGTransport->GetLorGaussModel_Norm_Cp0();
349  Float_t LorGaussModel_Norm_Cp1 = aDIGTransport->GetLorGaussModel_Norm_Cp1();
350 
351  Float_t lorlorgausModel_Norm1_Cp0 = aDIGTransport->GetlorlorgausModel_Norm1_Cp0() ;
352  Float_t lorlorgausModel_Norm1_Cp1 = aDIGTransport->GetlorlorgausModel_Norm1_Cp1() ;
353  Float_t lorlorgausModel_x01_Cp0 = aDIGTransport->GetlorlorgausModel_x01_Cp0() ;
354  Float_t lorlorgausModel_x01_Cp1 = aDIGTransport->GetlorlorgausModel_x01_Cp1() ;
355  Float_t lorlorgausModel_sigmax1_Cp0 = aDIGTransport->GetlorlorgausModel_sigmax1_Cp0();
356  Float_t lorlorgausModel_sigmax1_Cp1 = aDIGTransport->GetlorlorgausModel_sigmax1_Cp1();
357  Float_t lorlorgausModel_y01_Cp0 = aDIGTransport->GetlorlorgausModel_y01_Cp0();
358  Float_t lorlorgausModel_y01_Cp1 = aDIGTransport->GetlorlorgausModel_y01_Cp1() ;
359  Float_t lorlorgausModel_sigmay1_Cp0 = aDIGTransport->GetlorlorgausModel_sigmay1_Cp0();
360  Float_t lorlorgausModel_sigmay1_Cp1 = aDIGTransport->GetlorlorgausModel_sigmay1_Cp1();
361  Float_t lorlorgausModel_Gamma_Cp0 = aDIGTransport->GetlorlorgausModel_Gamma_Cp0() ;
362  Float_t lorlorgausModel_Gamma_Cp1 = aDIGTransport->GetlorlorgausModel_Gamma_Cp1() ;
363  Float_t lorlorgausModel_x0_Cp0 = aDIGTransport->GetlorlorgausModel_x0_Cp0() ;
364  Float_t lorlorgausModel_x0_Cp1 = aDIGTransport->GetlorlorgausModel_x0_Cp1() ;
365  Float_t lorlorgausModel_y0_Cp0 = aDIGTransport->GetlorlorgausModel_y0_Cp0() ;
366  Float_t lorlorgausModel_y0_Cp1 = aDIGTransport->GetlorlorgausModel_y0_Cp1() ;
367  Float_t lorlorgausModel_norm_Cp0 = aDIGTransport->GetlorlorgausModel_norm_Cp0() ;
368  Float_t lorlorgausModel_norm_Cp1 = aDIGTransport->GetlorlorgausModel_norm_Cp1() ;
369  Float_t lorlorgausModel_normgaus2_Cp0 = aDIGTransport->GetlorlorgausModel_normgaus2_Cp0();
370  Float_t lorlorgausModel_normgaus2_Cp1 = aDIGTransport->GetlorlorgausModel_normgaus2_Cp1() ;
371  Float_t lorlorgausModel_sigma2_Cp0 = aDIGTransport->GetlorlorgausModel_sigma2_Cp0() ;
372  Float_t lorlorgausModel_sigma2_Cp1 = aDIGTransport->GetlorlorgausModel_sigma2_Cp1() ;
373 
374  Float_t l1dimgauslor_Norm_g_1st = aDIGTransport->Getf1dimgauslor_Norm_g_1st();
375  Float_t l1dimgauslor_x0_g_1st = aDIGTransport->Getf1dimgauslor_x0_g_1st();
376  Float_t l1dimgauslor_sigma_g_1st = aDIGTransport->Getf1dimgauslor_sigma_g_1st();
377  Float_t l1dimgauslor_Gamma_lor_1st = aDIGTransport->Getf1dimgauslor_Gamma_lor_1st();
378  Float_t l1dimgauslor_x0_lor_1st = aDIGTransport->Getf1dimgauslor_x0_lor_1st();
379  Float_t l1dimgauslor_norm_lor_1st = aDIGTransport->Getf1dimgauslor_norm_lor_1st();
380  Float_t l1dimgauslor_Norm_g_2nd = aDIGTransport->Getf1dimgauslor_Norm_g_2nd();
381  Float_t l1dimgauslor_x0_g_2nd = aDIGTransport->Getf1dimgauslor_x0_g_2nd();
382  Float_t l1dimgauslor_sigma_g_2nd = aDIGTransport->Getf1dimgauslor_sigma_g_2nd();
383  Float_t l1dimgauslor_Gamma_lor_2nd = aDIGTransport->Getf1dimgauslor_Gamma_lor_2nd();
384  Float_t l1dimgauslor_x0_lor_2nd = aDIGTransport->Getf1dimgauslor_x0_lor_2nd();
385  Float_t l1dimgauslor_norm_lor_2nd = aDIGTransport->Getf1dimgauslor_norm_lor_2nd();
386 
387 
388  GlobalSeed++;
389  TRandom3 *r3 = new TRandom3(GlobalSeed);
390  r3->SetSeed(GlobalSeed);
391 
392  TF2 *mymodel2D=0;
393  TF1 *mymodel1D_1st=0;
394  TF1 *mymodel1D_2nd=0;
395  //-----------------------------------------------------
396  // chose model
397  //-----------------------------------------------------
398  if(ChargeModel==1){
399  //-----------------------------------------------------
400  // model Lorentz
401  //-----------------------------------------------------
402 
403  Double_t Cvalue= lorentz2Dmodel_Cp0 + pitchX * lorentz2Dmodel_Cp1 ;
404 
405  mymodel2D = new TF2("funlor2d",Lorentz2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
406  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,4);
407  mymodel2D->SetParNames("C","X_{0}","Y_{0}","Norm");
408  Double_t params1[] = {Cvalue,0.,0.,1.0};
409  mymodel2D->SetParameters(params1);
410 
411  }else if(ChargeModel==2){
412  //-----------------------------------------------------
413  // 2xGaussian model
414  //-----------------------------------------------------
415  //Double_t SumGaus2D(Double_t *x, Double_t *par){
416  //par[0] Norm_1
417  //par[1] x0_1
418  //par[2] sigma_x_1
419  //par[3] y0_1
420  //par[4] sigma_y_1
421  //par[5] Norm_2
422  //par[6] x0_2
423  //par[7] sigma_x_2
424  //par[8] y0_2
425  //par[9] sigma_y_2
426  Double_t Gsigma1 = Gauss2DModel_sigma1_Cp0 + pitchX * Gauss2DModel_sigma1_Cp1 ;
427  Double_t Gsigma2 = Gauss2DModel_sigma2_Cp0 + pitchX * Gauss2DModel_sigma2_Cp1 ;
428 
429  mymodel2D = new TF2("funlor2d",SumGaus2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
430  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,10);
431  mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Norm_2","x0_2","sigma_x_2","y0_2","sigma_y_2");
432  Double_t params1[] = {1.0,0.0,Gsigma1,0.0,Gsigma1,Gauss2DModel_weight,0.0,Gsigma2,0.0,Gsigma2,0.0};
433  mymodel2D->SetParameters(params1);
434  }else if(ChargeModel==3){
435  //-----------------------------------------------------
436  // Lorentz + Gauss model
437  //-----------------------------------------------------
438  //par[0] Norm_1
439  //par[1] x0_1
440  //par[2] sigma_x_1
441  //par[3] y0_1
442  //par[4] sigma_y_1
443  // par[5] = Gamma
444  // par[6] = x0
445  // par[7] = y0
446  // par[8] = norm
447  Double_t Norm1Value = LorGaussModel_Norm1_Cp0 + pitchX * LorGaussModel_Norm1_Cp1 + pitchX *pitchX * LorGaussModel_Norm1_Cp2;
448  Double_t sigmaValue = LorGaussModel_sigma_Cp0 + pitchX *LorGaussModel_sigma_Cp1;
449  Double_t Cvalue= LorGaussModel_C_Cp0 + pitchX *LorGaussModel_C_Cp1;
450  Double_t NormValue =LorGaussModel_Norm_Cp0 + pitchX * LorGaussModel_Norm_Cp1;
451  mymodel2D = new TF2("funlor2d",SumGausLorentz2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
452  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,9);
453  mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Gamma","x0","y0","norm");
454  Double_t params1[] = {Norm1Value,0.0,sigmaValue,0.0,sigmaValue,Cvalue,0.0,0.0,NormValue};
455  mymodel2D->SetParameters(params1);
456  }else if(ChargeModel==4){
457  //-----------------------------------------------------
458  // Lorentz + Gauss + Gauss model
459  //-----------------------------------------------------
460  //par[0] Norm_1
461  //par[1] x0_1
462  //par[2] sigma_x_1
463  //par[3] y0_1
464  //par[4] sigma_y_1
465  // par[5] = Gamma
466  // par[6] = x0
467  // par[7] = y0
468  // par[8] = norm
469  // par[9] = normgaus2
470  // par[10] = sigma_2
471  /*
472  SumGausLorentz2Dnew
473  */
474  Double_t Norm1Value = lorlorgausModel_Norm1_Cp0 + pitchX * lorlorgausModel_Norm1_Cp1;
475  Double_t x0_1Value = lorlorgausModel_x01_Cp0 + pitchX * lorlorgausModel_x01_Cp1;
476  Double_t sigma_x_1Value = lorlorgausModel_sigmax1_Cp0 + pitchX * lorlorgausModel_sigmax1_Cp1;
477  Double_t y0_1Value = lorlorgausModel_y01_Cp0 + pitchX * lorlorgausModel_y01_Cp1;
478  Double_t sigma_y_1Value = lorlorgausModel_sigmay1_Cp0 + pitchX * lorlorgausModel_sigmay1_Cp1;
479  Double_t GammaValue = lorlorgausModel_Gamma_Cp0 + pitchX * lorlorgausModel_Gamma_Cp1;
480  Double_t x0Value = lorlorgausModel_x0_Cp0 + pitchX * lorlorgausModel_x0_Cp1;
481  Double_t y0Value = lorlorgausModel_y0_Cp0 + pitchX * lorlorgausModel_y0_Cp1;
482  Double_t normValue = lorlorgausModel_norm_Cp0 + pitchX * lorlorgausModel_norm_Cp1;
483  Double_t normgaus2Value = lorlorgausModel_normgaus2_Cp0 + pitchX * lorlorgausModel_normgaus2_Cp1;
484  Double_t sigma_2Value = lorlorgausModel_sigma2_Cp0 + pitchX * lorlorgausModel_sigma2_Cp1;
485  mymodel2D = new TF2("funlor2d",SumGausLorentz2Dnew,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
486  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,11);
487  mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Gamma","x0","y0","norm","normgaus2","sigma_2");
488  Double_t params1[] = {Norm1Value,x0_1Value,sigma_x_1Value,y0_1Value,
489  sigma_y_1Value,GammaValue,x0Value,y0Value,normValue,normgaus2Value,sigma_2Value};
490  mymodel2D->SetParameters(params1);
491  }else if(ChargeModel==5){
492 
493  //-----------------------------------------------------
494  // 1 dimension model with 2 function for each squares around the hit
495  //-----------------------------------------------------
496  // "Norm_g","x0_g","sigma_g","Gamma_lor","x0_lor","norm_lor"
497  Double_t Norm_g_1st =l1dimgauslor_Norm_g_1st ;
498  Double_t x0_g_1st =l1dimgauslor_x0_g_1st ;
499  Double_t sigma_g_1st =l1dimgauslor_sigma_g_1st ;
500  Double_t Gamma_lor_1st =l1dimgauslor_Gamma_lor_1st ;
501  Double_t x0_lor_1st =l1dimgauslor_x0_lor_1st ;
502  Double_t norm_lor_1st =l1dimgauslor_norm_lor_1st ;
503  Double_t Norm_g_2nd =l1dimgauslor_Norm_g_2nd ;
504  Double_t x0_g_2nd =l1dimgauslor_x0_g_2nd ;
505  Double_t sigma_g_2nd =l1dimgauslor_sigma_g_2nd ;
506  Double_t Gamma_lor_2nd =l1dimgauslor_Gamma_lor_2nd ;
507  Double_t x0_lor_2nd =l1dimgauslor_x0_lor_2nd ;
508  Double_t norm_lor_2nd =l1dimgauslor_norm_lor_2nd ;
509 
510 
511  mymodel1D_1st = new TF1("namemodel1D_1st",SumGausLorentz1D,0,4.0*pitchX,6);
512  mymodel1D_2nd = new TF1("namemodel1D_2nd",SumGausLorentz1D,0,4.0*pitchX,6);
513  mymodel1D_1st->SetParNames("Norm_g","x0_g","sigma_g","Gamma_lor","x0_lor","norm_lor");
514  mymodel1D_2nd->SetParNames("Norm_g","x0_g","sigma_g","Gamma_lor","x0_lor","norm_lor");
515  Double_t params1[] = {Norm_g_1st,x0_g_1st,sigma_g_1st, Gamma_lor_1st,x0_lor_1st,norm_lor_1st};
516  Double_t params2[] = {Norm_g_2nd,x0_g_2nd,sigma_g_2nd, Gamma_lor_2nd,x0_lor_2nd,norm_lor_2nd};
517  mymodel1D_1st->SetParameters(params1);
518  mymodel1D_2nd->SetParameters(params2);
519 
520  //std::cout<<" LIST OF PARAMETERS 1st "<<Norm_g_1st<<" "<<x0_g_1st<<" "<<sigma_g_1st<<" "<<Gamma_lor_1st<<" "<<x0_lor_1st<<" "<<norm_lor_1st<<endl;
521  //std::cout<<" LIST OF PARAMETERS 2nd "<<Norm_g_2nd<<" "<<x0_g_2nd<<" "<<sigma_g_2nd<<" "<<Gamma_lor_2nd<<" "<<x0_lor_2nd<<" "<<norm_lor_2nd<<endl;
522  //std::cout<<" function proba "<<mymodel1D_1st->Eval(0)<<" "<<mymodel1D_1st->Eval(15.0)<<mymodel1D_1st->Eval(38.)<<mymodel1D_1st->Eval(70.)<<endl;
523  //std::cout<<" function proba "<<mymodel1D_2nd->Eval(0)<<" "<<mymodel1D_2nd->Eval(15.0)<<mymodel1D_2nd->Eval(38.)<<mymodel1D_2nd->Eval(70.)<<endl;
524 
525  }else{
526  Double_t Gsigma1 = Gauss2DModel_sigma1_Cp0 + pitchX * Gauss2DModel_sigma1_Cp1 ;
527  Double_t Gsigma2 = Gauss2DModel_sigma2_Cp0 + pitchX * Gauss2DModel_sigma2_Cp1 ;
528 
529  mymodel2D = new TF2("funlor2d",SumGaus2D,-rangelimit_inpitchunit*pitchX,rangelimit_inpitchunit*pitchX,
530  -rangelimit_inpitchunit*pitchY,rangelimit_inpitchunit*pitchY,10);
531  mymodel2D->SetParNames("Norm_1","x0_1","sigma_x_1","y0_1","sigma_y_1","Norm_2","x0_2","sigma_x_2","y0_2","sigma_y_2");
532  Double_t params1[] = {1.0,0.0,Gsigma1,0.0,Gsigma1,Gauss2DModel_weight,0.0,Gsigma2,0.0,Gsigma2,0.0};
533  mymodel2D->SetParameters(params1);
534 
535  }
536 
537 
538  //---------------------------------------------------------------------
539  // method A: lorentz function define a probability density function which gives where the charge is collected.
540  //---------------------------------------------------------------------
541  /*
542  Float_t Xdimension = aDIGPlane->GetXdimension();
543  Float_t Ydimension = aDIGPlane->GetYdimension();
544 
545  for (Int_t i=0 ; i<fNSegment ; i++){
546  //std::cout<<" i="<<i<<" x="<< fSegmentX[i]<<" y="<<fSegmentY[i] <<" z="<<fSegmentZ[i] <<" c="<< fSegmentCharge[i]<<endl;
547  Double_t xrandom;
548  Double_t yrandom;
549  gRandom = new TRandom3;
550  GlobalSeed++;
551  gRandom->SetSeed(GlobalSeed);
552  mylorentz2D->GetRandom2(xrandom,yrandom);
553 
554  Double_t xpos_aftertransport = fSegmentX[i] + xrandom;
555  Double_t ypos_aftertransport = fSegmentY[i] + yrandom;
556  //if charge reaches a x,y position outside the size of the detector, consider that the charge is lost.
557  if((xpos_aftertransport>=0.0)
558  &&(xpos_aftertransport<=Xdimension)
559  &&(ypos_aftertransport>=0.0)
560  &&(ypos_aftertransport<=Ydimension)){
561  Int_t PixelReached = GetPixelNumber(aDIGPlane, xpos_aftertransport , ypos_aftertransport );
562  UpdatePixel(fSegmentCharge[i],PixelReached);
563  }else{
564  // Int_t PixelReached = GetPixelNumber(aDIGPlane, xpos_aftertransport , ypos_aftertransport );
565  // Int_t xnpos, ynpos;
566  // GetXYPixelNumber(xnpos,ynpos, aDIGPlane, PixelReached);
567  // Float_t xdpos, ydpos;
568  // GetXYPixelCenter(xdpos,ydpos, aDIGPlane, PixelReached);
569  // std::cout<<"charge outside DIMENSIONS "<<Xdimension <<" x "<<Ydimension <<endl;
570  // std::cout<<" after ranx="<< xrandom<<" rany="<<yrandom<<" x="<<xpos_aftertransport<<" y="<<ypos_aftertransport
571  // <<" n="<<PixelReached<<" xn= "<<xnpos
572  // <<" yn="<<ynpos <<" xd="<<xdpos<<" yd="<<ydpos<<endl;
573 
574  }
575  }
576  */
577  //---------------------------------------------------------------------
578  // method B: more realistic model. The lorentz function allows to compute 25 probabilities (5x5 cluster)
579  // the random number is generated to see where the segment charge is deposited.
580  //---------------------------------------------------------------------
581 
582  for (Int_t i=0 ; i<fNSegment ; i++){
583 
584  Int_t PixelReached = GetPixelNumber(aDIGPlane, fSegmentX[i], fSegmentY[i]);
585  Float_t xdpos, ydpos;
586  GetXYPixelCenter(xdpos, ydpos, aDIGPlane, PixelReached);
587  Int_t xpixn,ypixn;
588  GetXYPixelNumber( xpixn,ypixn, aDIGPlane, PixelReached);
589 
590  Float_t TotalProba=0.0;
591  const Int_t Npix=25;
592  Float_t PixelposX[Npix];
593  Float_t PixelposY[Npix];
594  Float_t Pixelproba[Npix];
595  Int_t Pixelnumber[Npix];
596  for (Int_t j=0 ; j<Npix ; j++){
597  PixelposX[j]=0.0;
598  PixelposY[j]=0.0;
599  Pixelproba[j]=0.0;
600  Pixelnumber[j]=0;
601  }
602 
603  for (Int_t j=0 ; j<Npix ; j++){
604  Int_t xi = j%5;
605  Int_t yi = j/5;
606  Float_t xeval = -(fSegmentX[i]-(xpixn+0.5)*pitchX) + (float(xi)-2.0)*pitchX ;
607  Float_t yeval = -(fSegmentY[i]-(ypixn+0.5)*pitchY) + (float(yi)-2.0)*pitchY ;
608  Float_t xydist = TMath::Sqrt(xeval*xeval+yeval*yeval);
609  // force Lorentz function to be centered in the seed pixel for tests:
610  //Float_t xeval = (float(xi)-2.0)*pitchX ;
611  //Float_t yeval = (float(yi)-2.0)*pitchY ;
612 
613  if((ChargeModel==1)||(ChargeModel==2)||(ChargeModel==3)||(ChargeModel==4)){
614  Pixelproba[j]=mymodel2D->Eval(xeval,yeval);
615  }else if(ChargeModel==5){
616  //find if pixel is one of the four in the first square around the hit:
617  if((fabs(xeval)<=pitchX)&&(fabs(yeval)<=pitchY)){
618  Pixelproba[j]=mymodel1D_1st->Eval(xydist);
619  }else{
620  Pixelproba[j]=mymodel1D_2nd->Eval(xydist);
621  }
622  }
623  TotalProba+=Pixelproba[j];
624  PixelposX[j]= xdpos + (float(xi)-2.0)*pitchX;
625  PixelposY[j]= ydpos + (float(yi)-2.0)*pitchY;
626  if( (PixelposX[j]<=0.0)||
627  (PixelposX[j]>=aDIGPlane->GetXdimension())||
628  (PixelposY[j]<=0.0)||
629  (PixelposY[j]>=aDIGPlane->GetYdimension())){
630  Pixelnumber[j] = -1;
631  }else{
632  Pixelnumber[j] = GetPixelNumber(aDIGPlane, PixelposX[j], PixelposY[j]);
633  }
634  //std::cout<<"j ChargeModel PIXEL PROBA / TOT / xydist " <<j<<" "<<ChargeModel<<" "<<Pixelproba[j]<<" / "<<TotalProba<<" "<<xydist<<endl;
635  // <<mymodel1D_1st->Eval(xydist)<<" "<< mymodel1D_2nd->Eval(xydist) <<endl;
636  }
637 
638  //for (Int_t j=0 ; j<Npix ; j++){
639  // std::cout<<" PIXEL PROBA " <<j<<" "<<Pixelproba[j]<<" / "<<TotalProba<<endl;
640  // }
641 
642  Double_t rando = r3->Rndm()*TotalProba;
643  Bool_t reached = false;
644  Double_t incrprob = 0.0;
645  Int_t icounter=0;
646  while((!reached)&&(icounter<Npix)){
647  incrprob+=Pixelproba[icounter];
648  if(incrprob>rando){
649  reached=true;
650  }else{
651  icounter++;
652  }
653  }
654 
655  if(Pixelnumber[icounter]>=0){
656  UpdatePixel(fSegmentCharge[i],Pixelnumber[icounter]);
657  }
658  }
659  delete r3;
660 }
661 //______________________________________________________________________________
662 //
663 
664 void DIGParticle::AddPixel(Float_t AnalogCharge, Int_t PixelNumber){
665  fNpixels++;
666  fPixelMap.push_back(PixelNumber);
667  fAnalogChargeMap.push_back(AnalogCharge);
668  fDigitalChargeMap.push_back(0);
669 }
670 
671 //______________________________________________________________________________
672 //
673 void DIGParticle::UpdatePixel(Float_t AnalogCharge, Int_t PixelNumber){
674  Bool_t found = false;
675  Int_t i=0;
676  while((!found)&&(i<fNpixels)){
677  if(PixelNumber==fPixelMap[i]){
678  found=true;
679  fAnalogChargeMap[i]+=AnalogCharge;
680  }
681  i++;
682  }
683  if(!found){
684  AddPixel(AnalogCharge, PixelNumber);
685  }
686 
687 }
688 //______________________________________________________________________________
689 //
690 void DIGParticle::AddRandomNoise(DIGPlane *myDIGPlane){
691  Double_t Noisefix = myDIGPlane->GetNoiseElectrons();
692  Double_t Noise;
693  for (Int_t i = 0; i < fNpixels ; i++){
694  Noise = GaussianLaw(0.0, Noisefix);
695  fAnalogChargeMap[i]+=Noise;
696  }
697 }
698 //______________________________________________________________________________
699 //
700 void DIGParticle::AnalogToDigitalconversion(DIGADC *myDIGADC, DIGPlane *myDIGPlane ){
701  Float_t Noisefix = myDIGPlane->GetNoiseElectrons();
702  if(Noisefix<=0.0){
703  std::cout <<"<---- DIGParticle::AnalogToDigitalconversion --->"<<endl;
704  std::cout<<"WARNING negative or null Noise is not physical, please correct the input file"<<endl;
705  Noisefix = 1.0;
706  }
707  Float_t *myADCthreshold = 0;
708  myADCthreshold = myDIGADC->GetADC_thresholds();
709  for (Int_t i = 0; i < fNpixels ; i++){
710  if (fAnalogChargeMap[i]<=0.0){
711  fDigitalChargeMap[i]=0;
712  }else{
713  Bool_t thresholdfound = false;
714  Int_t ithres = 0;
715  fDigitalChargeMap[i]=0;
716  while((thresholdfound==false)&&(ithres< (myDIGADC->GetNThresholds()) )){
717  if( (fAnalogChargeMap[i]/Noisefix) < myADCthreshold[ithres] ){
718  thresholdfound = true;
719  }else{
720  fDigitalChargeMap[i]++;
721  ithres++;
722  }
723  }
724  }
725  }
726 }
727 //______________________________________________________________________________
728 //
729 Int_t DIGParticle::GetPixelNumber(DIGPlane *myDIGPlane, Float_t Xpos, Float_t Ypos){
730  Int_t XPixelnumber = int(Xpos/(myDIGPlane->GetPitchX()));
731  Int_t YPixelnumber = int(Ypos/(myDIGPlane->GetPitchY()));
732  Int_t PixelNumber = XPixelnumber + (myDIGPlane->GetNpixelsX())*YPixelnumber;
733 
734  if(
735  (Xpos<0.0)||
736  (Xpos>myDIGPlane->GetXdimension())||
737  (Ypos<0.0)||
738  (Ypos>myDIGPlane->GetYdimension())){
739 // cout<<" DIGParticle::GetPixelNumber WARNING charge is going outside the plane limits: x/y([0,"<<myDIGPlane->GetXdimension()<<"],[0,"<<myDIGPlane->GetYdimension()<<"])="<<Xpos<<"/"<<Ypos<<endl;
740  return 0;
741  }else{
742  return PixelNumber;
743  }
744 }
745 //______________________________________________________________________________
746 //
747 void DIGParticle::GetXYPixelNumber(Int_t &Xpix, Int_t &Ypix, DIGPlane *myDIGPlane, Int_t PixelNumber){
748 
749  Xpix = PixelNumber%(myDIGPlane->GetNpixelsX());
750  Ypix = PixelNumber/(myDIGPlane->GetNpixelsX());
751 
752 }
753 //______________________________________________________________________________
754 //
755 void DIGParticle::GetXYPixelCenter(Float_t &Xpix, Float_t &Ypix, DIGPlane *myDIGPlane, Int_t PixelNumber){
756 
757  Xpix = (myDIGPlane->GetPitchY()) * (0.5+PixelNumber%(myDIGPlane->GetNpixelsX()));
758  Ypix = (myDIGPlane->GetPitchX()) * (0.5+PixelNumber/(myDIGPlane->GetNpixelsX()));
759 
760 }
761 //______________________________________________________________________________
762 //
763 Double_t DIGParticle::GaussianLaw(Double_t mean, Double_t sigma)
764 {
765  Double_t x;
766  TRandom3 *r3 = new TRandom3(0);
767  GlobalSeed++;
768  r3->SetSeed(GlobalSeed);
769 
770  x = r3->Gaus(mean,sigma);
771 
772  delete r3;
773  return x;
774 }
775 //______________________________________________________________________________
776 //
777 Float_t DIGParticle::GetTotalAnalogCharge(){
778  Float_t TotalCharge = 0;
779  for (Int_t i=0 ; i<fNpixels ; i++){
780  TotalCharge+=fAnalogChargeMap[i];
781  }
782  return TotalCharge;
783 }
784 //______________________________________________________________________________
785 //
786 Int_t DIGParticle::GetTotalDigitalCharge(){
787  Int_t TotalCharge = 0;
788  for (Int_t i=0 ; i<fNpixels ; i++){
789  TotalCharge+=fDigitalChargeMap[i];
790  }
791  return TotalCharge;
792 }
793 //==============================================================================
DIGADC
Definition: digadc.h:36
Generated by   doxygen 1.8.5