StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Laser.C
1 // $Id: Laser.C,v 1.5 2007/12/10 19:54:02 fisyak Exp $
2 // $Log: Laser.C,v $
3 // Revision 1.5 2007/12/10 19:54:02 fisyak
4 // Add Id and Log, correct spelling error in README
5 //
6 #if !defined(__CINT__) || defined(__MAKECINT__)
7 //#include <ostream>
8 #include "Riostream.h"
9 #include <stdio.h>
10 #include "TROOT.h"
11 #include "TSystem.h"
12 #include "TMath.h"
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TH3.h"
16 #include "TStyle.h"
17 #include "TF1.h"
18 #include "TProfile.h"
19 #include "TTree.h"
20 #include "TChain.h"
21 #include "TFile.h"
22 #include "TNtuple.h"
23 #include "TCanvas.h"
24 #include "TString.h"
25 #include "TList.h"
26 #include "TLegend.h"
27 #include "TFileSet.h"
28 #include "TDataSetIter.h"
29 #include "TDirIter.h"
30 #include "TTreeIter.h"
31 #endif
32 static Int_t run = 0;
33 static Int_t date = 0;
34 static Int_t Time = 0;
35 // Delta/dv All,West,East
36 static Double_t DVAll[2][3];
37 static Double_t dDVAll[2][3];
38 //________________________________________________________________________________
39 void MakeTable() {
40  TString fOut = Form("tpcDriftVelocity.%8i.%06i.C",date,Time);
41  ofstream out;
42  cout << "Create " << fOut << endl;
43  out.open(fOut.Data());
44  out << "TDataSet *CreateTable() {" << endl;
45  out << " if (!gROOT->GetClass(\"St_tpcDriftVelocity\")) return 0;" << endl;
46  out << " St_tpcDriftVelocity *tableSet = new St_tpcDriftVelocity(\"tpcDriftVelocity\",1);" << endl;
47  out << " tpcDriftVelocity_st row;// Laser Run " << run << endl;
48  out << " row.laserDriftVelocityEast = " << DVAll[1][0] << "; // +/- " << dDVAll[1][0] << " cm/us All: East = " << DVAll[1][2] << " +/- " << dDVAll[1][2] << endl;
49  out << " row.laserDriftVelocityWest = " << DVAll[1][0] << "; // +/- " << dDVAll[1][0] << " cm/us All: West = " << DVAll[1][1] << " +/- " << dDVAll[1][1] << endl;
50  out << " row.cathodeDriftVelocityEast = 0; // cm/us : from cathode emission ;" << endl;
51  out << " row.cathodeDriftVelocityWest = 0; // cm/us : from cathode emission ;" << endl;
52  out << " tableSet->AddAt(&row);// 1e3*Delta: All = " << DVAll[0][0] << " +/- " << dDVAll[0][0] << endl;
53  out << " return (TDataSet *)tableSet;//"
54  << " West = " << DVAll[0][1] << " +/- " << dDVAll[0][1]
55  << " East = " << DVAll[0][2] << " +/- " << dDVAll[0][2] << endl;
56  out << "};" << endl;
57 }
58 //________________________________________________________________________________
59 void Slopes(Int_t n1 = 1, Int_t n2 = 9999) {
60  Int_t NF = 0;
61  TSeqCollection *files = gROOT->GetListOfFiles();
62  Int_t nn = files->GetSize();
63  if (! nn) return;
64  cout << "Oneped " << nn << " Files" << endl;
65  TFile **Files = new TFile *[nn];
66  // TString FileNames[nn];
67  if (! files) return;
68  TIter next(files);
69  TFile *f = 0;
70  Int_t l = 0;
71  TString cut("fFit.Prob>1e-2&&fFit.ndf>4");
72  while ( (f = (TFile *) next()) ) {
73  TTree *laser = (TTree *) f->Get("laser");
74  if (! laser) {
75  cout << "Skip file " << f->GetName() << " with no laser Tree" << endl;
76  delete f; continue;
77  }
78  Long64_t n = laser->Draw("fFit.slope",cut,"goff");
79  if (n < 100) {
80  cout << "Skip file " << f->GetName() << " with " << n << " Entries" << endl;
81  delete f;
82  } else {
83  l++;
84  if (l >= n1 && l <= n2) {
85  cout << l << " Add file " << f->GetName() << " with " << n << " Entries" << endl;
86  // FileNames[NF] = f->GetName();
87  Files[NF] = f;
88  NF++;
89  }
90  // delete f;
91  }
92  }
93  if (!NF) return;
94  Int_t ny = (Int_t) TMath::Sqrt(NF);
95  Int_t nx = NF/ny;
96  if (nx*ny != NF) nx++;
97  cout << "NF = " << NF << " nx " << nx << " ny " << ny << endl;
98  TCanvas *c[2];
99  c[0] = new TCanvas(Form("Slopes %i",n1),Form("Slopes %i",n1));
100  c[1] = new TCanvas(Form("LDV %i",n1),Form("Laser Drift Velocities %i",n1));
101  c[0]->Divide(nx,ny);
102  c[1]->Divide(nx,ny);
103  for (Int_t i = 0; i < NF; i++) {
104  memset(&DVAll[0][0], 0, 6*sizeof(Double_t));
105  memset(&dDVAll[0][0], 0, 6*sizeof(Double_t));
106  f = Files[i];
107  // f = TFile::Open(FileNames[i]);
108  if (! f) continue;
109  cout << "i = " << i << " " << f->GetName() << endl;
110  f->cd();
111  TTree *laser = (TTree *) f->Get("laser");
112  if (! laser) {
113  cout << "Did not find TTree laser in " << f->GetName() << endl;
114  continue;
115  }
116  laser->SetMarkerStyle(20);
117  laser->Draw("fEvtHdr.fRun:fEvtHdr.fDate:fEvtHdr.fTime","","goff",1);
118  run = (Int_t) laser->GetV1()[0];
119  date = (Int_t) laser->GetV2()[0];
120  Time = (Int_t) laser->GetV3()[0];
121  cout << "Run = " << run << " date = " << date << " Time = " << Time << endl;
122  // Int_t color = i+1;
123  // laser->SetMarkerColor(color);
124  TString histN, plot, plotW, Plot, Plot2, profN;
125  TString plotE;
126  TString CutW(cut); CutW += " && fFit.Sector < 13";
127  TString CutE(cut); CutE += " && fFit.Sector > 12";
128  Char_t *cutName[3] = {"All ","West","East"};
129  Char_t *plName[2] = {"SS", "dV"};
130  Char_t *FMT[2] = {"%s dV = (%7.2f +/- %7.2f) x 10^3","%s dV = (%10.5f +/- %10.5f) cm/mksec"};
131  for (Int_t l = 0; l < 2; l++) {
132  c[l]->cd(i+1)->SetLogz(1);
133  histN = Form("%s%i",plName[l],run);
134  if (l == 0) {
135  plotW = "1e3*fFit.slope:fFit.Sector>>";
136  plotE = plotW; plotE += "+";
137  // plotE = "1e3*(fFit.slope+3.63518e-04):fFit.Sector>>+";
138  } else {
139  plotW = "1e-6*fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/(1.+fFit.slope-fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/TMath::Ccgs()):fFit.Sector>>";
140  plotE = plotW; plotE += "+";
141  // plotE = "1e-6*fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/(1.+fFit.slope+3.63518e-04-fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/TMath::Ccgs()):fFit.Sector>>+";
142  }
143  Plot = plotW;
144  profN = histN;
145  profN += "_pfx";
146  Plot += profN;
147  Plot += "(24,1,25)";
148  Long64_t n = laser->Draw(Plot,CutW,"prof");
149  Plot = plotE; Plot += profN;
150  n += laser->Draw(Plot,CutE,"prof");
151  if (n < 10) {
152  cout << f->GetName() << " empty only " << n << " entries" << endl;
153  continue;
154  }
155  TProfile *SSp = (TProfile *) gDirectory->Get(profN);
156  if (! SSp) continue;
157  if (SSp->GetEntries() <= 10) {
158  cout << SSp->GetName() << " empty only " << n << " entries" << endl;
159  continue;
160  }
161  Plot2 = plotW;
162  Plot2 += histN;
163  SSp->SetMarkerColor(3);
164  Double_t ymin = SSp->GetMean(2) - 3*SSp->GetRMS(2);
165  Double_t ymax = SSp->GetMean(2) + 3*SSp->GetRMS(2);
166  Plot2 += Form("(24,1,25,100,%f,%f)",ymin,ymax);
167  laser->Draw(Plot2,CutW,"colz");
168  Plot2 = plotE;
169  Plot2 += histN;
170  laser->Draw(Plot2,CutE,"colz");
171  TH2F *SS = (TH2F *) gDirectory->Get(histN);
172  if (! SS) {
173  cout << "Did not find histogram " << histN << endl;
174  continue;
175  }
176  SS->FitSlicesY();
177  TString fitN(histN);
178  fitN += "_1";
179  TH1D *fit = (TH1D *) gDirectory->Get(fitN);
180  TLegend *leg = new TLegend(0.1,0.1,0.55,0.25,"");
181  if (fit) {
182  fit->SetMarkerStyle(20);
183  TF1 *pol0 = (TF1*) gROOT->GetFunction("pol0");
184  if (pol0) {
185  for (Int_t k = 0; k < 3; k++) {
186  pol0->SetLineColor(k+1);
187  if (k == 0) fit->Fit(pol0,"er","",1,25);
188  if (k == 1) fit->Fit(pol0,"er+","",1,13);
189  if (k == 2) fit->Fit(pol0,"er+","",13,25);
190  DVAll[l][k] = pol0->GetParameter(0);
191  dDVAll[l][k] = pol0->GetParError(0);
192  if (!l || dDVAll[l][k] > 0 && dDVAll[l][k]< 1e-3)
193  leg->AddEntry((TF1 *) fit->GetListOfFunctions()->Last(), Form(FMT[l],cutName[k],DVAll[l][k],dDVAll[l][k]));
194  }
195  }
196  }
197  SS->SetTitle(Form("Run: %i",run));
198  SS->SetXTitle("Sector");
199  SS->SetYTitle("Drift Velocity ");
200  SS->Draw("colz");
201  if (SSp) SSp->Draw("same");
202  if (fit) fit->Draw("same");
203  if (l && DVAll[l][0] > 0) {
204  if (DVAll[l][0] > 5.5 && dDVAll[l][0] > 0 && dDVAll[l][0]< 1e-3) MakeTable();
205  else {
206  cout << "File " << f->GetName() << " fails =============================" << endl;
207  leg->AddEntry(fit, "Rejected");
208  }
209  }
210  leg->Draw();
211  c[l]->Update();
212  }
213  }
214 }
215 #if 0
216 //________________________________________________________________________________
217 void Offsets() {
218  Int_t NF = 0;
219  TSeqCollection *files = gROOT->GetListOfFiles();
220  Int_t nn = files->GetSize();
221  if (! nn) return;
222  TFile **Files = new TFile *[nn];
223  if (! files) return;
224  TIter next(files);
225  TFile *f = 0;
226  while ( (f = (TFile *) next()) ) {
227  TTree *laser = (TTree *) f->Get("laser");
228  if (! laser) continue;
229  cout << "Add " << f->GetName() << endl;
230  Files[NF] = f; NF++;
231  }
232  Int_t ny = (Int_t) TMath::Sqrt(NF);
233  Int_t nx = NF/ny;
234  if (nx*ny != NF) nx++;
235  TCanvas *c1 = new TCanvas("Offsets","Offsets");
236  c1->Divide(nx,ny);
237  for (Int_t i = 0; i < NF; i++) {
238  c1->cd(i+1);
239  f = Files[i];
240  if (! f) continue;
241  TTree *laser = (TTree *) f->Get("laser");
242  if (! laser) continue;
243  laser->SetMarkerStyle(20);
244  // Int_t color = i+1;
245  // laser->SetMarkerColor(color);
246  laser->Draw(Form("abs(fFit.offset)-6.66274:fFit.Sector>>OS%i(24,1,25)",i),"fFit.doffset > 0 && fFit.doffset <5","prof");
247  TProfile *OS = (TProfile *) gDirectory->Get(Form("OS%i",i));
248  TString Name(gSystem->BaseName(f->GetName()));
249  Name.ReplaceAll(".root","");
250  Name.ReplaceAll("laser_","Run ");
251  cout << Name << " i = " << i << endl;
252  OS->SetTitle(Name);
253  OS->SetXTitle("Sector");
254  OS->SetYTitle("dZ(cm)");
255  OS->Draw();
256  c1->Update();
257  // Int_t j;
258  // cin >> j;
259  }
260 }
261 #endif
262 //________________________________________________________________________________
263 void dPhi() {
264  TTree *laser = (TTree *) gDirectory->Get("laser");
265  if (! laser) return;
266  // laser->Draw("fTracks.dPhi:fTracks.mSector>>dPhi(12,1,25,100,-.15,0.15)","fTracks.Flag>1","colz");
267  // laser->Draw("TMath::ATan2(fTracks.Laser.dirU.mX2, fTracks.Laser.dirU.mX2)-TMath::ATan2(fTracks.dirPU.mX2,fTracks.dirPU.mX1):fTracks.mSector>>dPhi(12,1,25,100,-.15,0.15)","fTracks.Flag>1","colz");
268  laser->Draw(
269  "atan(fTracks.Laser.dirU.mX2/fTracks.Laser.dirU.mX1)-atan(fTracks.dirPU.mX2/fTracks.dirPU.mX1):fTracks.mSector>>dPhi(12,1,25,100,-.15,0.15)",
270  "fTracks.Flag>1&&abs(fTracks.dU.mX1)<1&&abs(fTracks.dU.mX2)<1","colz");
271  TH2D *dPhi = (TH2D *) gDirectory->Get("dPhi");
272  if (! dPhi) return;
273  dPhi->FitSlicesY();
274  TH1D *dPhi_1 = (TH1D *) gDirectory->Get("dPhi_1");
275  dPhi_1->SetMarkerStyle(20);
276  dPhi_1->Draw("same");
277  for (Int_t bin = 1; bin <= 12; bin++) {
278  cout << "sector: " << 2*bin << " dPhi: " << 1e3*dPhi_1->GetBinContent(bin) << " +/- "
279  << 1e3*dPhi_1->GetBinError(bin) << "(mrad)" << endl;
280  }
281 }
282 //________________________________________________________________________________
283 void dTheta() {
284  TTree *laser = (TTree *) gDirectory->Get("laser");
285  if (! laser) return;
286  laser->Draw("atan(fTracks.Laser.dirU.mX3)-atan(fTracks.dirPU.mX3):fTracks.mSector>>dTheta(12,1,25,100,-.15,0.15)",
287  "fTracks.Flag>1&&abs(fTracks.dU.mX1)<1&&abs(fTracks.dU.mX2)<1","colz");
288  TH2D *dTheta = (TH2D *) gDirectory->Get("dTheta");
289  if (! dTheta) return;
290  dTheta->FitSlicesY();
291  TH1D *dTheta_1 = (TH1D *) gDirectory->Get("dTheta_1");
292  dTheta_1->SetMarkerStyle(20);
293  dTheta_1->Draw("same");
294  for (Int_t bin = 1; bin <= 12; bin++) {
295  cout << "sector: " << 2*bin << " dTheta: " << 1e3*dTheta_1->GetBinContent(bin) << " +/- "
296  << 1e3*dTheta_1->GetBinError(bin) << "(mrad)" << endl;
297  }
298 }
299 //________________________________________________________________________________
300 void dX() {
301  TTree *laser = (TTree *) gDirectory->Get("laser");
302  if (! laser) return;
303  laser->Draw("fTracks.XyzPU.mX1-fTracks.Laser.XyzU.mX1:7*(3*(fTracks.Laser.Sector-2)+fTracks.Laser.Bundle-1)+fTracks.Laser.Mirror-0.5>>dX(504,0,504,100,-.2,0.2)","fTracks.Flag>1","colz");
304  TH2D *dX = (TH2D *) gDirectory->Get("dX");
305  if (! dX) return;
306  dX->FitSlicesY();
307  TH1D *dX_1 = (TH1D *) gDirectory->Get("dX_1");
308  dX_1->SetMarkerStyle(20);
309  dX_1->Draw("same");
310  for (Int_t bin = 1; bin <= 504; bin++) {
311  Int_t mirror = (bin-1)%7 + 1;
312  Int_t bundle = ((bin-1)/7)%6 + 1;
313  Int_t sector = ((bin-1)/(7*6)+1)*2;
314  cout << "sector: " << sector << " bundle: " << bundle << " mirror: " << mirror
315  << " dX: " << dX_1->GetBinContent(bin) << " +/- "
316  << dX_1->GetBinError(bin) << "(cm)" << endl;
317  }
318 }
319 //________________________________________________________________________________
320 void dY() {
321  TTree *laser = (TTree *) gDirectory->Get("laser");
322  if (! laser) return;
323  laser->Draw("fTracks.XyzPU.mX2-fTracks.Laser.XyzU.mX2:7*(3*(fTracks.Laser.Sector-2)+fTracks.Laser.Bundle-1)+fTracks.Laser.Mirror-0.5>>dY(504,0,504,100,-.2,0.2)","fTracks.Flag>1","colz");
324  TH2D *dY = (TH2D *) gDirectory->Get("dY");
325  if (! dY) return;
326  dY->FitSlicesY();
327  TH1D *dY_1 = (TH1D *) gDirectory->Get("dY_1");
328  dY_1->SetMarkerStyle(20);
329  dY_1->Draw("same");
330  for (Int_t bin = 1; bin <= 504; bin++) {
331  Int_t mirror = (bin-1)%7 + 1;
332  Int_t bundle = ((bin-1)/7)%6 + 1;
333  Int_t sector = ((bin-1)/(7*6)+1)*2;
334  cout << "sector: " << sector << " bundle: " << bundle << " mirror: " << mirror
335  << " dY: " << dY_1->GetBinContent(bin) << " +/- "
336  << dY_1->GetBinError(bin) << "(cm)" << endl;
337  }
338 }
339 //________________________________________________________________________________
340 void FillHists(Int_t n1 = 0,const Char_t *Dir = ".") {
341  TFileSet *fileset = new TFileSet(Dir);
342  TDataSetIter iter(fileset);
343  Int_t NF = 0;
344  Int_t NFiles = 0;
345  const Int_t NperC = 64;
346  TString cut("fFit.Prob>1e-2&&fFit.ndf>4");
347  TDataSet *set = 0;
348  while ((set = iter())) {
349  TString Title(set->GetTitle());
350  if (Title != "file") continue;
351  TString file(set->GetName());
352  // cout << file << end;
353  if (!file.BeginsWith("laser") || !file.EndsWith(".root")) continue;
354  // cout << " +++++++ accepted" << endl;
355  TFile *f = new TFile(file);
356  TTree *laser = (TTree *) f->Get("laser");
357  if (! laser) {
358  cout << "Skip file " << f->GetName() << " with no laser TTree" << endl;
359  }
360  else {
361  Long64_t n = laser->Draw("fFit.slope",cut,"goff");
362  if (n < 100) {
363  cout << "Skip file " << f->GetName() << " with " << n << " Entries" << endl;
364  } else {
365  NFiles++;
366  if (n1 >= 0) {
367  if (NFiles <= NperC*n1) continue;
368  if (NFiles > NperC*(n1+1)) break;
369  }
370  cout << "Add file " << f->GetName() << " with " << n << " Entries" << endl;
371  set->Mark();
372  NF++;
373  }
374  }
375  delete f;
376  }
377  if (NF == 0) return;
378  Int_t ny = (Int_t) TMath::Sqrt(NF);
379  Int_t nx = NF/ny;
380  if (nx*ny != NF) nx++;
381  cout << "NF = " << NF << " nx " << nx << " ny " << ny << endl;
382  TCanvas *c[2] = {0,0};
383  TString OutFile("LaserOut");
384  if (n1 >= 0) OutFile += n1;
385  OutFile += ".root";
386  TFile *fOut = new TFile(OutFile,"recreate");
387  iter.Reset();
388  Int_t i = 0;
389  // Int_t i1 = 0;
390  if (n1 < 0) n1 = 0;
391  while ((set = iter())) {
392  if (! set->IsMarked()) continue;
393  memset(&DVAll[0][0], 0, 6*sizeof(Double_t));
394  memset(&dDVAll[0][0], 0, 6*sizeof(Double_t));
395  TString file(set->GetName());
396  cout << "Open " << file << endl;
397  TFile *f = new TFile(file);
398  TTree *laser = (TTree *) f->Get("laser");
399  if (! laser) {
400  cout << "Cannot find laser TTRee in marked file " << file << endl;
401  continue;
402  }
403  if (! c[0] || ! c[1]) {
404  c[0] = new TCanvas(Form("Slopes %i",n1),Form("Slopes %i",n1));
405  c[1] = new TCanvas(Form("LDV %i",n1),Form("Laser Drift Velocities %i",n1));
406  c[0]->Divide(nx,ny);
407  c[1]->Divide(nx,ny);
408  }
409  laser->SetMarkerStyle(20);
410  laser->Draw("fEvtHdr.fRun:fEvtHdr.fDate:fEvtHdr.fTime","","goff",1);
411  run = (Int_t) laser->GetV1()[0];
412  date = (Int_t) laser->GetV2()[0];
413  Time = (Int_t) laser->GetV3()[0];
414  cout << "Run = " << run << " date = " << date << " Time = " << Time << endl;
415  // Int_t color = i+1;
416  // laser->SetMarkerColor(color);
417  TString histN, plot, plotW, Plot, Plot2, profN;
418  TString plotE;
419  TString CutW(cut); CutW += " && fFit.Sector < 13";
420  TString CutE(cut); CutE += " && fFit.Sector > 12";
421  Char_t *cutName[3] = {"All ","West","East"};
422  Char_t *plName[2] = {"SS", "dV"};
423  Char_t *FMT[2] = {"%s dV = (%7.2f +/- %7.2f) x 10^3","%s dV = (%10.5f +/- %10.5f) cm/mksec"};
424  for (Int_t l = 0; l < 2; l++) {
425  c[l]->cd(i+1)->SetLogz(1);
426  histN = Form("%s%i",plName[l],run);
427  if (l == 0) {
428  plotW = "1e3*fFit.slope:fFit.Sector>>";
429  plotE = plotW; plotE += "+";
430  // plotE = "1e3*(fFit.slope+3.63518e-04):fFit.Sector>>+";
431  } else {
432  plotW = "1e-6*fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/(1.+fFit.slope-fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/TMath::Ccgs()):fFit.Sector>>";
433  plotE = plotW; plotE += "+";
434  // plotE = "1e-6*fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/(1.+fFit.slope+3.63518e-04-fEvtHdr.fOnlClock/fEvtHdr.fClock*fEvtHdr.fDriVel/TMath::Ccgs()):fFit.Sector>>+";
435  }
436  Plot = plotW;
437  profN = histN;
438  profN += "_pfx";
439  Plot += profN;
440  Plot += "(24,1,25)";
441  fOut->cd();
442  Long64_t n = laser->Draw(Plot,CutW,"prof");
443  Plot = plotE; Plot += profN;
444  n += laser->Draw(Plot,CutE,"prof");
445  if (n < 10) {
446  cout << f->GetName() << " empty only " << n << " entries" << endl;
447  continue;
448  }
449  TProfile *SSp = (TProfile *) gDirectory->Get(profN);
450  if (! SSp) continue;
451  if (SSp->GetEntries() <= 10) {
452  cout << SSp->GetName() << " empty only " << n << " entries" << endl;
453  continue;
454  }
455  Plot2 = plotW;
456  Plot2 += histN;
457  SSp->SetMarkerColor(3);
458  Double_t ymin = SSp->GetMean(2) - 3*SSp->GetRMS(2);
459  Double_t ymax = SSp->GetMean(2) + 3*SSp->GetRMS(2);
460  Plot2 += Form("(24,1,25,100,%f,%f)",ymin,ymax);
461  laser->Draw(Plot2,CutW,"colz");
462  Plot2 = plotE;
463  Plot2 += histN;
464  laser->Draw(Plot2,CutE,"colz");
465  TH2F *SS = (TH2F *) gDirectory->Get(histN);
466  if (! SS) {
467  cout << "Did not find histogram " << histN << endl;
468  continue;
469  }
470  SS->FitSlicesY();
471  TString fitN(histN);
472  fitN += "_1";
473  TH1D *fit = (TH1D *) gDirectory->Get(fitN);
474  TLegend *leg = new TLegend(0.1,0.1,0.55,0.25,"");
475  if (fit) {
476  fit->SetMarkerStyle(20);
477  TF1 *pol0 = (TF1*) gROOT->GetFunction("pol0");
478  if (pol0) {
479  for (Int_t k = 0; k < 3; k++) {
480  pol0->SetLineColor(k+1);
481  if (k == 0) fit->Fit(pol0,"er","",1,25);
482  if (k == 1) fit->Fit(pol0,"er+","",1,13);
483  if (k == 2) fit->Fit(pol0,"er+","",13,25);
484  DVAll[l][k] = pol0->GetParameter(0);
485  dDVAll[l][k] = pol0->GetParError(0);
486  if (!l || dDVAll[l][k] > 0 && dDVAll[l][k]< 1e-3)
487  leg->AddEntry((TF1 *) fit->GetListOfFunctions()->Last(), Form(FMT[l],cutName[k],DVAll[l][k],dDVAll[l][k]));
488  }
489  }
490  }
491  SS->SetTitle(Form("Run: %i",run));
492  SS->SetXTitle("Sector");
493  SS->SetYTitle("Drift Velocity ");
494  SS->Draw("colz");
495  if (SSp) SSp->Draw("same");
496  if (fit) fit->Draw("same");
497  if (l && DVAll[l][0] > 0) {
498  if (DVAll[l][0] > 5.5 && dDVAll[l][0] > 0 && dDVAll[l][0]< 1e-3) MakeTable();
499  else {
500  cout << "File " << f->GetName() << " fails =============================" << endl;
501  leg->AddEntry(fit, "Rejected");
502  }
503  }
504  leg->Draw();
505  c[l]->Update();
506  }
507  i++;
508 // if (n1 > 1 && i >= nx*ny) {
509 // i1++; i = 0;
510 // }
511  }
512  fOut->Write();
513 }
514 //________________________________________________________________________________
515 void Laser(Int_t n1 = -1){
516  // Slopes();
517  FillHists(n1);
518 }