StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
LanaTrees.C
1 /*
2  root.exe 'bfc.C(-1,"lana,nodefault")' LanaTrees.C+
3  */
4 // $Id: LanaTrees.C,v 1.4 2015/02/13 22:55:17 fisyak Exp $
5 // $Log: LanaTrees.C,v $
6 // Revision 1.4 2015/02/13 22:55:17 fisyak
7 // Force East == West drift velocities for Run XIV, because of bad East one
8 //
9 // Revision 1.3 2015/02/13 15:31:24 fisyak
10 // Check fit status
11 //
12 // Revision 1.2 2015/02/10 20:27:16 fisyak
13 // Adjust split style for ROOT_VERSION_CODE
14 //
15 // Revision 1.1 2015/01/20 19:41:54 fisyak
16 // Add new Laser drift velocity calculator
17 //
18 // Revision 1.1 2015/01/09 16:56:30 fisyak
19 // Freeze
20 //
21 // Revision 1.15 2014/03/19 20:52:37 fisyak
22 // Tight good velocity cut
23 //
24 // Revision 1.14 2014/03/13 21:58:43 fisyak
25 // Add more plots
26 //
27 // Revision 1.13 2013/05/28 20:08:49 fisyak
28 // Increase cut on membrane drift distance precision 1e-3 => 2e-3
29 //
30 // Revision 1.12 2013/05/28 06:59:56 fisyak
31 // Add drift velocity estimation from membrane cluster positions if standard procedure fails
32 //
33 // Revision 1.11 2013/05/20 12:35:12 fisyak
34 // Set time stamp at the begin of the laser run
35 //
36 // Revision 1.10 2013/05/16 15:25:39 fisyak
37 // Relax cut on drift velocity precision : 5e-4 => 1e-3
38 //
39 // Revision 1.9 2010/01/02 23:41:19 genevb
40 // Fix issues with dates starting with 2010
41 //
42 // Revision 1.8 2010/01/02 19:29:35 genevb
43 // switch to laser.root files as default
44 //
45 // Revision 1.7 2009/03/11 14:39:33 fisyak
46 // relax cuts on no. of good lasers and drift velocity precision
47 //
48 // Revision 1.6 2009/03/06 22:46:43 fisyak
49 // Increase acceptable drift velocity interval from [5.5,5.9] to [5.2,5.9]
50 //
51 // Revision 1.2 2008/04/25 15:25:15 fisyak
52 // Freeze macros
53 //
54 // Revision 1.4 2007/12/28 13:20:25 fisyak
55 // Use average drift velocity from East and West
56 //
57 // Revision 1.3 2007/12/10 19:54:02 fisyak
58 // Add Id and Log, correct spelling error in README
59 //
60 //#define ADJUSTABLE_BINNING
61 //#define __REFIT__
62 //#define INTEGRATE_OVER_HOURS
63 //#define SeparateWestandEast
64 //#define __Memberane__
65 #if !defined(__CINT__) || defined(__MAKECINT__)
66 //#include <ostream>
67 #include "Riostream.h"
68 #include <stdio.h>
69 #include "TROOT.h"
70 #include "TSystem.h"
71 #include "TMath.h"
72 #include "TH1.h"
73 #include "TH2.h"
74 #include "TH3.h"
75 #include "TAxis.h"
76 #include "TStyle.h"
77 #include "TF1.h"
78 #include "TFitResult.h"
79 #include "TProfile.h"
80 #include "TTree.h"
81 #include "TChain.h"
82 #include "TFile.h"
83 #include "TNtuple.h"
84 #include "TCanvas.h"
85 #include "TString.h"
86 #include "TList.h"
87 #include "TLegend.h"
88 #include "TFileSet.h"
89 #include "StLaserAnalysisMaker/LaserEvent.h"
90 #include "TDataSetIter.h"
91 #include "TDirIter.h"
92 #include "TTreeIter.h"
93 #endif
94 static Int_t run = 0;
95 static Int_t date = 0;
96 static Int_t Time = 0;
97 // Delta/dv All,West,East
98 static Double_t DVAll[2][3];
99 static Double_t dDVAll[2][3];
100 static Int_t _debug = 0;
101 static Double_t sigmaAcceptedDV = 5e-4; // maximum sigma for acceptable drift velocity
102 TH2D *dv = 0;
103 TH2D *slope = 0;
104 TH2D *memAdc = 0;
105 // io
106 TH2D *zMembrane[2];
107 // io we xy
108 TH2D *dMembraneY[2][2][2];
109 TNtuple *runNT = 0;
110 struct Run_t {
111  Float_t run, date, time, events, day, dvAll, ddvAll, dvWest, ddvWest, dvEast, ddvEast, slAll, dslAll, slWest, dslWest, slEast, dslEast;
112  Float_t vWest, vEast, zWI, dzWI, zEI, dzEI, zWO, dzWO, zEO, dzEO, utime, ok, dvSet, XWI, dXWI, XEI, dXEI, XWO, dXWO, XEO, dXEO, YWI, dYWI, YEI, dYEI, YWO, dYWO, YEO, dYEO;
113 };
114 const Char_t *vRun = "run:date:time:events:day:dvAll:ddvAll:dvWest:ddvWest:dvEast:ddvEast:slAll:dslAll:slWest:dslWest:slEast:dslEast:vWest:vEast:zWI:dzWI:zEI:dzEI:zWO:dzWO:zEO:dzEO:utime:ok:dvSet:XWI:dXWI:XEI:dXEI:XWO:dXWO:XEO:dXEO:YWI:dYWI:YEI:dYEI:YWO:dYWO:YEO:dYEO";
115 Run_t Run;
116 //________________________________________________________________________________
117 Double_t memDelay(Double_t *x, Double_t *p) {
118  static Double_t MembraneDelay[2][24] = {// delay in cm
119 // {397.377,366.758,341.826,332.002,341.826,366.758,397.377,425.086,443.949,450.599,443.949,425.086,448.749,455.399,448.749,429.886,402.177,371.558,346.626,336.802,346.626,371.558,402.177,429.886}, //Inner
120 // {423.378,372.639,328.529,310.094,328.529,372.639,423.378,467.23 ,496.294,506.418,496.294,467.23 ,501.094,511.218,501.094,461.02 ,416.682,377.439,333.329,314.894,333.329,377.439,428.178,472.03 }};//Outer
121  {398.122,367.481,342.294,332.015,341.383,366.056,396.65,424.512,443.64 ,450.607,444.276,425.678,448.44 ,455.407,449.076,430.478,402.922,372.281,347.094,336.815,346.183,370.856,401.45,429.312}, //Inner
122  {424.573,373.866,329.379,310.108,327.703,371.43 ,422.198,466.33,495.816,506.426,496.787,468.144,500.616,511.226,501.587,472.944,429.373,378.666,334.179,314.908,332.503,376.23,426.998,471.13}}; //Outer
123  Int_t sector = TMath::Nint(x[0]);
124  Double_t value = 0;
125  if (sector < 1 || sector > 24) return value;
126  Int_t inout = TMath::Nint(p[0]); // in = 0, out = 1
127  if (inout < 0 || inout > 1) return value;
128  return p[1] + p[2]*MembraneDelay[inout][sector-1];
129 }
130 //________________________________________________________________________________
131 TF1 *MemDelay(Int_t inout = 0) {
132  TF1 *f = new TF1("MemDelay",memDelay,-0.5,24.5,3);
133  f->SetParameters(inout,0,1);
134  f->SetParNames("InOut","Shift","Scale");
135  f->FixParameter(0,inout);
136  return f;
137 }
138 //________________________________________________________________________________
139 Double_t ScaleE2W(Double_t day) {// scale East to West drift velocity
140  // RunNT->Draw("1e3*(dvEast/dvWest-1):day>>diffE(30,90,180)","(ddvWest>0&&ddvWest<4e-5&&ddvEast>0&&ddvEast<4e-5)/((ddvWest/dvWest)**2+(ddvEast/dvEast)**2)","profw")
141  static Double_t par[2] = {1.57361e-01,-4.59752e-03};
142  return par[0] + par[1]*day;
143 }
144 //________________________________________________________________________________
145 void MakeTable() {
146  Double_t dv = DVAll[0][0];
147  Double_t ddv = dDVAll[0][0];
148  Double_t dvWest = DVAll[0][1];
149  Double_t ddvWest = dDVAll[0][1];
150  Double_t dvEast = DVAll[0][2];
151  Double_t ddvEast = dDVAll[0][2];
152  Run.ok = 0; // ok == 0 => use both: west and east; ok == 1 => use averaged drift velocities; ok > 1 ==> no. acceptable drift velocities
153  if (dvWest < 5.3 || dvWest > 5.9 || ddvWest <= 0 || ddvWest>sigmaAcceptedDV ||
154  dvEast < 5.3 || dvEast > 5.9 || ddvEast <= 0 || ddvEast>sigmaAcceptedDV) {
155  // if (! (dvWest < 5.3 && dvWest > 5.9 && ddvWest < 0 && ddvWest>sigmaAcceptedDV) ||
156  // ! (dvEast < 5.3 && dvEast > 5.9 && ddvEast < 0 && ddvEast>sigmaAcceptedDV)) {
157  cout << "Run " << run << " fails ============================= to make separated East and West drift velocities" << endl;
158  cout << "vWest = " << dvWest << " +/- " << ddvWest
159  << "\tvEast = " << dvEast << " +/- " << ddvEast << endl;
160  Run.ok = 1;
161  }
162 #ifndef SeparateWestandEast
163  Run.ok = 1;
164 #endif
165  // if (ok == 1 && ! (dv > 5.3 && dv < 5.9 && ddv > 0 && ddv<sigmaAcceptedDV)) {
166  if (Run.ok == 1 && (dv < 5.3 || dv > 5.9 || ddv <= 0 || ddv >sigmaAcceptedDV)) {
167  cout << "Run " << run << " fails ============================= to make averaged drift velocities" << endl;
168  cout << "v = " << dv << " +/- " << ddv << endl;
169  Run.ok = 2;
170  }
171  if (Run.ok == 2 && (Run.date < 130101 || Run.date > 140000 && Run.date < 150000)) return;
172  if (Run.ok == 2) { // try to drift length for Run XIII
173  if (Run.dzWO > 2e-3 || Run.dzEO > 2e-3 ||
174  Run.zWO + Run.zEO < 409.9 ||
175  Run.zWO + Run.zEO > 410.5) {Run.ok = 3; return;}
176  /*
177  RunNT->Draw("1e6*dvAll/vWest-1:zEO+zWO>>D(10,409.9,410.5)","dzWO<1e-3&&dzEO<1e-3&&ddvAll<1e-3&&ok<2","prof")
178  D->Fit("pol1","e")
179  FCN=4.90371 FROM MINOS STATUS=FAILURE 156 CALLS 389 TOTAL
180  EDM=5.83814e-09 STRATEGY= 1 ERR MATRIX NOT POS-DEF
181  EXT PARAMETER APPROXIMATE STEP FIRST
182  NO. NAME VALUE ERROR SIZE DERIVATIVE
183  1 p0 5.85670e-01 2.04554e-04 2.04554e-04 2.08675e-04
184  2 p1 -1.42815e-03 4.98660e-07 4.98660e-07 3.42967e+03
185  */
186  dv = 1e-6*Run.vWest*(1. + (5.85670e-01 -1.42815e-03*(Run.zWO + Run.zEO)));
187  }
188  TString fOut = Form("tpcDriftVelocity.%8i.%06i.C",date,Time);
189  ofstream out;
190  cout << "Create " << fOut << endl;
191  out.open(fOut.Data());
192  out << "TDataSet *CreateTable() {" << endl;
193  out << " if (!gROOT->GetClass(\"St_tpcDriftVelocity\")) return 0;" << endl;
194  out << " St_tpcDriftVelocity *tableSet = new St_tpcDriftVelocity(\"tpcDriftVelocity\",1);" << endl;
195  out << " tpcDriftVelocity_st row;// Laser Run " << run << endl;
196  out << " memset(&row, 0, tableSet->GetRowSize());"<< endl;
197  if (! Run.ok) {// ok == 0 => use both: west and east
198  Run.dvSet = dvWest;
199  out << " row.laserDriftVelocityEast = " << dvEast << "; // +/- " << ddvEast
200  << " cm/us East: Slope = " << DVAll[1][2] << " +/- " << dDVAll[1][2] << " DV = " << dvEast << " +/- " << ddvEast
201  << endl;
202  out << " row.laserDriftVelocityWest = " << dvWest << "; // +/- " << ddvWest
203  << " cm/us West: Slope = " << DVAll[1][1] << " +/- " << dDVAll[1][1] << " DV = " << dvWest << " +/- " << ddvWest<< endl;
204  out << " tableSet->AddAt(&row); " << endl;
205  out << " return (TDataSet *)tableSet; // 1e3*Delta: All = " << dv << " +/- " << ddv << endl;
206  } else { // averaged drif tvelocity
207  Run.dvSet = dv;
208  out << " row.laserDriftVelocityEast = " << dv << "; // +/- " << ddv
209  << " cm/us All: East = " << DVAll[1][2] << " +/- " << dDVAll[1][2];
210  if (Run.ok == 2) out << " From Membrane";
211  out << endl;
212  out << " row.laserDriftVelocityWest = " << dv << "; // +/- " << ddv
213  << " cm/us All: West = " << DVAll[1][1] << " +/- " << dDVAll[1][1] << endl;
214  out << " tableSet->AddAt(&row);// 1e3*Delta: All = " << dv << " +/- " << ddv << endl;
215  out << " return (TDataSet *)tableSet;//"
216  << " West = " << dvWest << " +/- " << ddvWest
217  << " East = " << dvEast << " +/- " << ddvEast << endl;
218  }
219  out << "};" << endl;
220 }
221 //________________________________________________________________________________
222 void Fit() {
223  static TDatime t0(2007,1,1,0,0,0);
224  static UInt_t u0 = t0.Convert();
225  date = 20000000 + (Int_t) Run.date;
226  Time = (Int_t) Run.time;
227  TDatime t(date, Time);
228  Run.utime = t.Convert();
229  Run.day = 1. + (Run.utime - u0)/(24.*60.*60.);
230  run = (Int_t) (1000000*((Int_t) (Run.date/1000000)) + Run.run);
231  memset(&DVAll[0][0], 0, 6*sizeof(Double_t));
232  memset(&dDVAll[0][0], 0, 6*sizeof(Double_t));
233  if (! slope) return;
234  Run.events = slope->GetEntries();
235  cout << "Run " << run << " has " << Run.events << " entries" << endl;
236  TH2D *hist[2] = {dv, slope};
237  Float_t *par = &Run.dvAll;
238  for (Int_t l = 0; l < 2; l++) {
239 #ifdef ADJUSTABLE_BINNING
240  hist[l]->BufferEmpty();
241 #endif
242  // hist[l]->Write();
243  hist[l]->FitSlicesY(0,1,0,10,"QNRI");
244  TString fitN(hist[l]->GetName());
245  fitN += "_1";
246  TH1D *fit = (TH1D *) gDirectory->Get(fitN);
247  if (! fit) continue;
248  // Disable sector 16 for Run 20
249  TAxis *xax = fit->GetXaxis();
250  Int_t sector16 = xax->FindBin(16.);
251  fit->SetBinError(sector16, 1.0);
252  Int_t sector8 = xax->FindBin(8.);
253  fit->SetBinError(sector8, 1.0);
254 
255  fit->SetMarkerStyle(20);
256  Double_t xmin = fit->GetXaxis()->GetXmin();
257  Double_t xmax = fit->GetXaxis()->GetXmax();
258  TF1 *pol0 = (TF1*) gROOT->GetFunction("pol0");
259  if (pol0) {
260  for (Int_t k = 0; k < 3; k++) {
261  pol0->SetLineColor(k+1);
262  Double_t x1 = xmin;
263  Double_t x2 = xmax;
264  TString opt("er");
265  if (k == 1) {
266  x2 = 12.5;
267  } else if (k == 2) {
268  opt += "+";
269  // x1 = 12.5;
270  x1 = 17.0;
271  }
272  DVAll[l][k] = -999;
273  dDVAll[l][k] = 999;
274  Int_t status = fit->Fit(pol0,opt,"",x1,x2);
275  if (! status) {
276  DVAll[l][k] = pol0->GetParameter(0);
277  dDVAll[l][k] = pol0->GetParError(0);
278  }
279  par[2*(k+3*l)] = DVAll[l][k];
280  par[2*(k+3*l)+1] = dDVAll[l][k];
281  }
282  }
283  fit->Write();
284  }
285  // Memberane
286  for (Int_t io = 0; io < 2; io++) {
287  zMembrane[io]->FitSlicesY(0,1,0,10,"QNRI");
288  TString fitN(zMembrane[io]->GetName());
289  fitN += "_1";
290  TH1D *fit = (TH1D *) gDirectory->Get(fitN);
291  if (! fit) continue;
292  fit->SetMarkerStyle(20);
293  Float_t *par = &Run.zWI;
294  for (Int_t we = 0; we < 2; we++) {
295  TF1 *pol0 = (TF1*) gROOT->GetFunction("pol0");
296  if (we == 0) fit->Fit(pol0,"er","",0.5,12.5);
297  else fit->Fit(pol0,"er+","",12.5,24.5);
298  par[4*io+2*we ] = pol0->GetParameter(0);
299  par[4*io+2*we+1] = pol0->GetParError(0);
300  }
301  fit->Write();
302  }
303  // X & Y slope from Memberane
304  for (Int_t xy = 0; xy < 2; xy++) {
305  for (Int_t io = 0; io < 2; io++) {
306  Float_t *par = &Run.XWI;
307  for (Int_t we = 0; we < 2; we++) {
308  dMembraneY[io][we][xy]->FitSlicesY(0,1,0,10,"QNRI");
309  TString fitN(dMembraneY[io][we][xy]->GetName());
310  fitN += "_1";
311  TH1D *fit = (TH1D *) gDirectory->Get(fitN);
312  if (! fit) continue;
313  fit->SetMarkerStyle(20);
314  TF1 *pol1 = (TF1*) gROOT->GetFunction("pol1");
315  Int_t status = fit->Fit(pol1);
316  if (! status) {
317  par[8*xy+4*io+2*we ] = pol1->GetParameter(1)/210;
318  par[8*xy+4*io+2*we+1] = pol1->GetParError(1)/210;
319  } else {
320  par[8*xy+4*io+2*we ] = -999;
321  par[8*xy+4*io+2*we+1] = 999;
322  }
323  fit->Write();
324  }
325  }
326  }
327  MakeTable();
328  runNT->Fill(&Run.run);
329 }
330 //________________________________________________________________________________
331 void LanaTrees(const Char_t *files="./st_laser_*.laser.root", const Char_t *Out = "LaserPlots.root") {
332  TDirIter Dir(files);
333  const Char_t *TreeName = "laser";
334  TChain *tree = new TChain(TreeName);
335  Char_t *file = 0;
336  while ((file = (Char_t *) Dir.NextFile())) {
337  tree->Add(file);
338  }
339  LaserEvent *event = new LaserEvent();
340  tree->SetBranchAddress("event", &event);
341 
342  TFile *fOut = new TFile(Out,"recreate");
343  runNT = new TNtuple("RunNT","Run date time",vRun);
344  static const Double_t EastWRTWestDiff = 0;//3.55700e-01; // +/- 1.77572e-01 3.38530e-01; // +/- 1.39566e-01 permill
345 
346  Double_t OnlFreq = 0;
347  Int_t oldRun = -1;
348  Double_t oldDate = -1;
349  TDatime t0(2007,1,1,0,0,0);
350  UInt_t ut0 = t0.Convert();
351  Long64_t nentries = tree->GetEntriesFast();
352  Long64_t nbytes = 0, nb = 0;
353  for (Long64_t jentry=0; jentry<nentries;jentry++) {
354  Long64_t ientry = tree->LoadTree(jentry);
355  if (ientry < 0) break;
356  nb = tree->GetEntry(jentry); nbytes += nb;
357 #if 0
358  Int_t iok = -1;
359  for (Int_t i = 0; i < NZFR; i++) {if (ZFieldRuns[i] == event->GetHeader()->fRun) {iok = i; break;}}
360  if (iok >= 0) continue;
361 #endif
362  if (event->GetHeader()->fRun != oldRun) {
363  OnlFreq = event->GetHeader()->fOnlClock;
364  fOut->cd();
365  TDatime t(event->GetHeader()->fDate,event->GetHeader()->fTime);
366  UInt_t ut = t.Convert();
367  Double_t Date = 1. + (ut - ut0)/(24.*60.*60.);
368 #ifdef INTEGRATE_OVER_HOURS
369  if (Date - oldDate > 0.125) { // 3/24 = 1/8
370 #endif
371  if (oldRun != -1) Fit();
372  cout << "New run " << event->GetHeader()->fRun << " Date Old/New " << oldDate << "/" << Date << endl;
373  Run.run = event->GetHeader()->fRun%1000000;
374  Run.date = event->GetHeader()->fDate%1000000;
375  Run.time = event->GetHeader()->fTime;
376  Run.vWest = event->GetHeader()->fDriVelWest;
377  Run.vEast = event->GetHeader()->fDriVelEast;
378  Run.events = 0;
379 #ifdef ADJUSTABLE_BINNING
380  dv = new TH2D(Form("DV%i",event->GetHeader()->fRun%1000000),Form("Drift Velocity for run %i",event->GetHeader()->fRun%1000000),12,1,25,400,1,-1);
381  dv->SetBuffer(100000);
382 #else
383  dv = new TH2D(Form("DV%i",event->GetHeader()->fRun%1000000),Form("Drift Velocity for run %i",event->GetHeader()->fRun%1000000),12,1,25,2000,5.3,5.9);
384 #endif
385  dv->SetXTitle("Sector");
386  dv->SetYTitle("Drift Velocity ");
387 #ifdef ADJUSTABLE_BINNING
388  slope = new TH2D(Form("SL%i",event->GetHeader()->fRun%1000000),Form("Slope for run %i",event->GetHeader()->fRun%1000000),12,1,25,400,1,-1);
389  slope->SetBuffer(100000);
390 #else
391  slope = new TH2D(Form("SL%i",event->GetHeader()->fRun%1000000),Form("Slope for run %i",event->GetHeader()->fRun%1000000),12,1,25,2000,-10.,10.);
392 #endif
393  slope->SetXTitle("Sector");
394  slope->SetYTitle("Difference wrt reference Drift Velocity in pemill");
395  memAdc = new TH2D(Form("memAdc%i", event->GetHeader()->fRun%1000000),
396  Form("Log(Adc) @ Membrane for West and East, Inner and Outer sectors for run %i",event->GetHeader()->fRun%1000000),
397  4, 0.5, 4.5, 100, 0, 10);
398  for (Int_t io = 0; io < 2; io++) {// 0 => Inner, 1 => Outer
399  TString name("Z");
400  TString title("Z[cm] of Membrane");
401  if (io == 0) {name += "I"; title += " Inner";}
402  else {name += "O"; title += " Outer";}
403  name += Form("%i",event->GetHeader()->GetRun()%1000000);
404  title += Form(" for run %i",event->GetHeader()->GetRun()%1000000);
405  // dMembraneY[io] = new TH2D(name,title,24,0.5,24.5,2000,-10.,10.);
406  zMembrane[io] = new TH2D(name,title,24,0.5,24.5,200,200,210);
407  for (Int_t we = 0; we < 2; we++) {
408  for (Int_t xy = 0; xy < 2; xy++) {
409  name = "R"; name += Form("%i",event->GetHeader()->GetRun()%1000000);
410  title = "Drift length for Membrane clusters versus";
411  if (xy == 0) {name += "X"; title += " X. ";}
412  else {name += "Y"; title += " Y. ";}
413  if (io == 0) {name += "I"; title += " Inner";}
414  else {name += "O"; title += " Outer";}
415  if (we == 0) {name += "W"; title += " West";}
416  else {name += "E"; title += " East";}
417  dMembraneY[io][we][xy] = new TH2D(name,title,100,-200,200,1000,-10,10);
418  }
419  }
420  }
421  oldRun = (Int_t) event->GetHeader()->GetRun();
422  oldDate = Date;
423 #ifdef INTEGRATE_OVER_HOURS
424  }
425 #endif
426  }
427 #if 1
428  for (Int_t i = 0; i < event->GetNhit(); i++) {
429  Hit *hit = (Hit *) event->Hits()->UncheckedAt(i);
430  Int_t io = 0; if (hit->row > 13) io = 1;
431  zMembrane[io]->Fill(hit->sector,hit->xyzL.z());
432  Int_t we = 0; if (hit->sector > 12) we = 1;
433  dMembraneY[io][we][0]->Fill(hit->xyz.x(),hit->xyzTpcL.z());
434  dMembraneY[io][we][1]->Fill(hit->xyz.y(),hit->xyzTpcL.z());
435  if (TMath::Abs(TMath::Abs(hit->xyzTpcL.z())-205) < 25 && hit->adc > 0) memAdc->Fill(2.*we + io + 1., TMath::Log(hit->adc));
436  }
437 #endif
438  Double_t dt = event->GetHeader()->fDate%1000000 + ((Double_t) event->GetHeader()->fTime)*1e-6;
439  Double_t DT = Run.date + Run.time*1e-6;
440  if (dt < DT) {
441  Run.date = event->GetHeader()->fDate%1000000;
442  Run.time = event->GetHeader()->fTime;
443  }
444  for (Int_t k = 0; k < 12; k++) {
445  FitDV *fit = (FitDV *) event->Fit()->UncheckedAt(k);
446  if (! fit) continue;
447 #ifndef __REFIT__
448  if (fit->ndf > 2 && fit->Prob > 1.e-2) {
449  Double_t Slope = 1e3*fit->slope;
450  if (fit->Sector > 12) Slope -= EastWRTWestDiff;
451  slope->Fill(fit->Sector, Slope);
452  Double_t Vm = OnlFreq/event->GetHeader()->fClock;//*event->GetHeader()->fDriVel;
453  if (fit->Sector <= 12) Vm *= event->GetHeader()->fDriVelWest;
454  else Vm *= event->GetHeader()->fDriVelEast;
455  Double_t V = 1e-6*Vm/(1.+1e-3*Slope-Vm/TMath::Ccgs());
456  dv->Fill(fit->Sector,V);
457  }
458 #else
459  Int_t N = fit->N;
460  if (N > 2) {
461  Double_t Flag[42];
462  memset (Flag, 0, 42*sizeof(Double_t));
463  Double_t YA = 0;
464  for (Int_t i = 0; i < N; i++) {
465  YA += fit->Y[i];
466  }
467  YA /= N;
468  for (Int_t i = 0; i < N; i++) {if (TMath::Abs(fit->Y[i] - YA) > 10) Flag[i] = 1;}
469  for (;;) {
470  TRVector Y;
471  TRMatrix A(0,2);
472  for (Int_t i = 0; i < N; i++) {
473  if (! fit->Flag[i]) {
474  Double_t dev = fit->Y[i]/sigma;
475  Y.AddRow(&dev);
476  Double_t a[2] = {1./sigma, fit->X[i]/sigma};
477  A.AddRow(a);
478  }
479  }
480  Int_t ndf = A.GetNrows() - 2;
481  if (ndf < 1) break;
482  TRSymMatrix S(A,TRArray::kATxA); if (_debug) cout << "S: " << S << endl;
483  TRVector B(A,TRArray::kATxB,Y); if (_debug) cout << "B: " << B << endl;
484  TRSymMatrix SInv(S,TRArray::kInverted); if (_debug) cout << "SInv: " << SInv << endl;
485  TRVector X(SInv,TRArray::kSxA,B); if (_debug) cout << "X: " << X << endl;
486  TRVector R(Y);
487  R -= TRVector(A,TRArray::kAxB,X); if (_debug) cout << "R: " << R << endl;
488  Double_t chisq = R*R;
489  Double_t prob = TMath::Prob(chisq,ndf);
490  if (prob > 0.01) {
491  Double_t offset = X[0];
492  Double_t slope = X[1];
493  Double_t chisq = chisq;
494  Double_t dslope = SInv[2];
495  Double_t doffset = SInv[0];
496  if (ndf > 4 && prob > 1.e-2) {
497  Double_t Slope = 1e3*slope;
498  if (fit->Sector > 12) Slope -= EastWRTWestDiff;
499  slope->Fill(fit->Sector, Slope);
500  Double_t Vm = OnlFreq/event->GetHeader()->fClock*event->GetHeader()->fDriVel;
501  Double_t V = 1e-6*Vm/(1.+1e-3*Slope-Vm/TMath::Ccgs());
502  dv->Fill(fit->Sector,V);
503  }
504  }
505  break;
506  } else {
507  Int_t j = -1;
508  Int_t imax = -1;
509  Double_t Rmax = 0;
510  for (Int_t i = 0; i < N; i++) {
511  if (! fit->Flag[i]) {
512  j++;
513  Double_t RR = TMath::Abs(R[j]);
514  if (RR > Rmax) {
515  imax = i;
516  Rmax = RR;
517  }
518  }
519  }
520  if (imax < 0) break;
521  Flag[imax] = 1;
522  }
523  }
524 #endif
525  }
526  }
527  Fit();
528  // runNT->Write();
529  fOut->Write();
530 }
531 //________________________________________________________________________________
532 /*
533  A.Lebedev 02/11/08
534  Delays for laser rafts.
535  T-zero is a moment, when a laser arrived to TPC wheel surface.
536  All other numbers corresponded to time for laser light to propagate to particular raft (TPC sector).
537  Estimated error in table 0.1ns
538  West: sector 2 4 6 8 12
539  Time(ns) 10.33 3.34 6.14 13.11 17.31
540  East: sector 14 16 18 20 22 24
541  Time(ns) 19.88 12.95 5.97 3.18 10.17 17.14
542 
543 */
544 //________________________________________________________________________________
545 Double_t OffSets(Double_t *x/*, Double_t *par = 0 */) {
546  // TF1 *f = new TF1("Off",OffSets,0,24,0)
547  // 2 4 6 8 10 12
548  static Double_t Toffsets[12] = {10.33, 3.34, 6.14, 13.11, -1., 17.31, // ns
549  19.88, 12.95, 5.97, 3.18, 10.17, 17.14};
550  static Double_t DV = 55.42840e-4; // cm/ns
551  Int_t sector2 = (Int_t) ((x[0]-1)/2.);
552  Double_t off = -1;
553  if (sector2 < 0 || sector2 > 12) return off;
554  return DV*Toffsets[sector2];
555 }
556 /*
557  RunNT->Draw("1e6*dvAll/vWest-1:zWO-zEO>>zO(10,7,7.5)","ok<2&&zEO<0","prof")
558 
559 
560 c1 = new TCanvas()
561 //TH1* frame = c1->DrawFrame(575e6,5.4,581.e6,5.6)
562 TH1* frame = c1->DrawFrame(632e6,5.4,642.e6,5.6)
563 frame->SetTitle("Drift velocitry")
564 frame->GetXaxis()->SetTimeDisplay(1)
565 c1->Update()
566 TLegend *l = new TLegend(0.2,0.2,0.4,0.4)
567 RunNT->Draw("dvWest:utime-788936400>>west","ok==0","same")
568 l->AddEntry(west,"West")
569 l->Draw()
570 RunNT->SetMarkerColor(2)
571 RunNT->Draw("dvEast:utime-788936400>>east","ok==0","same")
572 l->AddEntry(east,"East")
573 RunNT->SetMarkerColor(3)
574 RunNT->Draw("dvAll:utime-788936400>>all","ok==1","same")
575 l->AddEntry(all,"All")
576 RunNT->SetMarkerColor(4)
577 RunNT->Draw("dvSet:utime-788936400>>memb","ok==2","same")
578 l->AddEntry(memb,"Membrane")
579 
580 
581 TH1F *frame = c1->DrawFrame(632e6,5.46,644e6,5.57)
582 .x DrawTime.C(frame)
583 RunNT->Draw("dvAll:utime-788936400>>h2","ok<2","same")
584 
585 RunXIX
586 TH1* frame = c1->DrawFrame(758e6,5.4,767.e6,5.6)
587 frame->SetTitle("Drift velocitry")
588 frame->GetXaxis()->SetTimeDisplay(1)
589 c1->Update()
590 
591  TH1* frame = c1->DrawFrame(786e6,5.4,794.e6,5.8)
592 
593 Run XX
594  root.exe 2020F.root 2020P.root
595 TCanvas *c1 = new TCanvas("c1","c1");
596 TH1* frame = c1->DrawFrame(786e6,5.4,795e6,5.8);
597 frame->SetTitle("Tpc drift verlocity versus date");
598 .x DrawTime.C(frame)
599 
600 _file0->cd();
601 RunNT->SetMarkerColor(2);
602 RunNT->Draw("dvAll:utime-788936400>>W16","ok==1","same");
603 TLegend *l = new TLegend(0.7,0.3,0.8,0.4)
604 l->AddEntry(W16,"With Sector 16");
605 _file1->cd();
606 RunNT->Draw("dvAll:utime-788936400>>WO16","ok==1","same");
607 l->AddEntry(WO16,"Withot Sector 16");
608 l->Draw();
609  */
Run_t
Definition: LanaTrees.C:110
Generated by   doxygen 1.8.5