Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JEventProcessor_BCAL_TimeCalibration.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // File: JEventProcessor_BCAL_TimeCalibration.cc
4 // Created: Mon Apr 18 15:28:52 CST 2016
5 // Creator: semenov (on Linux selene.phys.uregina.ca 2.6.32-573.18.1.el6.x86_64 x86_64)
6 //
7 
8 #include "JEventProcessor_BCAL_TimeCalibration.h"
9 using namespace jana;
10 
11 #include <TLorentzVector.h>
12 #include "TMath.h"
13 
14 #include "BCAL/DBCALHit.h"
15 #include "BCAL/DBCALTDCHit.h"
16 #include "BCAL/DBCALCluster.h"
17 #include "BCAL/DBCALPoint.h"
18 #include "BCAL/DBCALUnifiedHit.h"
19 #include "BCAL/DBCALGeometry.h"
20 #include "PID/DChargedTrack.h"
21 #include "TRACKING/DTrackTimeBased.h"
22 #include "PID/DEventRFBunch.h"
23 #include "PID/DDetectorMatches.h"
24 #include "DAQ/Df250PulsePedestal.h" // Needed for pulse peak information
25 #include "DAQ/Df250PulseIntegral.h" // Needed for pulse integral !!! added
26 #include "BCAL/DBCALDigiHit.h"
27 
28 #include "DANA/DApplication.h"
29 #include "BCAL/DBCALShower.h"
30 #include "BCAL/DBCALTruthShower.h"
31 #include "FCAL/DFCALShower.h"
32 #include "TRACKING/DMCThrown.h"
33 #include "ANALYSIS/DAnalysisUtilities.h"
34 #include "PID/DVertex.h"
35 #include "TAGGER/DTAGHHit.h"
36 #include "TAGGER/DTAGMHit.h"
37 //#include "TRACKING/DTrackFinder.h"
38 #include "TRIGGER/DL1Trigger.h"
39 
40 #include <TFile.h>
41 #include <TH1.h>
42 #include <TH2.h>
43 #include <TCanvas.h>
44 #include <TF1.h>
45 #include <TStyle.h>
46 #include <TChain.h>
47 #include <TPostScript.h>
48 #include <TLatex.h>
49 #include <TMath.h>
50 #include <TPaveLabel.h>
51 #include <TGraphErrors.h>
52 #include <TLine.h>
53 
54 #include <iostream>
55 #include <fstream>
56 #include <stdlib.h>
57 #include <math.h>
58 #include <stdio.h>
59 #include <string.h>
60 #include "TROOT.h"
61 
62 #include <mutex>
63 #include <vector>
64 #include <deque>
65 #include <string>
66 #include <fstream>
67 #include <iostream>
68 #include <algorithm>
69 #include <stdio.h>
70 #include <stdlib.h>
71 
72 // Routine used to create our JEventProcessor
73 #include <JANA/JApplication.h>
74 #include <JANA/JFactory.h>
75 
76 static mutex mtx;
77 
78 static int ievent;
79 
80 static int iRunNumber;
81 
82 static double R0[5]={0., 65.089, 67.1464, 71.2612, 77.4334}; //layer-dependent radius of the beginning cell position
83 static double DD[5]={0., 2.0574, 4.1148, 6.1722, 8.23}; //layer-dependent cell thickness
84 
85 static int busyShower[1000];
86 
87 static double cvel = 29.9792458 ; //Speed-of-Light (cm/ns)
88 
89 static unsigned short int aU[48][4][3][5000],aD[48][4][3][5000],tU[48][4][3][5000],tD[48][4][3][5000];
90 static unsigned short int Z0[48][4][3][5000], zH[48][4][3][5000], RUD[48][4][3][5000];
91 static float Ene[48][4][3][5000];
92 
93 static double r2U,r2D,rUD,zUD, cosPHI,sinPHI,PHIp,PHIs,PHIdiff;
94 static unsigned short int ipo[48][4][3];
95 static double itU, itD;
96 static double Ze[2]={17.,407.}; //CORRECT Z-position of BCAL ends (cm) May 19, 2016
97 static int imo,ise,ila,ip;
98 static int pdfps = 0; //PDF(0) or PostScript(1) picture
99 
100 static TGraphErrors * gr01;
101 static TGraphErrors * gr11;
102 static TGraphErrors * gr02;
103 static TGraphErrors * gr12;
104 
105 static TGraph * gr1;
106 static TGraph * gr2;
107 static TGraph * gr3;
108 
109 static TH1F * h1;
110 static TH1F * h2;
111 static TH1F * h3;
112 static TH1F * h4;
113 static TH1F * his4;
114 static TH1F * his44;
115 
116 static TF1 * fgg10;
117 static TF1 * fgg20;
118 static TF1 * fgau;
119 static TF1 * fgau2;
120 static TF1 * fz1;
121 static TF1 * fz2;
122 
123 static char * ftxt;
124 static char * ftxt2;
125 static char * ftxt3;
126 static char * ftxt4;
127 static char * hisname;
128 static char * hisname2;
129 
130 static TLatex * txt1;
131 static TLatex * txt2;
132 static TLatex * txt3;
133 static TLatex * txt4;
134 
135 static TLine * lin1;
136 
137 static double tADC_U, tADC_D;
138 static double E1[100], E3[100], E13[100];
139 
140 static double dZs1[3]={-11., -9.5, -23.};
141 static double dZs2[3]={ 9.0, 10.5, 37.};
142 
143 //*********************************************************************************
144 static Double_t tzfunc (Double_t *x, Double_t *par) {
145  Double_t x1=x[0]-212.;
146  Double_t func = par[0]+par[1]*x1;
147  return func;
148 }
149 //*********************************************************************************
150 static Double_t twfunc2 (Double_t *x, Double_t *par) {
151  Double_t t0=par[0];
152  Double_t x1=x[0]-t0;
153 
154  Double_t func = par[1] + par[2]*pow(x1,par[3]) + par[4]*pow(x1,par[5]);
155  return func;
156 }
157 //*********************************************************************************
158 static Double_t gfunc (Double_t *x, Double_t *par) {
159  Double_t x1=x[0]-par[1];
160 
161  Double_t func;
162 
163  func = par[0]*exp(-x1*x1/(2.*par[2]*par[2]));
164  return func;
165 }
166 //*********************************************************************************
167 
168 extern "C"{
169 void InitPlugin(JApplication *app){
170  InitJANAPlugin(app);
171  app->AddProcessor(new JEventProcessor_BCAL_TimeCalibration());
172 }
173 } // "C"
174 
175 //------------------
176 // JEventProcessor_BCAL_TimeCalibration (Constructor)
177 //------------------
178 JEventProcessor_BCAL_TimeCalibration::JEventProcessor_BCAL_TimeCalibration()
179 {
180 
181 }
182 
183 //------------------
184 // ~JEventProcessor_BCAL_TimeCalibration (Destructor)
185 //------------------
186 JEventProcessor_BCAL_TimeCalibration::~JEventProcessor_BCAL_TimeCalibration()
187 {
188 
189 }
190 
191 //------------------
192 // init
193 //------------------
194 jerror_t JEventProcessor_BCAL_TimeCalibration::init(void)
195 {
196  // This is called once at program startup. If you are creating
197  // and filling historgrams in this plugin, you should lock the
198  // ROOT mutex like this:
199  //
200  // japp->RootWriteLock();
201  // ... fill historgrams or trees ...
202  // japp->RootUnLock();
203  //
204 
205  ievent = 0;
206 
207  for (imo=0; imo<48; imo++) {
208  for (ise=0; ise<4; ise++) {
209  for (ila=0; ila<3; ila++) {
210  ipo[imo][ise][ila] = 0; // Null the counter
211  }
212  }
213  }
214 
215  return NOERROR;
216 }
217 
218 //------------------
219 // brun
220 //------------------
221 jerror_t JEventProcessor_BCAL_TimeCalibration::brun(JEventLoop *locEventLoop, int32_t runnumber)
222 {
223  // This is called whenever the run number changes
224  iRunNumber = runnumber;
225  return NOERROR;
226 }
227 
228 //------------------
229 // evnt
230 //------------------
231 jerror_t JEventProcessor_BCAL_TimeCalibration::evnt(JEventLoop *locEventLoop, uint64_t eventnumber)
232 {
233  // This is called for every event. Use of common resources like writing
234  // to a file or filling a histogram should be mutex protected. Using
235  // loop->Get(...) to get reconstructed objects (and thereby activating the
236  // reconstruction algorithm) should be done outside of any mutex lock
237  // since multiple threads may call this method at the same time.
238  // Here's an example:
239  //
240  // vector<const MyDataClass*> mydataclasses;
241  // loop->Get(mydataclasses);
242  //
243  // japp->RootWriteLock();
244  // ... fill historgrams or trees ...
245  // japp->RootUnLock();
246 
247 static int goodFlag = 0;
248 static double ttU, ttD;
249 static int aaU,aaD, hend,hmodule,hsector,hlayer;
250 static double t0,x0,y0,z0, tt0;
251 static double Xs,Ys,Zs, Xh,Yh,Zh, Rh,rU,rD,zU,zD ;
252 static double x,y,z,R, dPhi,dZ, E;
253 
254  goodFlag = 0;
255 
256 //****************
257 //CHECK TRIGGER TYPE
258 // const DL1Trigger* locTrigger = NULL;
259 // locEventLoop->GetSingle(locTrigger);
260 
261  uint32_t trig_mask=0;
262  const DL1Trigger *trig_words = NULL;
263  try {
264  locEventLoop->GetSingle(trig_words);
265  } catch(...) {};
266  if (trig_words) {
267  trig_mask = trig_words->trig_mask;
268  } else {
269  trig_mask=0;
270  }
271 
272 
273 
274 //****************
275 
276  vector<const DTrackTimeBased*> locTrackTimeBased;
277  locEventLoop->Get(locTrackTimeBased);
278  vector<const DChargedTrackHypothesis*> locTrackHypothesis;
279  locEventLoop->Get(locTrackHypothesis);
280 
281  vector<const DBCALShower*> locBCALShowers;
282  locEventLoop->Get(locBCALShowers);
283 
284  vector<const DBCALShower*> matchedShowers;
285 
286  vector<const DBCALCluster*> locBCALClusters;
287 // locEventLoop->Get(locBCALClusters);
288 
289  vector<const DBCALPoint*> matchedBCALPoints;
290 
291  vector<const DBCALPoint*> locBCALPoints;
292  locEventLoop->Get(locBCALPoints);
293 
294  vector<const DVertex*> kinfitVertex;
295  locEventLoop->Get(kinfitVertex);
296 
297  vector<const DEventRFBunch*> RFbunch;
298  locEventLoop->Get(RFbunch);
299  int Nvotes;
300 
301  for (unsigned int ij=0; ij<RFbunch.size(); ij++) {
302  Nvotes = RFbunch[ij]->dNumParticleVotes;
303  }
304 
305  vector<const DBCALUnifiedHit *> locBCALUnifiedHit;
306  locEventLoop->Get(locBCALUnifiedHit);
307 
308  vector<const DBCALUnifiedHit *> matchedBCALUnifiedHit;
309 
310  const DBCALHit * thisADCHit;
311 
312  const DBCALTDCHit * thisTDCHit;
313 
314 // Here we want to match tracks to the BCAL so that we can throw away charged particles in the BCAL to decrease combinatorics. To do this we look over time based tracks and swim each track in an event out to each BCAL shower in each event. Then check to see if the BCAL shower and charged track are within some common z and phi region of the calorimeter. If they are then we fill a vector with those showers that are geometrically matched to a charged track. This vector will be called later to skip such showers.
315 
316  DVector3 mypos(0.0,0.0,0.0);
317 
318  // Only one thread at a time from here to end of method
319  lock_guard<mutex> lck(mtx);
320 
321  for (int ij=0; ij<1000; ij++) busyShower[ij]=0;
322 
323  for (unsigned int i=0; i < locTrackTimeBased.size() ; ++i){
324  tt0 = locTrackHypothesis[i]->t0(); if (tt0<-10000.||tt0>10000.) goodFlag = 1;
325  for (unsigned int j=0; j< locBCALShowers.size(); ++j){
326 
327  x = locBCALShowers[j]->x;
328  y = locBCALShowers[j]->y;
329  z = locBCALShowers[j]->z;
330  DVector3 pos_bcal(x,y,z); // Vector that has the BCAL shower centroid position so we know where to try and match charged track
331  R = pos_bcal.Perp();
332  //locTrackTimeBased[i]->rt->GetIntersectionWithRadius(R, mypos); // This will swim a track out of the drift chambers to some radius that is defined by pos_bcal and returns a vector mypos that will give the location of the charged track at the redius we input (R)
333  locTrackTimeBased[i]->GetProjection(SYS_BCAL, mypos); // This will swim a track out of the drift chambers to some radius that is defined by pos_bcal and returns a vector mypos that will give the location of the charged track at the redius we input (R)
334 
335  dPhi = TMath::Abs(mypos.Phi()-pos_bcal.Phi()); // find how far away in phi and z that the shower and track are
336  dZ = TMath::Abs(mypos.Z() - z);
337  if(dZ<60.0&&dPhi<0.50&&busyShower[j]==0) { // if the track and shower are close enough in phi and z then fill the matchedshowers vector
338  busyShower[j]=1; // Mark the shower to be NOT used
339  }
340 
341  }
342  }
343 
344  for (unsigned int j=0; j< locBCALShowers.size(); ++j){
345  if (busyShower[j]==0) matchedShowers.push_back(locBCALShowers[j]); //Showers that are NOT matched with tracks
346  }
347 
348  t0 = 99999999.;
349  for (unsigned int ij=0; ij<kinfitVertex.size(); ij++) {
350  t0= kinfitVertex[ij]->dSpacetimeVertex.T();
351  x0= kinfitVertex[ij]->dSpacetimeVertex.X();
352  y0= kinfitVertex[ij]->dSpacetimeVertex.Y();
353  z0= kinfitVertex[ij]->dSpacetimeVertex.Z();
354  }
355 
356 //Check on the signal in the first layer and E1/E3 ratio
357  for (unsigned int ij=0; ij< matchedShowers.size(); ++ij){ //<------------------------- ij loop on Neutral Showers
358  E1[ij]=0;
359  E3[ij]=0;
360  E13[ij] = -999999.;
361 
362  if (goodFlag==0) {
363  matchedShowers[ij]->Get(matchedBCALPoints);
364  for (unsigned int il=0; il<matchedBCALPoints.size(); il++) { //<------------------------- il loop on Points
365  matchedBCALPoints[il]->Get(matchedBCALUnifiedHit);
366  for (unsigned int ik=0; ik<matchedBCALUnifiedHit.size(); ik++) { //<---------------------- ik loop on Hits
367  matchedBCALUnifiedHit[ik]->GetSingle(thisADCHit);
368  hlayer = thisADCHit->layer;
369  aaU = thisADCHit->pulse_peak;
370  }//<------------------------------------------------------------------------ End of ik loop on Hits
371 
372  if (hlayer==1) E1[ij] += matchedBCALPoints[il]->E();
373  if (hlayer==3) E3[ij] += matchedBCALPoints[il]->E();
374 
375  }//<------------------------------------------------------------------------ End of il loop on Points
376  }
377  if (E3[ij]>0) {
378  E13[ij] = E1[ij]/E3[ij];
379  } else {
380  E13[ij] = 999999.;
381  }
382  } //<----------------------------------------------------------------------------- End of ij loop on Neutral Showers
383 //End of check on the signal in the first layer
384 
385  for (unsigned int ij=0; ij< matchedShowers.size(); ++ij){ //<------------------------- ij loop on Neutral Showers
386  if (goodFlag==0) {
387  Xs = matchedShowers[ij]->x;
388  Ys = matchedShowers[ij]->y;
389  Zs = matchedShowers[ij]->z;
390  PHIs = atan2(Ys,Xs);
391 
392  matchedShowers[ij]->Get(matchedBCALPoints);
393 
394  for (unsigned int il=0; il<matchedBCALPoints.size(); il++) { //<------------------------- il loop on Points
395 
396  Rh = matchedBCALPoints[il]->r();
397  PHIp = matchedBCALPoints[il]->phi();
398  cosPHI = cos(PHIp);
399  sinPHI = sin(PHIp);
400 
401  Zh = matchedBCALPoints[il]->z();
402 
403  E = matchedBCALPoints[il]->E();
404 
405  matchedBCALPoints[il]->Get(matchedBCALUnifiedHit);
406 
407  ttU=-9999.;
408  ttD=-9999.;
409  aaU=-9999.;
410  aaD=-9999.;
411  tADC_U = -9999.;
412  tADC_D = -9999.;
413 
414  for (unsigned int ik=0; ik<matchedBCALUnifiedHit.size(); ik++) { //<---------------------- ik loop on Hits
415  matchedBCALUnifiedHit[ik]->GetSingle(thisADCHit);
416  matchedBCALUnifiedHit[ik]->GetSingle(thisTDCHit);
417 
418  hend = thisADCHit->end;
419  hmodule = thisADCHit->module;
420  hsector = thisADCHit->sector;
421  hlayer = thisADCHit->layer;
422 
423  if (hend==0) tADC_U = matchedBCALUnifiedHit[ik]->t_ADC;
424  if (hend==1) tADC_D = matchedBCALUnifiedHit[ik]->t_ADC;
425 
426  if (hend==0) {
427  aaU = thisADCHit->pulse_peak;
428  rU=0.;
429  if (aaU>0) rU = R0[hlayer]+DD[hlayer]*2./aaU; //<---- Threshold of 2 counts
430  zU = z0 + (Zh+65.-z0) * rU/Rh;
431  Xh = rU*cosPHI;
432  Yh = rU*sinPHI;
433  r2U = (Xh-x0)*(Xh-x0)+(Yh-y0)*(Yh-y0);
434  }
435 
436  if (hend==1) {
437  aaD = thisADCHit->pulse_peak;
438  rD=0.;
439  if (aaD>0) rD = R0[hlayer]+DD[hlayer]*2./aaD; //<---- Threshold of 2 counts
440  zD = z0 + (Zh+65.-z0) * rD/Rh;
441  Xh = rD*cosPHI;
442  Yh = rD*sinPHI;
443  r2D = (Xh-x0)*(Xh-x0)+(Yh-y0)*(Yh-y0);
444  }
445 
446  if (matchedBCALUnifiedHit[ik]->has_TDC_hit && Nvotes>=2) {
447  if (hend==0) ttU = thisTDCHit->t;
448  if (hend==1) ttD = thisTDCHit->t;
449  }
450  } //end of ik loop
451 
452  imo = hmodule-1;
453  ise = hsector-1;
454  ila = hlayer -1;
455  ip = ipo[imo][ise][ila];
456 
457  PHIdiff = abs(PHIp-PHIs);
458  if (PHIdiff>3.1415927) PHIdiff = 6.2831854 - PHIdiff;
459  itU = 100.*(ttU - tt0);
460  itD = 100.*(ttD - tt0);
461  zUD = 0.5*(zU+zD);
462  rUD = 0.5*(sqrt(r2U)+sqrt(r2D));
463 
464 //japp->RootFillLock(this);
465  if (ip<5000 // To not overflow array
466  && itU>0. && itD>0. && itU<60000. && ttD<60000.
467  && aaU>5. && aaD>5. && aaU<60000. && aaD<60000. // Both ends made signals
468  && matchedBCALPoints.size()>2 && matchedBCALPoints.size()<30 // Reasonable multiplicity
469  && hlayer>0 && hlayer<4 // Layers 1-3
470  && trig_mask==4 // BCAL trigger
471  && zUD>(Ze[0]+10.) && zUD<(Ze[1]-10.) // Z inside BCAL
472  && rUD>60. && rUD<100.
473  && E13[ij]>0.1 && E13[ij]<10000. // Suppressed neutrons as on 170706
474  && abs(z0-65.)<15. // Vertex is inside the target
475  && PHIdiff<0.035
476  && (zUD-Zs)>dZs1[ila] && (zUD-Zs)<dZs2[ila]) { // Hit is not far away from the Shower
477 
478  aU[imo][ise][ila][ip] = 1*aaU;
479  aD[imo][ise][ila][ip] = 1*aaD;
480  tU[imo][ise][ila][ip] = 1*itU;
481  tD[imo][ise][ila][ip] = 1*itD;
482  Z0[imo][ise][ila][ip] = 100*z0;
483  zH[imo][ise][ila][ip] = 100*zUD;
484  RUD[imo][ise][ila][ip] = 100*rUD;
485  Ene[imo][ise][ila][ip] = E;
486 
487  ipo[imo][ise][ila]++;
488  ievent++;
489  }
490 //japp->RootFillUnLock(this);
491 
492  } //end of il loop
493 
494  } //if(GoodFlag end
495  } //end of ij loop
496 
497  return NOERROR;
498 }
499 
500 //------------------
501 // erun
502 //------------------
503 jerror_t JEventProcessor_BCAL_TimeCalibration::erun(void)
504 {
505  // This is called whenever the run number changes, before it is
506  // changed to give you a chance to clean up before processing
507  // events from the next run number.
508  return NOERROR;
509 }
510 
511 //------------------
512 // fini
513 //------------------
514 //*********************************************************************************
515 jerror_t JEventProcessor_BCAL_TimeCalibration::fini(void)
516 {
517  // Called before program exit after event processing is finished.
518 
519 //*********************************************************************************
520 static double f1p0,f1p1,f1p2,f1p3,f1p4,f1p5;
521 static double ff1p0,ff1p1,ff1p2,ff1p3,ff1p4,ff1p5;
522 
523 static double f2p0,f2p1,f2p2,f2p3,f2p4,f2p5;
524 static double ff2p0,ff2p1,ff2p2,ff2p3,ff2p4,ff2p5;
525 
526 static double f3p0,f3p1,f3p2,f3p3,f3p4,f3p5;
527 static double ff3p0,ff3p1,ff3p2,ff3p3,ff3p4,ff3p5;
528 
529 static double f4p0,f4p1,f4p2,f4p3,f4p4,f4p5;
530 static double ff4p0,ff4p1,ff4p2,ff4p3,ff4p4,ff4p5;
531 
532 static double xp10,xp11,xp12, xp20,xp21,xp22;
533 static double xp1,xp2,xp3,xp4,xp5, ytmp, ytmp2, Zh, pL;
534 static double veff2, zp0,zp1, veff22, zzp0,zzp1, vef3,vef33;
535 
536 int ipo3;
537 //************************************************************************************
538 //double cvel = 29.9792458 ; //Speed-of-Light (cm/ns)
539 //double veff = 16.5; //Init veff (cm/ns)
540 static double vef1 = 15.9; //Init veff (cm/ns)
541 static double vef2 = 16.9; //Init veff (cm/ns)
542 //double Ze[2]={17.,407.}; //CORRECT Z-position of BCAL ends (cm) May 19, 2016
543 
544 static double y1[5000],x1[5000],ey1[5000],ex1[5000],z1[5000],y10[5000],y20[5000];
545 static double yy1[5000],yy2[5000],xx1[5000],eyy1[5000],exx1[5000];
546 static double y2[5000],x2[5000],ey2[5000],ex2[5000];
547 static int ipo2 = 0;
548 
549 static double xaU,xaD,xtU,xtD,xZ0,xZH,xRUD;
550 
551 static double gme,gsi,egme,egsi, gme2,gsi2,egme2,egsi2, cont1,cont4;
552 
553 static float itt1, itt2, itt3, itt4, xkoe;
554 //************************************************************************************
555 
556 static char *outcalname = new char[90];
557  sprintf(outcalname, "TDCcalib_run%d.vec",iRunNumber);
558  ofstream dat100(outcalname);
559 
560 static char *outcalname2 = new char[90];
561  sprintf(outcalname2, "ccdb-BCAL_effective_velocities_regina-run%d.dat",iRunNumber);
562  ofstream datveff(outcalname2);
563 
564 static char *outcalname3 = new char[90];
565  sprintf(outcalname3, "ccdb-BCAL_timewalk_constants_regina-run%d.dat",iRunNumber);
566  ofstream datwalk(outcalname3);
567 
568 static char *outcalname4 = new char[90];
569  sprintf(outcalname4, "WARNINGS-run%d.vec",iRunNumber);
570  ofstream datwarn(outcalname4);
571 
572 //************************************************************************************
573 //************************************************************************************
574 
575 //Canvas Definition
576 static TCanvas *c1 = new TCanvas("c1", " c1",0,0,699,499);
577 
578 //************************************************************************************
579 
580 static char *pdfname = new char[90];
581 static char *pdfname1 = new char[90];
582 static char *pdfname2 = new char[90];
583  if (pdfps==0) {
584  sprintf(pdfname, "TDCcalib_run%d.pdf",iRunNumber);
585  sprintf(pdfname1, "TDCcalib_run%d.pdf[",iRunNumber);
586  sprintf(pdfname2, "TDCcalib_run%d.pdf]",iRunNumber);
587  } else {
588  sprintf(pdfname, "TDCcalib_run%d.ps",iRunNumber);
589  sprintf(pdfname1, "TDCcalib_run%d.ps[",iRunNumber);
590  sprintf(pdfname2, "TDCcalib_run%d.ps]",iRunNumber);
591  }
592  c1->Update(); c1->Print(pdfname1);
593 //***************************************************
594 
595 for (int imo=0; imo<48; imo++) { //Start of imo loop
596 for (int ila=0; ila<3; ila++) { //Start of ila loop
597 for (int ise=0; ise<4; ise++) { //Start of ise loop
598 //************************************************************************************
599 
600 h1 = new TH1F("h1","Middle U",150,-50.,150.);
601 h2 = new TH1F("h2","High U", 150,-50.,150.);
602 
603 h3 = new TH1F("h3","Middle D",150,-50.,150.);
604 h4 = new TH1F("h4","High D", 150,-50.,150.);
605 
606  hisname = new char[90];
607  sprintf(hisname, "#gamma Hits: Run %d, Module=%02i, Sector=%i, Layer=%i, End=0",iRunNumber,imo+1,ise+1,ila+1);
608  hisname2 = new char[90];
609  sprintf(hisname2, "#gamma Hits: Run %d, Module=%02i, Sector=%i, Layer=%i, End=1",iRunNumber,imo+1,ise+1,ila+1);
610 
611 //*******
612  his4 = new TH1F("his4"," ",200,-8.,12.);
613  his4->GetXaxis()->SetTitle("Corrected #DeltaTime (ns)");
614  his4->GetYaxis()->SetTitle("Counts");
615  his4->SetTitle(hisname);
616 // his4->SetMaximum(20.);
617  his4->SetLineColor(1);
618  his4->SetLineWidth(1.9);
619  his4->SetFillStyle(1001);
620  his4->SetFillColor(kGreen-7);
621 
622 //*******
623  his44 = new TH1F("his44"," ",200,-8.,12.);
624  his44->GetXaxis()->SetTitle("Corrected #DeltaTime (ns)");
625  his44->GetYaxis()->SetTitle("Counts");
626  his44->SetTitle(hisname2);
627 // his44->SetMaximum(20.);
628  his44->SetLineColor(1);
629  his44->SetLineWidth(1.9);
630  his44->SetFillStyle(1001);
631  his44->SetFillColor(kGreen-7);
632 
633 //************************************************************************************
634 //************************************************************************************
635 //1st Iteration
636 //************************************************************************************
637 //************************************************************************************
638 { cout<<endl<<"************ Imo="<<imo<<" Ise="<<ise<<" Ila="<<ila<<" *******************************"<<endl<<endl;
639 
640  ipo2 = 0;
641  for (unsigned short int ip=0; ip<ipo[imo][ise][ila]; ip++) {
642 
643  xaU = 1. *aU[imo][ise][ila][ip];
644  xaD = 1. *aD[imo][ise][ila][ip];
645  xtU = 0.01 *(tU[imo][ise][ila][ip]);
646  xtD = 0.01 *(tD[imo][ise][ila][ip]);
647  xZ0 = 0.01 *Z0[imo][ise][ila][ip];
648  xZH = 0.01 *zH[imo][ise][ila][ip];
649  xRUD= 0.01 *RUD[imo][ise][ila][ip];
650 
651  pL =(xZH-xZ0)*(xZH-xZ0);
652  pL =sqrt(pL+xRUD*xRUD);
653 
654  y1[ipo2]=xtU-pL/cvel-abs(xZH-Ze[0])/vef1;
655  y2[ipo2]=xtD-pL/cvel-abs(xZH-Ze[1])/vef2;
656  if (y1[ipo2]>-20.&&y1[ipo2]<30.&&y2[ipo2]>-20.&&y2[ipo2]<30.&&Ene[imo][ise][ila][ip]>0.
657  &&xZH>Ze[0]&&xZH<Ze[1]) {
658  x1[ipo2]=xaU;
659 // ey1[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
660  ey1[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
661  ex1[ipo2]= 0.1;
662  z1[ipo2] = xZH;
663 
664  x2[ipo2]=xaD;
665 // ey2[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
666  ey2[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
667  ex2[ipo2]= 0.1;
668 
669  ipo2++;
670  }
671  }
672 
673  itt1 = 1.*ipo2;
674 
675  for (int ip=0; ip<ipo2; ip++) {
676  if (x1[ip]>10.&&x1[ip]<20.) h1->Fill(y1[ip]);
677  if (x1[ip]>50.&&x1[ip]<100.) h2->Fill(y1[ip]);
678 
679  if (x2[ip]>10.&&x2[ip]<20.) h3->Fill(y2[ip]);
680  if (x2[ip]>50.&&x2[ip]<100.) h4->Fill(y2[ip]);
681  }
682 
683  cout<<endl<<"1st Iteration: "<<ipo2<<" events in the plot."<<endl<<endl;
684  cont1 = ipo2;
685 
686  xp10 = (h2->GetMean()) - 7.;
687  xp12 = log((h2->GetMean()-xp10)/(h1->GetMean()-xp10))/log(60./15.) -0.1;
688  xp11 = (h2->GetMean()-xp10)/pow(60.,xp12);
689 
690  cout<<endl<<"xp10="<<xp10<<" xp11="<<xp11<<" xp12="<<xp12<<endl;
691 
692  xp20 = (h4->GetMean()) - 7.;
693  xp22 = log((h4->GetMean()-xp20)/(h3->GetMean()-xp20))/log(60./15.);
694  xp21 = (h4->GetMean()-xp20)/pow(60.,xp22);
695 
696  cout<<"xp20="<<xp20<<" xp21="<<xp21<<" xp22="<<xp22<<endl;
697 
698 
699 //*******
700 //********
701  fgg10 = new TF1("fgg10",twfunc2,5.,2500.,6);
702 
703  fgg10->SetParameter(0,2.5); // Pole parameter
704  fgg10->SetParLimits(0,0.,4.9); // Pole parameter
705 
706  xp1=xp10 ;
707  fgg10->SetParameter(1,xp1); // Level
708  fgg10->SetParLimits(1,xp1-4.,xp1+4.);
709 
710  xp2=30.;
711  fgg10->SetParameter(2,xp2); // Amp1
712  fgg10->SetParLimits(2,10.,50.);
713 
714  xp3=-0.4;
715  fgg10->SetParameter(3,xp3); // Pow1
716  fgg10->SetParLimits(3,-0.55,-0.2);
717 
718  xp4 = 10000.;
719  fgg10->SetParameter(4,xp4); // Amp2
720  fgg10->SetParLimits(4,100.,1000000.);
721 
722  xp5 = -3.5;
723  fgg10->SetParameter(5,xp5); // Pow2
724  fgg10->SetParLimits(5,-5.5,-2.5);
725 
726  fgg10->SetLineColor(4);
727  fgg10->SetLineWidth(1.);
728 
729  gr01 = new TGraphErrors(ipo2,x1,y1,ex1,ey1);
730  gr01->Fit("fgg10","MER");
731  f1p0 = fgg10->GetParameter(0);
732  f1p1 = fgg10->GetParameter(1);
733  f1p2 = fgg10->GetParameter(2);
734  f1p3 = fgg10->GetParameter(3);
735  f1p4 = fgg10->GetParameter(4);
736  f1p5 = fgg10->GetParameter(5);
737 
738  cout<<"-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-1-"<<endl;
739 //********
740  fgg20 = new TF1("fgg20",twfunc2,5.,2500.,6);
741 
742  fgg20->SetParameter(0,2.5); // Pole parameter
743  fgg20->SetParLimits(0,0.,4.9); // Pole parameter
744 
745  xp1=xp20 ;
746  fgg20->SetParameter(1,xp1); // Level
747  fgg20->SetParLimits(1,xp1-4.,xp1+4.);
748 
749  xp2=30.;
750  fgg20->SetParameter(2,xp2); // Amp1
751  fgg20->SetParLimits(2,10.,50.);
752 
753  xp3=-0.4;
754  fgg20->SetParameter(3,xp3); // Pow1
755  fgg20->SetParLimits(3,-0.55,-0.2);
756 
757  xp4 = 10000.;
758  fgg20->SetParameter(4,xp4); // Amp2
759  fgg20->SetParLimits(4,100.,1000000.);
760 
761  xp5 = -3.5;
762  fgg20->SetParameter(5,xp5); // Pow2
763  fgg20->SetParLimits(5,-5.5,-2.5);
764 
765  fgg20->SetLineColor(4);
766  fgg20->SetLineWidth(1.);
767 
768  gr11 = new TGraphErrors(ipo2,x2,y2,ex2,ey2);
769  gr11->Fit("fgg20","MER");
770  ff1p0 = fgg20->GetParameter(0);
771  ff1p1 = fgg20->GetParameter(1);
772  ff1p2 = fgg20->GetParameter(2);
773  ff1p3 = fgg20->GetParameter(3);
774  ff1p4 = fgg20->GetParameter(4);
775  ff1p5 = fgg20->GetParameter(5);
776 
777 //********
778 //********
779 
780  int ipo3 = 0;
781  for (int ii=0; ii<ipo2; ii++) {
782  ytmp=y1[ii]-f1p1-f1p2*pow((x1[ii]-f1p0),f1p3)-f1p4*pow((x1[ii]-f1p0),f1p5);
783  ytmp2=y2[ii]-ff1p1-ff1p2*pow((x2[ii]-ff1p0),ff1p3)-ff1p4*pow((x2[ii]-ff1p0),ff1p5);
784 
785  if (abs(ytmp)<3.&&abs(ytmp2)<3.) {
786  y1[ipo3]=ytmp;
787  y2[ipo3]=ytmp2;
788  ey1[ipo3]=ey1[ii];
789  x1[ipo3]=z1[ii];
790  ex1[ipo3]=0.01;
791  ipo3++;
792  }
793  }
794 
795  fz1 = new TF1("fz1",tzfunc, 10.,420.,2);
796  fz1->SetLineColor(2);
797  fz1->SetLineWidth(1.);
798 
799  fz1->SetParameter(0,0.); // Extra-Shift
800  fz1->SetParLimits(0,-0.9,0.9);
801  fz1->SetParameter(1,0.0); // Slope
802  fz1->SetParLimits(1,-0.02,0.02);
803 
804  gr02 = new TGraphErrors(ipo3,x1,y1,ex1,ey1);
805  gr02->Fit("fz1","MER");
806  zp0 = fz1->GetParameter(0);
807  zp1 = fz1->GetParameter(1);
808  veff2 = vef1/(1+vef1*zp1);
809 
810 //********
811 
812  fz2 = new TF1("fz2",tzfunc, 10.,420.,2);
813  fz2->SetLineColor(2);
814  fz2->SetLineWidth(1.);
815 
816  fz2->SetParameter(0,0.); // Extra-Shift
817  fz2->SetParLimits(0,-0.9,0.9);
818  fz2->SetParameter(1,0.0); // Slope
819  fz2->SetParLimits(1,-0.02,0.02);
820 
821  gr12 = new TGraphErrors(ipo3,x1,y2,ex1,ey1);
822  gr12->Fit("fz2","MER");
823  zzp0 = fz2->GetParameter(0);
824  zzp1 = fz2->GetParameter(1);
825  veff22 = vef2/(1-vef2*zzp1);
826 
827 
828 delete gr01;
829 delete gr02;
830 delete gr11;
831 delete gr12;
832 delete fgg10;
833 delete fgg20;
834 delete fz1;
835 delete fz2;
836 
837 }
838  if (veff2<15.7) veff2=16.2;
839  if (veff2>16.7) veff2=16.2;
840  if (veff22<16.2) veff22=16.8;
841  if (veff22>17.2) veff22=16.8;
842 
843  cout<<endl<<"1st veff2="<<veff2<<" veff22="<<veff22<<endl<<endl;
844 
845 //*********************************************************************************
846 //*********************************************************************************
847 // 2nd Iteration
848 //*********************************************************************************
849 //*********************************************************************************
850 {
851  xkoe = 1.;
852 metka1: ipo2 = 0;
853 
854  for (unsigned short int ip=0; ip<ipo[imo][ise][ila]; ip++) {
855 
856  xaU = 1. *aU[imo][ise][ila][ip];
857  xaD = 1. *aD[imo][ise][ila][ip];
858  xtU = 0.01 *(tU[imo][ise][ila][ip]);
859  xtD = 0.01 *(tD[imo][ise][ila][ip]);
860  xZ0 = 0.01 *Z0[imo][ise][ila][ip];
861  xZH = 0.01 *zH[imo][ise][ila][ip];
862  xRUD= 0.01 *RUD[imo][ise][ila][ip];
863 
864  ytmp = xtU-f1p1-f1p2*pow((xaU-f1p0),f1p3)-f1p4*pow((xaU-f1p0),f1p5);
865  ytmp2= xtD-ff1p1-ff1p2*pow((xaD-ff1p0),ff1p3)-ff1p4*pow((xaD-ff1p0),ff1p5);
866 
867  Zh = ((ytmp-ytmp2)*veff2*veff22 + Ze[0]*veff22 + Ze[1]*veff2)/(veff2 + veff22);
868  pL = (Zh-xZ0)*(Zh-xZ0);
869  pL =sqrt(pL + xRUD*xRUD);
870 
871  ytmp = ytmp-pL/cvel-abs(Zh-Ze[0])/veff2;
872  ytmp2 = ytmp2-pL/cvel-abs(Zh-Ze[1])/veff22;
873 
874  y1[ipo2]=xtU-pL/cvel-abs(Zh-Ze[0])/veff2;
875  y2[ipo2]=xtD-pL/cvel-abs(Zh-Ze[1])/veff22;
876 
877 
878  if (y1[ipo2]>-20.&&y1[ipo2]<30.&&y2[ipo2]>-20.&&y2[ipo2]<30.&&Ene[imo][ise][ila][ip]>0.&&abs(ytmp)<8.*xkoe&&abs(ytmp2)<8.*xkoe
879  &&Zh>Ze[0]&&Zh<Ze[1]) {
880  x1[ipo2]=xaU;
881 // ey1[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
882  ey1[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
883  ex1[ipo2]= 0.1;
884  z1[ipo2] = Zh;
885 
886  x2[ipo2]=xaD;
887 // ey2[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
888  ey2[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
889  ex2[ipo2]= 0.1;
890 
891  ipo2++;
892  }
893  }
894 
895  itt2 = 1.*ipo2;
896  if ((itt2/itt1)<0.9) {
897  f2p0 = 2.5;
898  f2p1 = xp10;
899  f2p2 = 30.;
900  f2p3 = -0.4;
901  f2p4 = 100000.;
902  f2p5 = -4.0;
903 
904  ff2p0 = 2.5;
905  ff2p1 = xp20;
906  ff2p2 = 30.;
907  ff2p3 = -0.4;
908  ff2p4 = 100000.;
909  ff2p5 = -4.0;
910  xkoe = xkoe * 1.5;
911  cout<<"Itt2 reset: xkoe = "<<xkoe<<endl;
912  goto metka1;
913  }
914  cout<<endl<<"2nd Iteration: "<<ipo2<<" events in the plot."<<endl<<endl;
915 
916 //*******
917  fgg10 = new TF1("fgg10",twfunc2,5.,2500.,6);
918 
919  fgg10->SetParameter(0,2.5); // Pole parameter
920  fgg10->SetParLimits(0,0.,4.9); // Pole parameter
921 
922  xp1=f1p1 ;
923  fgg10->SetParameter(1,xp1); // Level
924  fgg10->SetParLimits(1,xp1-1.,xp1+1.);
925 
926  xp2=f1p2;
927  fgg10->SetParameter(2,xp2); // Amp1
928  fgg10->SetParLimits(2,0.8*xp2,1.2*xp2);
929 
930  xp3=f1p3;
931  fgg10->SetParameter(3,xp3); // Pow1
932  fgg10->SetParLimits(3,xp3-0.1,xp3+0.1);
933 
934  xp4 = f1p4;
935  fgg10->SetParameter(4,xp4); // Amp2
936  fgg10->SetParLimits(4,0.8*xp4,1.2*xp4);
937 
938  xp5 = f1p5;
939  fgg10->SetParameter(5,xp5); // Pow2
940  fgg10->SetParLimits(5,xp5-0.1,xp5+0.1);
941 
942  fgg10->SetLineColor(4);
943  fgg10->SetLineWidth(1.);
944 
945  gr01 = new TGraphErrors(ipo2,x1,y1,ex1,ey1);
946  gr01->Fit("fgg10","MER");
947  f2p0 = fgg10->GetParameter(0);
948  f2p1 = fgg10->GetParameter(1);
949  f2p2 = fgg10->GetParameter(2);
950  f2p3 = fgg10->GetParameter(3);
951  f2p4 = fgg10->GetParameter(4);
952  f2p5 = fgg10->GetParameter(5);
953 
954  cout<<"-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-2-"<<endl;
955 //********
956  fgg20 = new TF1("fgg20",twfunc2,5.,2500.,6);
957 
958  fgg20->SetParameter(0,2.5); // Pole parameter
959  fgg20->SetParLimits(0,0.,4.9); // Pole parameter
960 
961  xp1=ff1p1 ;
962  fgg20->SetParameter(1,xp1); // Level
963  fgg20->SetParLimits(1,xp1-1.,xp1+1.);
964 
965  xp2=ff1p2;
966  fgg20->SetParameter(2,xp2); // Amp1
967  fgg20->SetParLimits(2,0.8*xp2,1.2*xp2);
968 
969  xp3=ff1p3;
970  fgg20->SetParameter(3,xp3); // Pow1
971  fgg20->SetParLimits(3,xp3-0.1,xp3+0.1);
972 
973  xp4 = ff1p4;
974  fgg20->SetParameter(4,xp4); // Amp2
975  fgg20->SetParLimits(4,0.8*xp4,1.2*xp4);
976 
977  xp5 = ff1p5;
978  fgg20->SetParameter(5,xp5); // Pow2
979  fgg20->SetParLimits(5,xp5-0.1,xp5+0.1);
980 
981  fgg20->SetLineColor(4);
982  fgg20->SetLineWidth(1.);
983 
984  gr11 = new TGraphErrors(ipo2,x2,y2,ex2,ey2);
985  gr11->Fit("fgg20","MER");
986  ff2p0 = fgg20->GetParameter(0);
987  ff2p1 = fgg20->GetParameter(1);
988  ff2p2 = fgg20->GetParameter(2);
989  ff2p3 = fgg20->GetParameter(3);
990  ff2p4 = fgg20->GetParameter(4);
991  ff2p5 = fgg20->GetParameter(5);
992 
993 //********
994 //********
995 
996  ipo3 = 0;
997  for (int ii=0; ii<ipo2; ii++) {
998  ytmp=y1[ii]-f2p1-f2p2*pow((x1[ii]-f2p0),f2p3)-f2p4*pow((x1[ii]-f2p0),f2p5);
999  ytmp2=y2[ii]-ff2p1-ff2p2*pow((x2[ii]-ff2p0),ff2p3)-ff2p4*pow((x2[ii]-ff2p0),ff2p5);
1000 
1001  if (abs(ytmp)<3. && abs(ytmp2)<3.) {
1002  yy1[ipo3]=ytmp;
1003  yy2[ipo3]=ytmp2;
1004  eyy1[ipo3]=ey1[ii];
1005  xx1[ipo3]=z1[ii];
1006  exx1[ipo3]=0.1;
1007  ipo3++;
1008  }
1009  }
1010 
1011  fz1 = new TF1("fz1",tzfunc, 10.,420.,2);
1012  fz1->SetLineColor(2);
1013  fz1->SetLineWidth(1.);
1014 
1015  fz1->SetParameter(0,0.); // Extra-Shift
1016  fz1->SetParLimits(0,-0.9,0.9);
1017  fz1->SetParameter(1,0.0); // Slope
1018  fz1->SetParLimits(1,-0.1,0.1);
1019 
1020  gr02 = new TGraphErrors(ipo3,xx1,yy1,exx1,eyy1);
1021  gr02->Fit("fz1","MER");
1022  zp0 = fz1->GetParameter(0);
1023  zp1 = fz1->GetParameter(1);
1024  veff2 = veff2/(1+veff2*zp1);
1025 
1026 //********
1027 
1028  fz2 = new TF1("fz2",tzfunc, 10.,420.,2);
1029  fz2->SetLineColor(2);
1030  fz2->SetLineWidth(1.);
1031 
1032  fz2->SetParameter(0,0.); // Extra-Shift
1033  fz2->SetParLimits(0,-0.9,0.9);
1034  fz2->SetParameter(1,0.0); // Slope
1035  fz2->SetParLimits(1,-0.02,0.02);
1036 
1037  gr12 = new TGraphErrors(ipo3,xx1,yy2,exx1,eyy1);
1038  gr12->Fit("fz2","MER");
1039  zzp0 = fz2->GetParameter(0);
1040  zzp1 = fz2->GetParameter(1);
1041  veff22 = veff22/(1-veff22*zzp1);
1042 
1043 delete gr01;
1044 delete gr02;
1045 delete gr11;
1046 delete gr12;
1047 delete fgg10;
1048 delete fgg20;
1049 delete fz1;
1050 delete fz2;
1051 
1052 }
1053 
1054  if (veff2<15.7) veff2=16.2;
1055  if (veff2>16.7) veff2=16.2;
1056  if (veff22<16.2) veff22=16.8;
1057  if (veff22>17.2) veff22=16.8;
1058 
1059  cout<<endl<<"2nd veff2="<<veff2<<" veff22="<<veff22<<endl<<endl;
1060 
1061 
1062 //*********************************************************************************
1063 //*********************************************************************************
1064 // 3rd Iteration
1065 //*********************************************************************************
1066 //*********************************************************************************
1067 {
1068  xkoe = 1.;
1069 metka2: ipo2 = 0;
1070 
1071  for (unsigned short int ip=0; ip<ipo[imo][ise][ila]; ip++) {
1072 
1073  xaU = 1. *aU[imo][ise][ila][ip];
1074  xaD = 1. *aD[imo][ise][ila][ip];
1075  xtU = 0.01 *(tU[imo][ise][ila][ip]);
1076  xtD = 0.01 *(tD[imo][ise][ila][ip]);
1077  xZ0 = 0.01 *Z0[imo][ise][ila][ip];
1078  xZH = 0.01 *zH[imo][ise][ila][ip];
1079  xRUD= 0.01 *RUD[imo][ise][ila][ip];
1080 
1081  ytmp = xtU-f2p1-f2p2*pow((xaU-f2p0),f2p3)-f2p4*pow((xaU-f2p0),f2p5);
1082  ytmp2= xtD-ff2p1-ff2p2*pow((xaD-ff2p0),ff2p3)-ff2p4*pow((xaD-ff2p0),ff2p5);
1083 
1084  Zh = ((ytmp-ytmp2)*veff2*veff22 + Ze[0]*veff22 + Ze[1]*veff2)/(veff2 + veff22);
1085  pL = (Zh-xZ0)*(Zh-xZ0);
1086  pL =sqrt(pL + xRUD*xRUD);
1087 
1088  ytmp = ytmp-pL/cvel-abs(Zh-Ze[0])/veff2;
1089  ytmp2 = ytmp2-pL/cvel-abs(Zh-Ze[1])/veff22;
1090 
1091  y1[ipo2]=xtU-pL/cvel-abs(Zh-Ze[0])/veff2;
1092  y2[ipo2]=xtD-pL/cvel-abs(Zh-Ze[1])/veff22;
1093 
1094 
1095  if (y1[ipo2]>-20.&&y1[ipo2]<30.&&y2[ipo2]>-20.&&y2[ipo2]<30.&&Ene[imo][ise][ila][ip]>0.&&abs(ytmp)<5.*xkoe&&abs(ytmp2)<5.*xkoe
1096  &&Zh>Ze[0]&&Zh<Ze[1]) {
1097  x1[ipo2]=xaU;
1098 // ey1[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
1099  ey1[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
1100  ex1[ipo2]= 0.1;
1101  z1[ipo2] = Zh;
1102 
1103  x2[ipo2]=xaD;
1104 // ey2[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
1105  ey2[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
1106  ex2[ipo2]= 0.1;
1107 
1108  ipo2++;
1109  }
1110  }
1111 
1112  itt3 = 1.*ipo2;
1113  if ((itt3/itt2)<0.9) {
1114  xkoe = xkoe * 1.5;
1115  f2p0 = f1p0;
1116  f2p1 = f1p1;
1117  f2p2 = f1p2;
1118  f2p3 = f1p3;
1119  f2p4 = f1p4;
1120  f2p5 = f1p5;
1121  ff2p0 = ff1p0;
1122  ff2p1 = ff1p1;
1123  ff2p2 = ff1p2;
1124  ff2p3 = ff1p3;
1125  ff2p4 = ff1p4;
1126  ff2p5 = ff1p5;
1127  cout<<"Itt3 reset: xkoe = "<<xkoe<<endl;
1128  goto metka2;
1129  }
1130  cout<<endl<<"3rd Iteration: "<<ipo2<<" events in the plot."<<endl<<endl;
1131 
1132 //*******
1133  fgg10 = new TF1("fgg10",twfunc2,5.,2500.,6);
1134 
1135  fgg10->SetParameter(0,f2p0); // Pole parameter
1136  fgg10->SetParLimits(0,0.,4.9); // Pole parameter
1137 
1138  xp1=f2p1 ;
1139  fgg10->SetParameter(1,xp1); // Level
1140  fgg10->SetParLimits(1,xp1-1.,xp1+1.);
1141 
1142  xp2=f2p2;
1143  fgg10->SetParameter(2,xp2); // Amp1
1144  fgg10->SetParLimits(2,0.9*xp2,1.1*xp2);
1145 
1146  xp3=f2p3;
1147  fgg10->SetParameter(3,xp3); // Pow1
1148  fgg10->SetParLimits(3,xp3-0.1,xp3+0.1);
1149 
1150  xp4 = f2p4;
1151  fgg10->SetParameter(4,xp4); // Amp2
1152  fgg10->SetParLimits(4,0.9*xp4,1.1*xp4);
1153 
1154  xp5 = f2p5;
1155  fgg10->SetParameter(5,xp5); // Pow2
1156  fgg10->SetParLimits(5,xp5-0.1,xp5+0.1);
1157 
1158  fgg10->SetLineColor(4);
1159  fgg10->SetLineWidth(1.);
1160 
1161  gr01 = new TGraphErrors(ipo2,x1,y1,ex1,ey1);
1162  gr01->Fit("fgg10","MER");
1163  f3p0 = fgg10->GetParameter(0);
1164  f3p1 = fgg10->GetParameter(1);
1165  f3p2 = fgg10->GetParameter(2);
1166  f3p3 = fgg10->GetParameter(3);
1167  f3p4 = fgg10->GetParameter(4);
1168  f3p5 = fgg10->GetParameter(5);
1169 
1170  cout<<"-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-3-"<<endl;
1171 //********
1172  fgg20 = new TF1("fgg20",twfunc2,5.,2500.,6);
1173 
1174  fgg20->SetParameter(0,ff2p0); // Pole parameter
1175  fgg20->SetParLimits(0,0.,4.9); // Pole parameter
1176 
1177  xp1=ff2p1 ;
1178  fgg20->SetParameter(1,xp1); // Level
1179  fgg20->SetParLimits(1,xp1-1.,xp1+1.);
1180 
1181  xp2=ff2p2;
1182  fgg20->SetParameter(2,xp2); // Amp1
1183  fgg20->SetParLimits(2,0.9*xp2,1.1*xp2);
1184 
1185  xp3=ff2p3;
1186  fgg20->SetParameter(3,xp3); // Pow1
1187  fgg20->SetParLimits(3,xp3-0.1,xp3+0.1);
1188 
1189  xp4 = ff2p4;
1190  fgg20->SetParameter(4,xp4); // Amp2
1191  fgg20->SetParLimits(4,0.9*xp4,1.1*xp4);
1192 
1193  xp5 = ff2p5;
1194  fgg20->SetParameter(5,xp5); // Pow2
1195  fgg20->SetParLimits(5,xp5-0.1,xp5+0.1);
1196 
1197  fgg20->SetLineColor(4);
1198  fgg20->SetLineWidth(1.);
1199 
1200  gr11 = new TGraphErrors(ipo2,x2,y2,ex2,ey2);
1201  gr11->Fit("fgg20","MER");
1202  ff3p0 = fgg20->GetParameter(0);
1203  ff3p1 = fgg20->GetParameter(1);
1204  ff3p2 = fgg20->GetParameter(2);
1205  ff3p3 = fgg20->GetParameter(3);
1206  ff3p4 = fgg20->GetParameter(4);
1207  ff3p5 = fgg20->GetParameter(5);
1208 
1209 //********
1210 //********
1211 
1212  ipo3 = 0;
1213  for (int ii=0; ii<ipo2; ii++) {
1214  ytmp=y1[ii]-f3p1-f3p2*pow((x1[ii]-f3p0),f3p3)-f3p4*pow((x1[ii]-f3p0),f3p5);
1215  ytmp2=y2[ii]-ff3p1-ff3p2*pow((x2[ii]-ff3p0),ff3p3)-ff3p4*pow((x2[ii]-ff3p0),ff3p5);
1216 
1217  if (abs(ytmp)<2.5 && abs(ytmp2)<2.5) {
1218  yy1[ipo3]=ytmp;
1219  yy2[ipo3]=ytmp2;
1220  eyy1[ipo3]=ey1[ii];
1221  xx1[ipo3]=z1[ii];
1222  exx1[ipo3]=0.1;
1223  ipo3++;
1224  }
1225  }
1226 
1227  fz1 = new TF1("fz1",tzfunc, 10.,420.,2);
1228  fz1->SetLineColor(2);
1229  fz1->SetLineWidth(1.);
1230 
1231  fz1->SetParameter(0,0.); // Extra-Shift
1232  fz1->SetParLimits(0,-0.9,0.9);
1233  fz1->SetParameter(1,0.0); // Slope
1234  fz1->SetParLimits(1,-0.1,0.1);
1235 
1236  gr02 = new TGraphErrors(ipo3,xx1,yy1,exx1,eyy1);
1237  gr02->Fit("fz1","MER");
1238  zp0 = fz1->GetParameter(0);
1239  zp1 = fz1->GetParameter(1);
1240  veff2 = veff2/(1+veff2*zp1);
1241 
1242 //********
1243 
1244  fz2 = new TF1("fz2",tzfunc, 10.,420.,2);
1245  fz2->SetLineColor(2);
1246  fz2->SetLineWidth(1.);
1247 
1248  fz2->SetParameter(0,0.); // Extra-Shift
1249  fz2->SetParLimits(0,-0.9,0.9);
1250  fz2->SetParameter(1,0.0); // Slope
1251  fz2->SetParLimits(1,-0.02,0.02);
1252 
1253  gr12 = new TGraphErrors(ipo3,xx1,yy2,exx1,eyy1);
1254  gr12->Fit("fz2","MER");
1255  zzp0 = fz2->GetParameter(0);
1256  zzp1 = fz2->GetParameter(1);
1257  veff22 = veff22/(1-veff22*zzp1);
1258 
1259 delete gr01;
1260 delete gr02;
1261 delete gr11;
1262 delete gr12;
1263 delete fgg10;
1264 delete fgg20;
1265 delete fz1;
1266 delete fz2;
1267 
1268 }
1269 
1270  if (veff2<15.7) veff2=16.2;
1271  if (veff2>16.7) veff2=16.2;
1272  if (veff22<16.2) veff22=16.8;
1273  if (veff22>17.2) veff22=16.8;
1274 
1275  cout<<endl<<"3rd veff2="<<veff2<<" veff22="<<veff22<<endl<<endl;
1276 
1277  vef3 = vef2;
1278  vef33= veff22;
1279 
1280 //*********************************************************************************
1281 //*********************************************************************************
1282 // 4th Iteration
1283 //*********************************************************************************
1284 //*********************************************************************************
1285 {
1286 
1287  xkoe = 1.;
1288 metka3: ipo2 = 0;
1289 
1290  for (unsigned short int ip=0; ip<ipo[imo][ise][ila]; ip++) {
1291 
1292  xaU = 1. *aU[imo][ise][ila][ip];
1293  xaD = 1. *aD[imo][ise][ila][ip];
1294  xtU = 0.01 *(tU[imo][ise][ila][ip]);
1295  xtD = 0.01 *(tD[imo][ise][ila][ip]);
1296  xZ0 = 0.01 *Z0[imo][ise][ila][ip];
1297  xZH = 0.01 *zH[imo][ise][ila][ip];
1298  xRUD= 0.01 *RUD[imo][ise][ila][ip];
1299 
1300  ytmp = xtU-f3p1-f3p2*pow((xaU-f3p0),f3p3)-f3p4*pow((xaU-f3p0),f3p5);
1301  ytmp2= xtD-ff3p1-ff3p2*pow((xaD-ff3p0),ff3p3)-ff3p4*pow((xaD-ff3p0),ff3p5);
1302 
1303  Zh = ((ytmp-ytmp2)*veff2*veff22 + Ze[0]*veff22 + Ze[1]*veff2)/(veff2 + veff22);
1304  pL = (Zh-xZ0)*(Zh-xZ0);
1305  pL =sqrt(pL + xRUD*xRUD);
1306 
1307  ytmp = ytmp-pL/cvel-abs(Zh-Ze[0])/veff2;
1308  ytmp2 = ytmp2-pL/cvel-abs(Zh-Ze[1])/veff22;
1309 
1310  y1[ipo2]=xtU-pL/cvel-abs(Zh-Ze[0])/veff2;
1311  y2[ipo2]=xtD-pL/cvel-abs(Zh-Ze[1])/veff22;
1312 
1313 
1314  if (y1[ipo2]>-20.&&y1[ipo2]<30.&&y2[ipo2]>-20.&&y2[ipo2]<30.&&Ene[imo][ise][ila][ip]>0.&&abs(ytmp)<5.*xkoe&&abs(ytmp2)<5.*xkoe
1315  &&Zh>Ze[0]&&Zh<Ze[1]) {
1316  x1[ipo2]=xaU;
1317 // ey1[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
1318  ey1[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
1319  ex1[ipo2]= 0.1;
1320  z1[ipo2] = Zh;
1321 
1322  x2[ipo2]=xaD;
1323 // ey2[ipo2]=sqrt(0.2*0.2 + 0.09*0.09/Ene[ip]);
1324  ey2[ipo2]=0.1+0.2/sqrt(Ene[imo][ise][ila][ip]);
1325  ex2[ipo2]= 0.1;
1326 
1327  ipo2++;
1328  }
1329  }
1330 
1331  itt4 = 1.*ipo2;
1332  if ((itt4/itt3)<0.9) {
1333  xkoe = xkoe * 1.5;
1334  f3p0 = f2p0;
1335  f3p1 = f2p1;
1336  f3p2 = f2p2;
1337  f3p3 = f2p3;
1338  f3p4 = f2p4;
1339  f3p5 = f2p5;
1340  ff3p0 = ff2p0;
1341  ff3p1 = ff2p1;
1342  ff3p2 = ff2p2;
1343  ff3p3 = ff2p3;
1344  ff3p4 = ff2p4;
1345  ff3p5 = ff2p5;
1346  cout<<"Itt4 reset: xkoe = "<<xkoe<<endl;
1347  goto metka3;
1348  }
1349  cout<<endl<<"4th Iteration: "<<ipo2<<" events in the plot."<<endl<<endl;
1350  cont4 = ipo2;
1351 
1352 //*******
1353  fgg10 = new TF1("fgg10",twfunc2,5.,2500.,6);
1354 
1355  fgg10->SetParameter(0,f3p0); // Pole parameter
1356  fgg10->SetParLimits(0,0.,4.9); // Pole parameter
1357 
1358  xp1=f3p1 ;
1359  fgg10->SetParameter(1,xp1); // Level
1360  fgg10->SetParLimits(1,xp1-1.,xp1+1.);
1361 
1362  xp2=f3p2;
1363  fgg10->SetParameter(2,xp2); // Amp1
1364  fgg10->SetParLimits(2,0.9*xp2,1.1*xp2);
1365 
1366  xp3=f3p3;
1367  fgg10->SetParameter(3,xp3); // Pow1
1368  fgg10->SetParLimits(3,xp3-0.05,xp3+0.05);
1369 
1370  xp4 = f3p4;
1371  fgg10->SetParameter(4,xp4); // Amp2
1372  fgg10->SetParLimits(4,0.9*xp4,1.1*xp4);
1373 
1374  xp5 = f3p5;
1375  fgg10->SetParameter(5,xp5); // Pow2
1376  fgg10->SetParLimits(5,xp5-0.05,xp5+0.05);
1377 
1378  fgg10->SetLineColor(4);
1379  fgg10->SetLineWidth(1.);
1380 
1381  gr01 = new TGraphErrors(ipo2,x1,y1,ex1,ey1);
1382  gr01->Fit("fgg10","MER");
1383  f4p0 = fgg10->GetParameter(0);
1384  f4p1 = fgg10->GetParameter(1);
1385  f4p2 = fgg10->GetParameter(2);
1386  f4p3 = fgg10->GetParameter(3);
1387  f4p4 = fgg10->GetParameter(4);
1388  f4p5 = fgg10->GetParameter(5);
1389 
1390  cout<<"-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-4-"<<endl;
1391 //********
1392  fgg20 = new TF1("fgg20",twfunc2,5.,2500.,6);
1393 
1394  fgg20->SetParameter(0,ff3p0); // Pole parameter
1395  fgg20->SetParLimits(0,0.,4.9); // Pole parameter
1396 
1397  xp1=ff3p1 ;
1398  fgg20->SetParameter(1,xp1); // Level
1399  fgg20->SetParLimits(1,xp1-1.,xp1+1.);
1400 
1401  xp2=ff3p2;
1402  fgg20->SetParameter(2,xp2); // Amp1
1403  fgg20->SetParLimits(2,0.9*xp2,1.1*xp2);
1404 
1405  xp3=ff3p3;
1406  fgg20->SetParameter(3,xp3); // Pow1
1407  fgg20->SetParLimits(3,xp3-0.05,xp3+0.05);
1408 
1409  xp4 = ff3p4;
1410  fgg20->SetParameter(4,xp4); // Amp2
1411  fgg20->SetParLimits(4,0.9*xp4,1.1*xp4);
1412 
1413  xp5 = ff3p5;
1414  fgg20->SetParameter(5,xp5); // Pow2
1415  fgg20->SetParLimits(5,xp5-0.05,xp5+0.05);
1416 
1417  fgg20->SetLineColor(4);
1418  fgg20->SetLineWidth(1.);
1419 
1420  gr11 = new TGraphErrors(ipo2,x2,y2,ex2,ey2);
1421  gr11->Fit("fgg20","MER");
1422  ff4p0 = fgg20->GetParameter(0);
1423  ff4p1 = fgg20->GetParameter(1);
1424  ff4p2 = fgg20->GetParameter(2);
1425  ff4p3 = fgg20->GetParameter(3);
1426  ff4p4 = fgg20->GetParameter(4);
1427  ff4p5 = fgg20->GetParameter(5);
1428 
1429 //********
1430 //********
1431 
1432  ipo3 = 0;
1433  for (int ii=0; ii<ipo2; ii++) {
1434  ytmp=y1[ii]-f4p1-f4p2*pow((x1[ii]-f4p0),f4p3)-f4p4*pow((x1[ii]-f4p0),f4p5);
1435  ytmp2=y2[ii]-ff4p1-ff4p2*pow((x2[ii]-ff4p0),ff4p3)-ff4p4*pow((x2[ii]-ff4p0),ff4p5);
1436 
1437  if (abs(ytmp)<2.5 && abs(ytmp2)<2.5) {
1438  yy1[ipo3]=ytmp;
1439  yy2[ipo3]=ytmp2;
1440  eyy1[ipo3]=ey1[ii];
1441  xx1[ipo3]=z1[ii];
1442  exx1[ipo3]=0.1;
1443  ipo3++;
1444  }
1445  }
1446 
1447  fz1 = new TF1("fz1",tzfunc, 10.,420.,2);
1448  fz1->SetLineColor(2);
1449  fz1->SetLineWidth(1.);
1450 
1451  fz1->SetParameter(0,0.); // Extra-Shift
1452  fz1->SetParLimits(0,-0.9,0.9);
1453  fz1->SetParameter(1,0.0); // Slope
1454  fz1->SetParLimits(1,-0.1,0.1);
1455 
1456  gr02 = new TGraphErrors(ipo3,xx1,yy1,exx1,eyy1);
1457  gr02->Fit("fz1","MER");
1458  zp0 = fz1->GetParameter(0);
1459  zp1 = fz1->GetParameter(1);
1460  veff2 = veff2/(1+veff2*zp1);
1461 
1462 //********
1463 
1464  fz2 = new TF1("fz2",tzfunc, 10.,420.,2);
1465  fz2->SetLineColor(2);
1466  fz2->SetLineWidth(1.);
1467 
1468  fz2->SetParameter(0,0.); // Extra-Shift
1469  fz2->SetParLimits(0,-0.9,0.9);
1470  fz2->SetParameter(1,0.0); // Slope
1471  fz2->SetParLimits(1,-0.02,0.02);
1472 
1473  gr12 = new TGraphErrors(ipo3,xx1,yy2,exx1,eyy1);
1474  gr12->Fit("fz2","MER");
1475  zzp0 = fz2->GetParameter(0);
1476  zzp1 = fz2->GetParameter(1);
1477  veff22 = veff22/(1-veff22*zzp1);
1478 
1479  cout<<endl<<"4th veff2="<<veff2<<" veff22="<<veff22<<endl<<endl;
1480 
1481 //************************************************************************************
1482 
1483  for (int ii=0; ii<ipo2; ii++) {
1484 // y7[ii]=y1[ii]-zp0-zp1*(z1[ii]-212.)-fp1-fp2*pow((x1[ii]-fp0),fp3)-fp4*pow((x1[ii]-fp0),fp5);
1485  y10[ii]=y1[ii]-f4p1-f4p2*pow((x1[ii]-f4p0),f4p3)-f4p4*pow((x1[ii]-f4p0),f4p5);
1486  his4->Fill(y10[ii]);
1487  y20[ii]=y2[ii]-ff4p1-ff4p2*pow((x2[ii]-ff4p0),ff4p3)-ff4p4*pow((x2[ii]-ff4p0),ff4p5);
1488  his44->Fill(y20[ii]);
1489  }
1490 //*******************************
1491 
1492  fgau = new TF1("fgau",gfunc,-0.55,0.85,3);
1493  fgau->SetNpx(100);
1494  fgau->SetParameter(1,0.); // Mean parameter
1495  fgau->SetParameter(2,0.4); // Sigma parameter
1496  his4->Fit("fgau","","MERQ",-0.75,0.75);
1497  gme = fgau->GetParameter(1);
1498  gsi = fgau->GetParameter(2);
1499  egme = fgau->GetParError(1);
1500  egsi = fgau->GetParError(2);
1501 
1502  fgau2 = new TF1("fgau2",gfunc,-0.55,0.85,3);
1503  fgau2->SetNpx(100);
1504  fgau2->SetParameter(1,0.); // Mean parameter
1505  fgau2->SetParameter(2,0.4); // Sigma parameter
1506  his44->Fit("fgau2","","MERQ",-0.75,0.75);
1507  gme2 = fgau2->GetParameter(1);
1508  gsi2 = fgau2->GetParameter(2);
1509  egme2 = fgau2->GetParError(1);
1510  egsi2 = fgau2->GetParError(2);
1511 //************************************************************************************
1512  dat100<<imo+1<<" "<<ise+1<<" "<<ila+1<<" 0 "
1513  <<f4p0<<" "<<f4p1<<" "<<f4p2<<" "<<f4p3<<" "<<f4p4<<" "<<f4p5<<" "
1514  <<zp0<<" "<<zp1<<" "
1515  <<gme<<" "<<egme<<" "<<gsi<<" "<<egsi<<" "
1516  <<vef3<<" "<<veff2<<endl;
1517  dat100<<imo+1<<" "<<ise+1<<" "<<ila+1<<" 1 "
1518  <<ff4p0<<" "<<ff4p1<<" "<<ff4p2<<" "<<ff4p3<<" "<<ff4p4<<" "<<ff4p5<<" "
1519  <<zzp0<<" "<<zzp1<<" "
1520  <<gme2<<" "<<egme2<<" "<<gsi2<<" "<<egsi2<<" "
1521  <<vef33<<" "<<veff22<<endl;
1522 
1523  datveff<<imo+1<<" "<<ila+1<<" "<<ise+1<<" "<<veff2<<" "<<veff22<<endl;
1524 
1525  datwalk<<imo+1<<" "<<ise+1<<" "<<ila+1<<" 0 "
1526  <<f4p0<<" "<<f4p1+gme<<" "<<f4p2<<" "<<f4p3<<" "<<f4p4<<" "<<f4p5<<endl;
1527  datwalk<<imo+1<<" "<<ise+1<<" "<<ila+1<<" 1 "
1528  <<ff4p0<<" "<<ff4p1+gme2<<" "<<ff4p2<<" "<<ff4p3<<" "<<ff4p4<<" "<<ff4p5<<endl;
1529 
1530  if (abs(gme)>0.250) datwarn<<"Shifted Mean: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=0 "<<endl;
1531  if (abs(gme2)>0.250) datwarn<<"Shifted Mean: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=1 "<<endl;
1532 
1533  if (abs(gsi)>0.60) datwarn<<"Big Sigma: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=0 "<<endl;
1534  if (abs(gsi2)>0.60) datwarn<<"Big Sigma: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=1 "<<endl;
1535 
1536  if (abs(gsi)<0.10) datwarn<<"Small Sigma: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=0 "<<endl;
1537  if (abs(gsi2)<0.10) datwarn<<"Small Sigma: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<" E=1 "<<endl;
1538 
1539  if (cont4<0.5*cont1) datwarn<<"Lost Statistics: M="<<imo+1<<" S="<<ise+1<<" L="<<ila+1<<endl;
1540 
1541 //************************************************************************************
1542  c1->Clear();
1543  c1->Divide(2,2);
1544  gStyle->SetOptFit(1110);
1545  c1->Range(0.1875,-81.7688,3.3125,735.919);
1546  c1->SetBorderSize(2);
1547  c1->SetFrameFillColor(0);
1548  c1->SetGridx();
1549  c1->SetGridy();
1550  c1->SetTitle("Datafile: ");
1551 // gStyle->SetPadTopMargin(.5);
1552 // gStyle->SetPadLeftMargin(2.5);
1553 // gStyle->SetPadRightMargin(2.5);
1554  gStyle->SetCanvasBorderMode(0);
1555  gStyle->SetOptFit(220);
1556  gStyle->SetOptStat(2220);
1557  gStyle->SetPalette(1);
1558 // gStyle->SetStatFontSize(0.5);
1559 
1560 // gStyle->SetStatW(0.35);
1561 // gStyle->SetStatH(0.10);
1562  gStyle->SetStatW(0.25);
1563  gStyle->SetStatH(0.12);
1564  gStyle->SetStatX(0.93);
1565  gStyle->SetStatY(0.92);
1566 
1567  gStyle->SetTitleW(0.9);
1568  gStyle->SetTitleH(0.08);
1569 
1570 // Set a bunch of parameters to make the plot look nice
1571 
1572 /* canvas->SetFillColor(0);
1573  canvas->UseCurrentStyle();
1574  canvas->SetBorderMode(0); // still leaves red frame bottom and right
1575  canvas->SetFrameBorderMode(0); // need this to turn off red hist frame!
1576  gROOT->SetStyle("Plain");
1577  canvas->UseCurrentStyle();
1578  gROOT->ForceStyle(); */
1579 
1580 // gStyle->SetOptStat(0);
1581 // gStyle->SetTitleBorderSize(0);
1582  gStyle->SetTitleSize(0.6,"t");
1583 // gStyle->SetTitleFontSize(0.3);
1584 // gStyle->SetTitleFont(42, "hxy"); // for histogram and axis titles
1585 // gStyle->SetLabelFont(42, "xyz"); // for axis labels (values)
1586 
1587 
1588 
1589 
1590  c1->Update();
1591  c1->cd(1);
1592  gPad->SetGridy(0);
1593  gPad->SetGridx(0);
1594  gPad->SetFillColor(0);
1595  gPad->SetTicks();
1596  gPad->SetTopMargin(0.13);
1597  gPad->SetBottomMargin(0.17);
1598  gPad->SetLeftMargin(0.1);
1599  gPad->SetRightMargin(0.1);
1600  gPad->SetLogx(1);
1601 
1602  gr1 = new TGraph(ipo2,x1,y1);
1603  gr1->SetTitle(hisname);
1604  gr1->GetXaxis()->SetTitle("Pulse Peak");
1605  gr1->GetXaxis()->SetLimits(4.5,3000.0);
1606  gr1->SetMinimum(-15.0);
1607  gr1->SetMaximum(25.0);
1608  gr1->GetYaxis()->SetNdivisions(810);
1609  gr1->GetYaxis()->SetTitle("Uncorrected #DeltaTime (ns)");
1610  gr1->SetMarkerColor(2);
1611 // gr1->SetLineColor(2);
1612  gr1->SetMarkerStyle(20);
1613  gr1->SetMarkerSize(0.3);
1614  gr1->Draw("AP");
1615 // fgg5->Draw("Same");
1616  fgg10->Draw("Same");
1617 //hl1->Draw();
1618 
1619 
1620  ftxt = new char[90];
1621  sprintf(ftxt, "f(X) = %5.2f + %4.2f#upoint(X-%4.2f)^{%5.2f}",f4p1,f4p2,f4p0,f4p3);
1622  ftxt2 = new char[90];
1623  sprintf(ftxt2, " + %4.2f#upoint(X-%4.2f)^{%5.2f}",f4p4,f4p0,f4p5);
1624  txt1 = new TLatex(40.,21.,ftxt);
1625  txt1->SetTextFont(52);
1626  txt1->SetTextColor(4);
1627  txt1->Draw("Same");
1628  txt2 = new TLatex(70.,17.0,ftxt2);
1629  txt2->SetTextFont(52);
1630  txt2->SetTextColor(4);
1631  txt2->Draw("Same");
1632 
1633  c1->cd(2);
1634  gPad->SetGridy(0);
1635  gPad->SetGridx(0);
1636  gPad->SetFillColor(0);
1637  gPad->SetTicks();
1638  gPad->SetTopMargin(0.13);
1639  gPad->SetBottomMargin(0.17);
1640  gPad->SetLeftMargin(0.1);
1641  gPad->SetRightMargin(0.1);
1642  gPad->SetLogx(1);
1643 
1644  gr3 = new TGraph(ipo2,x1,y10);
1645  gr3->SetTitle(hisname);
1646  gr3->GetXaxis()->SetTitle("Pulse Peak");
1647  gr3->GetXaxis()->SetLimits(4.5,3000.);
1648  gr3->SetMinimum(-10.0);
1649  gr3->SetMaximum(14.0);
1650  gr3->GetYaxis()->SetNdivisions(810);
1651  gr3->GetYaxis()->SetTitle("Corrected #DeltaTime (ns)");
1652  gr3->SetMarkerColor(4);
1653 // gr3->SetLineColor(2);
1654  gr3->SetMarkerStyle(20);
1655  gr3->SetMarkerSize(0.3);
1656  gr3->Draw("AP");
1657 
1658  lin1 = new TLine(5.,0.,2000.,0.);
1659  lin1->SetLineColor(2);
1660  lin1->Draw();
1661 // his1->Draw();
1662 //hl2->Draw();
1663 
1664  c1->cd(3);
1665  gPad->SetGridy(0);
1666  gPad->SetGridx(0);
1667  gPad->SetFillColor(0);
1668  gPad->SetTicks();
1669  gPad->SetTopMargin(0.13);
1670  gPad->SetBottomMargin(0.17);
1671  gPad->SetLeftMargin(0.1);
1672  gPad->SetRightMargin(0.1);
1673  gPad->SetLogx(0);
1674 
1675  gr2 = new TGraph(ipo3,xx1,yy1);
1676  gr2->SetTitle(hisname);
1677  gr2->GetXaxis()->SetTitle("Z-Coordinate (cm)");
1678 // gr2->GetXaxis()->SetLimits(-5.0,5.0);
1679  gr2->SetMinimum(-5.0);
1680  gr2->SetMaximum(5.0);
1681  gr2->GetYaxis()->SetNdivisions(810);
1682  gr2->GetYaxis()->SetTitle("Corrected #DeltaTime (ns)");
1683  gr2->SetMarkerColor(4);
1684 // gr2->SetLineColor(2);
1685  gr2->SetMarkerStyle(20);
1686  gr2->SetMarkerSize(0.3);
1687  gr2->Draw("AP");
1688  fz1->Draw("Same");
1689 //hl3->Draw();
1690 
1691  ftxt3 = new char[90];
1692  sprintf(ftxt3, "Optimal v_{eff} = %5.2f cm/ns",veff2);
1693  txt3 = new TLatex(30.,3.0,ftxt3);
1694  txt3->SetTextFont(52);
1695  txt3->SetTextColor(2);
1696  txt3->Draw("Same");
1697 
1698  ftxt4 = new char[90];
1699  sprintf(ftxt4, "Initial v_{eff} = %5.2f cm/ns",vef1);
1700  txt4 = new TLatex(30.,4.0,ftxt4);
1701  txt4->SetTextFont(52);
1702  txt4->SetTextColor(4);
1703  txt4->Draw("Same");
1704 
1705  c1->cd(4);
1706 // his4->Fit("gaus","","MER",-1.,0.8);
1707  his4->Draw();
1708 // gPad->SetLogy(1);
1709 //hl0->Draw();
1710 
1711  c1->Print(pdfname);
1712 
1713 //***************************************************
1714  c1->Update();
1715 
1716  c1->cd(1);
1717  gPad->SetGridy(0);
1718  gPad->SetGridx(0);
1719  gPad->SetFillColor(0);
1720  gPad->SetTicks();
1721  gPad->SetTopMargin(0.13);
1722  gPad->SetBottomMargin(0.17);
1723  gPad->SetLeftMargin(0.1);
1724  gPad->SetRightMargin(0.1);
1725  gPad->SetLogx(1);
1726 
1727  gr1 = new TGraph(ipo2,x2,y2);
1728  gr1->SetTitle(hisname2);
1729  gr1->GetXaxis()->SetTitle("Pulse Peak");
1730  gr1->GetXaxis()->SetLimits(4.5,3000.0);
1731  gr1->SetMinimum(-15.0);
1732  gr1->SetMaximum(25.0);
1733  gr1->GetYaxis()->SetNdivisions(810);
1734  gr1->GetYaxis()->SetTitle("Uncorrected #DeltaTime (ns)");
1735  gr1->SetMarkerColor(2);
1736 // gr1->SetLineColor(2);
1737  gr1->SetMarkerStyle(20);
1738  gr1->SetMarkerSize(0.3);
1739  gr1->Draw("AP");
1740 // fgg5->Draw("Same");
1741  fgg20->Draw("Same");
1742 //hl1->Draw();
1743 
1744 
1745  ftxt = new char[90];
1746  sprintf(ftxt, "f(X) = %5.2f + %4.2f#upoint(X-%4.2f)^{%5.2f}",ff4p1,ff4p2,ff4p0,ff4p3);
1747  ftxt2 = new char[90];
1748  sprintf(ftxt2, " + %4.2f#upoint(X-%4.2f)^{%5.2f}",ff4p4,ff4p0,ff4p5);
1749  txt1 = new TLatex(40.,21.,ftxt);
1750  txt1->SetTextFont(52);
1751  txt1->SetTextColor(4);
1752  txt1->Draw("Same");
1753  txt2 = new TLatex(70.,17.,ftxt2);
1754  txt2->SetTextFont(52);
1755  txt2->SetTextColor(4);
1756  txt2->Draw("Same");
1757 
1758  c1->cd(2);
1759  gPad->SetGridy(0);
1760  gPad->SetGridx(0);
1761  gPad->SetFillColor(0);
1762  gPad->SetTicks();
1763  gPad->SetTopMargin(0.13);
1764  gPad->SetBottomMargin(0.17);
1765  gPad->SetLeftMargin(0.1);
1766  gPad->SetRightMargin(0.1);
1767  gPad->SetLogx(1);
1768 
1769  TGraph *gr3 = new TGraph(ipo2,x2,y20);
1770  gr3->SetTitle(hisname2);
1771  gr3->GetXaxis()->SetTitle("Pulse Peak");
1772  gr3->GetXaxis()->SetLimits(4.5,3000.0);
1773  gr3->SetMinimum(-10.0);
1774  gr3->SetMaximum(14.0);
1775  gr3->GetYaxis()->SetNdivisions(810);
1776  gr3->GetYaxis()->SetTitle("Corrected #DeltaTime (ns)");
1777  gr3->SetMarkerColor(4);
1778 // gr3->SetLineColor(2);
1779  gr3->SetMarkerStyle(20);
1780  gr3->SetMarkerSize(0.3);
1781  gr3->Draw("AP");
1782 
1783  lin1 = new TLine(5.,0.,2000.,0.);
1784  lin1->SetLineColor(2);
1785  lin1->Draw();
1786 
1787  c1->cd(3);
1788  gPad->SetGridy(0);
1789  gPad->SetGridx(0);
1790  gPad->SetFillColor(0);
1791  gPad->SetTicks();
1792  gPad->SetTopMargin(0.13);
1793  gPad->SetBottomMargin(0.17);
1794  gPad->SetLeftMargin(0.1);
1795  gPad->SetRightMargin(0.1);
1796  gPad->SetLogx(0);
1797 
1798  gr2 = new TGraph(ipo3,xx1,yy2);
1799  gr2->SetTitle(hisname2);
1800  gr2->GetXaxis()->SetTitle("Z-Coordinate (cm)");
1801 // gr2->GetXaxis()->SetLimits(-5.0,5.0);
1802  gr2->SetMinimum(-5.0);
1803  gr2->SetMaximum(5.0);
1804  gr2->GetYaxis()->SetNdivisions(810);
1805  gr2->GetYaxis()->SetTitle("Corrected #DeltaTime (ns)");
1806  gr2->SetMarkerColor(4);
1807 // gr2->SetLineColor(2);
1808  gr2->SetMarkerStyle(20);
1809  gr2->SetMarkerSize(0.3);
1810  gr2->Draw("AP");
1811  fz1->Draw("Same");
1812 //hl3->Draw();
1813 
1814  ftxt3 = new char[90];
1815  sprintf(ftxt3, "Optimal v_{eff} = %5.2f cm/ns",veff22);
1816  txt3 = new TLatex(30.,3.0,ftxt3);
1817  txt3->SetTextFont(52);
1818  txt3->SetTextColor(2);
1819  txt3->Draw("Same");
1820 
1821  ftxt4 = new char[90];
1822  sprintf(ftxt4, "Initial v_{eff} = %5.2f cm/ns",vef2);
1823  txt4 = new TLatex(30.,4.0,ftxt4);
1824  txt4->SetTextFont(52);
1825  txt4->SetTextColor(4);
1826  txt4->Draw("Same");
1827 
1828  c1->cd(4);
1829 // his4->Fit("gaus","","MER",-1.,0.8);
1830  his44->Draw();
1831 // gPad->SetLogy(1);
1832 //hl0->Draw();
1833 
1834  c1->Print(pdfname);
1835 
1836 delete gr01;
1837 delete gr02;
1838 delete gr11;
1839 delete gr12;
1840 delete fgg10;
1841 delete fgg20;
1842 delete fz1;
1843 delete fz2;
1844 
1845 
1846 delete fgau;
1847 delete fgau2;
1848 delete h1;
1849 delete h2;
1850 delete h3;
1851 delete h4;
1852 delete his4;
1853 delete his44;
1854 delete gr1;
1855 delete gr2;
1856 delete gr3;
1857 delete ftxt;
1858 delete ftxt2;
1859 delete ftxt3;
1860 delete ftxt4;
1861 delete txt1;
1862 delete txt2;
1863 delete txt3;
1864 delete txt4;
1865 delete hisname;
1866 delete hisname2;
1867 delete lin1;
1868 
1869 } //End of 4th Iteration
1870 
1871 } //End of ise loop
1872 } //End of ila loop
1873 } //End of imo loop
1874 
1875 //*************************************************************
1876  c1->Update();
1877  c1->Print(pdfname2);
1878  dat100.close();
1879  datveff.close();
1880  datwalk.close();
1881  datwarn.close();
1882 //************************************************************************************
1883 
1884  return NOERROR;
1885 }
Ze
static double Ze[2]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:96
twfunc2
static Double_t twfunc2(Double_t *x, Double_t *par)
Definition: JEventProcessor_BCAL_TimeCalibration.cc:150
DEventRFBunch.h
ftxt2
static char * ftxt2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:124
E1
static double E1[100]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:138
Df250PulsePedestal.h
aD
static unsigned short int aD[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:89
DBCALTDCHit.h
DTAGMHit.h
DBCALShower.h
fgau2
static TF1 * fgau2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:119
DTrackTimeBased.h
DDetectorMatches.h
DVector3
TVector3 DVector3
Definition: DVector3.h:14
JEventProcessor_BCAL_TimeCalibration.h
DL1Trigger::trig_mask
uint32_t trig_mask
Definition: DL1Trigger.h:18
DBCALHit::module
int module
Definition: DBCALHit.h:25
x
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
Z0
static unsigned short int Z0[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:90
sprintf
sprintf(text,"Post KinFit Cut")
y
#define y
txt2
static TLatex * txt2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:131
gr3
static TGraph * gr3
Definition: JEventProcessor_BCAL_TimeCalibration.cc:107
JEventProcessor_BCAL_TimeCalibration::erun
jerror_t erun(void)
Called everytime run number changes, provided brun has been called.
Definition: JEventProcessor_BCAL_TimeCalibration.cc:503
ftxt4
static char * ftxt4
Definition: JEventProcessor_BCAL_TimeCalibration.cc:126
DL1Trigger.h
DApplication.h
dZs1
static double dZs1[3]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:140
fgau
static TF1 * fgau
Definition: JEventProcessor_BCAL_TimeCalibration.cc:118
DBCALUnifiedHit.h
ievent
static int ievent
Definition: JEventProcessor_BCAL_TimeCalibration.cc:78
DBCALHit::layer
int layer
Definition: DBCALHit.h:26
txt4
static TLatex * txt4
Definition: JEventProcessor_BCAL_TimeCalibration.cc:133
cvel
static double cvel
Definition: JEventProcessor_BCAL_TimeCalibration.cc:87
fz1
static TF1 * fz1
Definition: JEventProcessor_BCAL_TimeCalibration.cc:120
DBCALTruthShower.h
SYS_BCAL
Definition: GlueX.h:19
c1
Double_t c1[2][NMODULES]
Definition: tw_corr.C:68
gr02
static TGraphErrors * gr02
Definition: JEventProcessor_BCAL_TimeCalibration.cc:102
gr12
static TGraphErrors * gr12
Definition: JEventProcessor_BCAL_TimeCalibration.cc:103
gr01
static TGraphErrors * gr01
Definition: JEventProcessor_BCAL_TimeCalibration.cc:100
DBCALHit::end
DBCALGeometry::End end
Definition: DBCALHit.h:28
h4
static TH1F * h4
Definition: JEventProcessor_BCAL_TimeCalibration.cc:112
gr2
static TGraph * gr2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:106
r2U
static double r2U
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
PHIp
static double PHIp
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
E3
static double E3[100]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:138
fgg10
static TF1 * fgg10
Definition: JEventProcessor_BCAL_TimeCalibration.cc:116
hisname2
static char * hisname2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:128
JEventProcessor_BCAL_TimeCalibration::init
jerror_t init(void)
Called once at program start.
Definition: JEventProcessor_BCAL_TimeCalibration.cc:194
DBCALPoint.h
iRunNumber
static int iRunNumber
Definition: JEventProcessor_BCAL_TimeCalibration.cc:80
h3
static TH1F * h3
Definition: JEventProcessor_BCAL_TimeCalibration.cc:111
tzfunc
static Double_t tzfunc(Double_t *x, Double_t *par)
Definition: JEventProcessor_BCAL_TimeCalibration.cc:144
txt3
static TLatex * txt3
Definition: JEventProcessor_BCAL_TimeCalibration.cc:132
his44
static TH1F * his44
Definition: JEventProcessor_BCAL_TimeCalibration.cc:114
InitPlugin
InitPlugin_t InitPlugin
aU
static unsigned short int aU[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:89
hisname
static char * hisname
Definition: JEventProcessor_BCAL_TimeCalibration.cc:127
DTAGHHit.h
PHIdiff
static double PHIdiff
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
tADC_U
static double tADC_U
Definition: JEventProcessor_BCAL_TimeCalibration.cc:137
pdfps
static int pdfps
Definition: JEventProcessor_BCAL_TimeCalibration.cc:98
DBCALHit.h
r2D
static double r2D
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
DBCALGeometry.h
h2
static TH1F * h2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:110
cosPHI
static double cosPHI
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
fgg20
static TF1 * fgg20
Definition: JEventProcessor_BCAL_TimeCalibration.cc:117
gr1
static TGraph * gr1
Definition: JEventProcessor_BCAL_TimeCalibration.cc:105
mtx
static mutex mtx
Definition: JEventProcessor_BCAL_TimeCalibration.cc:76
JEventProcessor_BCAL_TimeCalibration::~JEventProcessor_BCAL_TimeCalibration
~JEventProcessor_BCAL_TimeCalibration()
Definition: JEventProcessor_BCAL_TimeCalibration.cc:186
Df250PulseIntegral.h
JEventProcessor_BCAL_TimeCalibration::brun
jerror_t brun(jana::JEventLoop *eventLoop, int32_t runnumber)
Called everytime a new run number is detected.
Definition: JEventProcessor_BCAL_TimeCalibration.cc:221
DBCALHit::sector
int sector
Definition: DBCALHit.h:27
ftxt
static char * ftxt
Definition: JEventProcessor_BCAL_TimeCalibration.cc:123
dZs2
static double dZs2[3]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:141
zH
static unsigned short int zH[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:90
E13
static double E13[100]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:138
sqrt
double sqrt(double)
sin
double sin(double)
h1
static TH1F * h1
Definition: JEventProcessor_BCAL_TimeCalibration.cc:109
DFCALShower.h
his4
static TH1F * his4
Definition: JEventProcessor_BCAL_TimeCalibration.cc:113
JEventProcessor_BCAL_TimeCalibration::evnt
jerror_t evnt(jana::JEventLoop *eventLoop, uint64_t eventnumber)
Called every event.
Definition: JEventProcessor_BCAL_TimeCalibration.cc:231
gr11
static TGraphErrors * gr11
Definition: JEventProcessor_BCAL_TimeCalibration.cc:101
busyShower
static int busyShower[1000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:85
fz2
static TF1 * fz2
Definition: JEventProcessor_BCAL_TimeCalibration.cc:121
DAnalysisUtilities.h
ftxt3
static char * ftxt3
Definition: JEventProcessor_BCAL_TimeCalibration.cc:125
Ene
static float Ene[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:91
JEventProcessor_BCAL_TimeCalibration::fini
jerror_t fini(void)
Called after last event of last event source has been processed.
Definition: JEventProcessor_BCAL_TimeCalibration.cc:515
ip
static int ip
Definition: JEventProcessor_BCAL_TimeCalibration.cc:97
DBCALHit::pulse_peak
int pulse_peak
Definition: DBCALHit.h:29
func
double func(double *x, double *par)
Definition: timewalk_fits.C:16
imo
static int imo
Definition: JEventProcessor_BCAL_TimeCalibration.cc:97
JEventProcessor_BCAL_TimeCalibration
Definition: JEventProcessor_BCAL_TimeCalibration.h:27
ila
static int ila
Definition: JEventProcessor_BCAL_TimeCalibration.cc:97
RUD
static unsigned short int RUD[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:90
DVertex.h
txt1
static TLatex * txt1
Definition: JEventProcessor_BCAL_TimeCalibration.cc:130
lin1
static TLine * lin1
Definition: JEventProcessor_BCAL_TimeCalibration.cc:135
tD
static unsigned short int tD[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:89
ipo
static unsigned short int ipo[48][4][3]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:94
tADC_D
static double tADC_D
Definition: JEventProcessor_BCAL_TimeCalibration.cc:137
zUD
static double zUD
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
itD
static double itD
Definition: JEventProcessor_BCAL_TimeCalibration.cc:95
PHIs
static double PHIs
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
DMCThrown.h
gfunc
static Double_t gfunc(Double_t *x, Double_t *par)
Definition: JEventProcessor_BCAL_TimeCalibration.cc:158
tU
static unsigned short int tU[48][4][3][5000]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:89
R0
static double R0[5]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:82
DChargedTrack.h
DD
static double DD[5]
Definition: JEventProcessor_BCAL_TimeCalibration.cc:83
rUD
static double rUD
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
DBCALDigiHit.h
itU
static double itU
Definition: JEventProcessor_BCAL_TimeCalibration.cc:95
JEventProcessor_BCAL_TimeCalibration::JEventProcessor_BCAL_TimeCalibration
JEventProcessor_BCAL_TimeCalibration()
Definition: JEventProcessor_BCAL_TimeCalibration.cc:178
sinPHI
static double sinPHI
Definition: JEventProcessor_BCAL_TimeCalibration.cc:93
ise
static int ise
Definition: JEventProcessor_BCAL_TimeCalibration.cc:97
DBCALCluster.h