StRoot  1
StVertexSeedMaker.cxx
1 // //
3 // StVertexSeedMaker class //
4 // Author: G. Van Buren, BNL //
5 // Description: Calculates mean primary vertex positions from //
6 // suitable events to use as seeds in finding better //
7 // primary vertex positions (helpful for low //
8 // multiplicity events like pp collisions). //
9 // See the header file for more information. //
10 // //
12
13 #include "Stiostream.h"
14 #include "StVertexSeedMaker.h"
15 #include "StDbLib/StDbManager.hh"
16 #include "StDbLib/StDbConfigNode.hh"
17 #include "StDbLib/StDataBaseI.hh"
18 #include "StDetectorDbMaker/StDetectorDbTriggerID.h"
19 #include "StMessMgr.h"
20 #include "StVertexId.h"
21 #include "tables/St_vertexSeed_Table.h"
22 #include "tables/St_vertexSeedTriggers_Table.h"
23 #include "tables/St_beamInfo_Table.h"
24 #include "TSystem.h"
25 #include "TFile.h"
26 #include "TVirtualFitter.h"
27 #include "TNtuple.h"
28 #include "TNtupleD.h"
29 #include "TTimeStamp.h"
30 #include "TEventList.h"
31 #include "TArrayF.h"
32 #include "StTree.h"
33 #include "TMath.h"
34 //_____________________________________________________________________________
35 // C variables and functions for fit/minimization
36 //_____________________________________________________________________________
37 static TArrayF xVert, yVert, zVert, multA, exVert, eyVert;
38 int nverts,nsize;
39 double beamWidth;
40 double funcX(float z,Double_t *par) {
41  double x = par[0] + par[1]*z;
42  return x;
43 }
44 double funcY(float z,Double_t *par) {
45  double y = par[2] + par[3]*z;
46  return y;
47 }
48 void fnch(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
49  //calculate chisquare
50  double chisq = 0;
51  for (int i=0;i<nverts; i++) {
52  if (exVert[i]==0) {
53
54  double delta_sq, error_sq;
55  // Error1 set such that 5 tracks => ~7mm, proportional to 1/::sqrt(mult)
56  // This was determined by trying to get the chisq/dof distribution
57  // to peak near 1.0
58  // => Error1 = (0.7cm) / ::sqrt(mult/5)
59  // => Error1^2 = (2.45) cm^2 / mult
60  // Beam spot size (sigma) Error2 ~= 0.04 cm (400 microns)
61  // => Error2^2 = 0.0016 cm^2
62  // The total error should be ::sqrt(Error1^2 + Error2^2)
63  error_sq = 0.0016 + (2.45 / multA[i]);
64  delta_sq = TMath::Power(xVert[i]-funcX(zVert[i],par),2) +
65  TMath::Power(yVert[i]-funcY(zVert[i],par),2);
66  chisq += (delta_sq/error_sq);
67
68  } else {
69
70  //chisq += TMath::Power((xVert[i]-funcX(zVert[i],par))/exVert[i],2) +
71  // TMath::Power((yVert[i]-funcY(zVert[i],par))/eyVert[i],2);
72
73  double errx2 = exVert[i]*exVert[i] + beamWidth*beamWidth;
74  double erry2 = eyVert[i]*eyVert[i] + beamWidth*beamWidth;
75  chisq += TMath::Power(xVert[i]-funcX(zVert[i],par),2)/errx2 +
76  TMath::Power(yVert[i]-funcY(zVert[i],par),2)/erry2;
77
78  }
79  }
80  f = chisq;
81 }
82 void setArraySize(int n) {
83  xVert.Set(n);
84  yVert.Set(n);
85  zVert.Set(n);
86  multA.Set(n);
87  exVert.Set(n);
88  eyVert.Set(n);
89  nsize = n;
90 }
91 void addVert(float x, float y, float z, float m, float ex, float ey) {
92  if (nverts >= nsize) setArraySize(nsize*2);
93  xVert[nverts] = x;
94  yVert[nverts] = y;
95  zVert[nverts] = z;
96  multA[nverts] = m;
97  exVert[nverts] = ex;
98  eyVert[nverts] = ey;
99  nverts++;
100 }
101 //_____________________________________________________________________________
102
103
104 ClassImp(StVertexSeedMaker)
105 //_____________________________________________________________________________
106 StVertexSeedMaker::StVertexSeedMaker(const char *name,
107  const char* defaultDir):StMaker(name){
108  // items which should be initialized only once
109  xdist = new TH1F("xdist","xdist",1000,HIST_MIN,HIST_MAX);
110  ydist = new TH1F("ydist","ydist",1000,HIST_MIN,HIST_MAX);
111  xerr = new TH1F("xerr","x measured - x guess",1000,HIST_MIN,HIST_MAX);
112  yerr = new TH1F("yerr","y measured - y guess",1000,HIST_MIN,HIST_MAX);
113  defDir = defaultDir;
114  mTempOut = 0;
115  resNtuple = 0;
116  parsNtuple = 0;
117  nsize = 0;
118  foffset = 0;
119  noclobber = kTRUE;
120  setArraySize(512);
121  UseEventDateTime(); // By default, use the data & time from the first event
122  useAllTriggers = kFALSE;
123  dbTriggersTable = 0;
124  minEntries = 100; //require 100 valid verts for a seed determination
125  maxX0Err = 0.05; //vertex x should be good to 500 microns
126  maxY0Err = 0.05; //vertex y should be good to 500 microns
127  mHistOut=kTRUE;
128  zVertexMax = 100.0;
129  zVertexMin = -100.0;
130  r2VertexMax = 15.0;
131  HIST_MIN = -3.0;
132  HIST_MAX = 3.0;
133  Reset();
134 }
135 //_____________________________________________________________________________
136 void StVertexSeedMaker::Reset() {
137  // items which should be initialized for every new calibration measurement
138  Clear("");
139  xdist->Reset();
140  ydist->Reset();
141  xerr->Reset();
142  yerr->Reset();
143
144  if (mTempOut) delete mTempOut;
145  else delete resNtuple;
146  mTempOut = new TFile(Form("%s/vertexseedhist.%d.root",
147  gSystem->TempDirectory(),
148  gSystem->GetPid()),"RECREATE");
149  resNtuple = new TNtuple("resNtuple","resNtuple","event:x:y:z:mult:trig:run:fill:zdc:rank:itpc:otpc:detmap:ex:ey:index:bmatch:ematch:tmatch:cmatch:hmatch:pmatch:pct:vpdz:tDay:tFill");
150  LOG_INFO << "Opening new temp file at " << mTempOut->GetName() << endm;
151
152  date = 0;
153  time = 0;
154  nverts = 0;
155  sumzdc = 0;
156
157  a[0] = -888.0;
158  //a[0] = 0.0;
159  a[1] = 0.0;
160  a[2] = 0.0;
161  a[3] = 0.0;
162  chi = 0.0;
163 }
164 //_____________________________________________________________________________
165 StVertexSeedMaker::~StVertexSeedMaker(){
166 }
167 //_____________________________________________________________________________
168 Int_t StVertexSeedMaker::Init(){
173  return StMaker::Init();
174 }
175 //_____________________________________________________________________________
176 void StVertexSeedMaker::Clear(Option_t *option){
177  // items which should be initialized for every event
178  fill = -1;
179  run = -1;
180  zdc = -1;
181  xguess = 0;
182  yguess = 0;
183  zvertex = -999.0;
184  yvertex = -999.0;
185  xvertex = -999.0;
186  exvertex = 0;
187  eyvertex = 0;
188  vpd_zvertex = -999.0;
189  rank = 0;
190  pvn = 0;
191  itpc = 0;
192  otpc = 0;
193  detmap = 0;
194  bmatch = 0;
195  ematch = 0;
196  tmatch = 0;
197  cmatch = 0;
198  hmatch = 0;
199  pmatch = 0;
200  pct = 0;
201  timeEvent = -1;
202  timeFill = -1;
203 }
204 //_____________________________________________________________________________
206  if (date==0) GetADateTime();
207
208  // Currently only finds database values for first event
209  if (a[0] == -888) {
210  LOG_INFO << "Reading db values at the start..." << endm;
211  int status = FillAssumed();
212  if (status != kStOk) return status;
213  status = GetVertexSeedTriggers();
214  if (status != kStOk) return status;
215  }
216
217  if (CheckTriggers()) { // returns false for good events
218  LOG_INFO << "event does not satisfy triggers" << endm;
219  return kStOk;
220  }
221
222  int eventResult = GetEventData();
223  if (eventResult != kStOk) return eventResult;
224
225  if (zvertex<-998) {
226  LOG_INFO << "No primary vertex" << endm;
227  return kStOk;
228  }
229
230  // Calculate guessed vertex x & y from assumed params and measured z
231  xguess = a[0] + (a[1] * zvertex);
232  yguess = a[2] + (a[3] * zvertex);
233  LOG_INFO << "x guess = " << xguess << endm;
234  LOG_INFO << "y guess = " << yguess << endm;
235
236  // Check to see if vertex is good for use in the fit
237  float r2vertex = xvertex*xvertex + yvertex*yvertex;
238  if ((zvertex > zVertexMin) && (zvertex < zVertexMax) &&
239  (mult >= 5) && (r2vertex < r2VertexMax)){
240  xdist->Fill(xvertex);
241  xerr ->Fill(xvertex-xguess);
242  ydist->Fill(yvertex);
243  yerr ->Fill(yvertex-yguess);
244
245  float XX[26];
246  XX[0] = (float) GetEventNumber();
247  XX[1] = xvertex;
248  XX[2] = yvertex;
249  XX[3] = zvertex;
250  XX[4] = mult;
251  XX[5] = trig;
252  XX[6] = (float) (run % 1000000);
253  XX[7] = (float) fill;
254  XX[8] = zdc;
255  XX[9] = rank;
256  XX[10] = (float) itpc;
257  XX[11] = (float) otpc;
258  XX[12] = (float) detmap; // likely to be valid only up to ~23 bits
259  XX[13] = exvertex;
260  XX[14] = eyvertex;
261  XX[15] = (float) pvn;
262  XX[16] = (float) bmatch;
263  XX[17] = (float) ematch;
264  XX[18] = (float) tmatch;
265  XX[19] = (float) cmatch;
266  XX[20] = (float) hmatch;
267  XX[21] = (float) pmatch;
268  XX[22] = (float) pct;
269  XX[23] = vpd_zvertex;
270  XX[24] = (float) (timeEvent % 86400);
271  XX[25] = (float) (timeFill >=0 ? timeEvent - timeFill : timeFill);
272  resNtuple->Fill(XX);
274  sumzdc += zdc;
275  }
276
277  return kStOk;
278 }
279 //_____________________________________________________________________________
281  FindResult();
282  return StMaker::Finish();
283 }
284 //_____________________________________________________________________________
285 bool StVertexSeedMaker::ValidTrigger(unsigned int tid) {
286  // Determine if trigger id is among valid set
287  if (!dbTriggersTable) return kTRUE; // running without DB access
288  vertexSeedTriggers_st* trigsTable = dbTriggersTable->GetTable();
289  int nTrigs = (int) (dbTriggersTable->GetNRows());
290  for (int i = 0; i < nTrigs; i++, trigsTable++) {
291  unsigned int dbTid = trigsTable->trigid;
292  if (useAllTriggers || dbTid == 9999999 || (dbTid > 0 && tid == dbTid)) {
293  trig = (float) tid;
294  return kTRUE;
295  }
296  }
297  return kFALSE;
298 }
299 //_____________________________________________________________________________
300 void StVertexSeedMaker::FindResult(bool checkDb) {
301  bool writeIt = kFALSE;
302  if (nverts >= minEntries){
303  FitData();
304  if (ep[0] > maxX0Err){
305  LOG_ERROR << "x unstable. x0 error = " << ep[0] << " cm." << endm;
306  }
307  if (ep[2] > maxY0Err){
308  LOG_ERROR << "y unstable. y0 error = " << ep[2] << " cm." << endm;
309  }
310  if ((ep[0] <= maxX0Err) && (ep[2] <= maxY0Err)) {
311  if (checkDb) {
312  // Do comparison of this data with data from database to see if
313  // values have changed or improved.
314  LOG_INFO << "Reading db values at the end..." << endm;
315  int status = FillAssumed();
316  if (status == kStOk) {
317  if (ChangedValues() || BetterErrors()) writeIt = kTRUE;
318  else { LOG_INFO << "Values have not changed/improved." << endm; }
319  } else {
320  LOG_WARN << "Could not get db values." << endm;
321  writeIt = kTRUE;
322  }
323  } else {
324  writeIt = kTRUE;
325  }
326  }
327  } else {
328  LOG_ERROR << "Insufficient stats for " <<
329  "mean vertex determination.\n Only " << nverts << " entries." << endm;
330  }
331
332  if (writeIt) WriteTableToFile();
333  else { LOG_WARN << "Not writing table!!!!!" << endm; }
334
335  if (mHistOut) WriteHistFile(writeIt);
336 }
337 //_____________________________________________________________________________
338 void StVertexSeedMaker::PrintInfo() {
339  LOG_INFO << "\n**************************************************************"
340  << "\n* \$Id: StVertexSeedMaker.cxx,v 1.64 2017/08/08 03:58:20 genevb Exp \$"
341  << "\n**************************************************************" << endm;
342
343  if (Debug()) StMaker::PrintInfo();
344 }
345 //_____________________________________________________________________________
346 void StVertexSeedMaker::WriteTableToFile(){
347  TString fileName = NameFile("table","vertexSeed","C");
348  ofstream *out = new ofstream(fileName.Data());
349  St_vertexSeed* vertexSeedTable = VertexSeedTable();
350  vertexSeedTable->SavePrimitive(*out,"");
351  if (parsNtuple) AddResults(parsNtuple);
352  delete out;
353  delete vertexSeedTable;
354  return;
355 }
356 //_____________________________________________________________________________
357 void StVertexSeedMaker::AddResults(TNtupleD* ntup){
358  double datetime = ((double) date) + 1e-6*((double) time);
359  TTimeStamp dt1(date,time,0,true,0);
360  dt1.SetSec(dt1.GetSec() - 4*60*60); // convert GMT -> EDT
361  int tm = dt1.GetTime();
362  int se = tm%100; int hm = (tm-se)/100;
363  int mn = hm%100; int hr = (hm-mn)/100;
364  double days = ((double) dt1.GetDayOfYear()) +
365  ((double) ((hr*60+mn)*60+se))/(24.*60.*60.);
366  ntup->Fill(days,p[0],ep[0],p[2],ep[2],p[1],ep[1],p[3],ep[3],
367  (double) nverts, datetime, (double) fill,
368  sumzdc/((double) nverts),chi,beamWidth);
369 }
370 //_____________________________________________________________________________
371 St_vertexSeed* StVertexSeedMaker::VertexSeedTable(){
372  // up to the user of this function to delete the table
373  St_vertexSeed* table = new St_vertexSeed("vertexSeed",1);
374  vertexSeed_st* row = table->GetTable();
375  row->x0 = p[0];
376  row->dxdz = p[1];
377  row->y0 = p[2];
378  row->dydz = p[3];
379  row->err_x0 = ep[0];
380  row->err_dxdz = ep[1];
381  row->err_y0 = ep[2];
382  row->err_dydz = ep[3];
383  row->chisq_dof = chi;
384  row->weight = 100.0; // Fixed for now!!!
385  row->stats = (float) nverts;
386  table->SetNRows(1);
387  return table;
388 }
389 //_____________________________________________________________________________
390 void StVertexSeedMaker::WriteHistFile(bool writeFit){
391  if (resNtuple->GetEntries() == 0) {
392  LOG_INFO << "Not writing histograms - no entries!!!" << endm;
393  return;
394  }
395  // .ROOT is NOT a typo !!!
396  TString fileName = NameFile("histograms","vertexseedhist","ROOT");
397  if (mTempOut) {
398  mTempOut->Write();
399  mTempOut->Close();
400  if (gSystem->CopyFile(mTempOut->GetName(),fileName.Data()) ||
402  LOG_ERROR << "Could not copy and/or delete temp vertexseedhist file!" << endm;
403  }
404  delete mTempOut;
405  mTempOut = 0;
406  // resNtuple disappears if & when mTempOut->Close() is called
407  } else delete resNtuple;
408  resNtuple = 0;
409
410  TFile out(fileName.Data(),"UPDATE");
411  GetHistList()->Write();
412  if (writeFit) {
413  // no need to garbage collect from newBLpars(): out.Close() takes care of it
415  out.Write();
416  }
417  out.Close();
418 }
419 //_____________________________________________________________________________
420 TString StVertexSeedMaker::NameFile(const char* type, const char* prefix, const char* suffix) {
421  int fdate = date;
422  int ftime = time;
423  if (foffset) { // apply any time offsets
424  TDatime fdatime(date,time);
425  fdatime.Set(fdatime.Convert()+foffset);
426  fdate = fdatime.GetDate();
427  ftime = fdatime.GetTime();
428  }
429  TString fileNameBase = Form("%s.%08d.%06d.%s",prefix,fdate,ftime,suffix);
430
431  if (defDir.Length()>0 && !defDir.EndsWith("/")) defDir.Append("/");
432  TString fileName = defDir;
433  fileName += fileNameBase;
434  LOG_INFO << "Writing new " << type << " to:\n "
435  << fileName << endm;
436  TString dirname = gSystem->DirName(fileName.Data());
437  if (gSystem->OpenDirectory(dirname.Data())==0) {
438  if (gSystem->mkdir(dirname.Data())) {
439  LOG_WARN << "Directory creation failed for:\n " << dirname
440  << "\n Putting " << type << " file in current directory" << endm;
441  fileName = fileNameBase;
442  }
443  }
444  TString searchFile = fileName;
445  if (gSystem->FindFile(".",searchFile)) {
446  if (noclobber) {
447  foffset++;
448  LOG_WARN << "Existing file: trying 1 second later (offset="
449  << foffset << ")..." << endm;
450  fileName = NameFile(type,prefix,suffix);
451  foffset--;
452  } else {
453  LOG_WARN << "Existing file: overwriting!" << endm;
454  }
455  }
456  return fileName;
457 }
458 //_____________________________________________________________________________
459 int StVertexSeedMaker::FillAssumed(){
460  TDataSet* dbDataSet = GetDataBase("Calibrations/rhic/vertexSeed");
461  memset( a,0,4*sizeof(double));
462  memset(ea,0,4*sizeof(double));
463  if (!dbDataSet) {
464  LOG_WARN << "Could not find Calibrations/rhic/vertexSeed in database" << endm;
465  } else {
466  St_vertexSeed* dbTableC =
467  static_cast<St_vertexSeed*>(dbDataSet->FindObject("vertexSeed"));
468  if (!dbTableC) {
469  LOG_WARN << "Could not find vertexSeed in database" << endm;
470  } else {
471  vertexSeed_st* dbTable = dbTableC->GetTable();
472  a[0] = dbTable->x0;
473  a[1] = dbTable->dxdz;
474  a[2] = dbTable->y0;
475  a[3] = dbTable->dydz;
476  ea[0] = dbTable->err_x0;
477  ea[1] = dbTable->err_dxdz;
478  ea[2] = dbTable->err_y0;
479  ea[3] = dbTable->err_dydz;
480  }
481  }
482  LOG_INFO << "Assumed values:"
483  << "\n x0 assumed = " << a[0] << " +/- " << ea[0]
484  << "\n dxdz assumed = " << a[1] << " +/- " << ea[1]
485  << "\n y0 assumed = " << a[2] << " +/- " << ea[2]
486  << "\n dydz assumed = " << a[3] << " +/- " << ea[3]
487  << endm;
488  return kStOk;
489 }
490 //_____________________________________________________________________________
491 int StVertexSeedMaker::GetVertexSeedTriggers(){
492  TDataSet* dbDataSet = GetDataBase("Calibrations/rhic/vertexSeedTriggers");
493  if (!dbDataSet) {
494  LOG_ERROR << "Could not find Calibrations/rhic/vertexSeedTriggers in database" << endm;
495  return kStErr;
496  }
497  dbTriggersTable =
498  static_cast<St_vertexSeedTriggers*>(dbDataSet->FindObject("vertexSeedTriggers"));
499  if (!dbTriggersTable) {
500  LOG_ERROR << "Could not find vertexSeedTriggers in database" << endm;
501  return kStErr;
502  }
503  return kStOk;
504 }
505 //_____________________________________________________________________________
506 bool StVertexSeedMaker::BetterErrors(){
507  bool better = kFALSE;
508  if ((ep[0] < ea[0]) || (ep[1] < ea[1]) ||
509  (ep[2] < ea[2]) || (ep[3] < ea[3])) better = kTRUE;
510  if (better) { LOG_INFO << "Values have improved" << endm; }
511  return better;
512 }
513 //_____________________________________________________________________________
514 bool StVertexSeedMaker::ChangedValues(){
515  bool changed = kFALSE;
516  for (int i = 0; i<4; i++) {
517  double diff = TMath::Abs(p[i] - a[i]);
518  if ((diff >= ep[i]) || (diff >= ea[i])) changed = kTRUE;
519  }
520  if (changed) { LOG_INFO << "Values have changed" << endm; }
521  return changed;
522 }
523 //_____________________________________________________________________________
524 void StVertexSeedMaker::GetADateTime() {
525  date = GetDate();
526  time = GetTime();
527  LOG_INFO << "event date = " << date << endm;
528  LOG_INFO << "event time = " << time << endm;
529  if (!useEventDateTime) GetFillDateTime();
530 }
531 //_____________________________________________________________________________
532 void StVertexSeedMaker::GetFillDateTime() {
533
535  StDbConfigNode* node=mgr->initConfig("RunLog_onl");
537  StDataBaseI* db=mgr->findDb("RunLog_onl");
538  unsigned int ts;
539  char queryStr[128];
540
541  // Find beamInfo entry for this date time
542  sprintf(queryStr,"%08d %06d",date,time);
543  TString tdStr = queryStr;
544  tdStr.Insert(4,'-').Insert(7,'-').Insert(13,':').Insert(16,':');
545  const char* tdstr = tdStr.Data();
546  sprintf(queryStr,
547  " where beginTime<='%s' and deactive=0 order by beginTime desc limit 1",
548  tdstr);
549  ts = db->QueryDb(tab,queryStr);
550
551  if (ts) {
552  // Find earliest entry for this fill
553  run = *(static_cast<int*>(tab->getDataValue("runNumber",0)));
554  LOG_INFO << tdstr << " is from run " << run << endm;
555  float thisFill = *(static_cast<float*>(tab->getDataValue("blueFillNumber",0)));
556  sprintf(queryStr,
557  " where blueFillNumber=%f and deactive=0 order by beginTime asc limit 1",
558  thisFill);
559  db->QueryDb(tab,queryStr);
560
561  // Extract date and time at start of fill
562  char* start = tab->getBeginDateTime();
563  fill = (int) thisFill;
564  time = atoi(&(start[8]));
565  start[8] = 0;
566  date = atoi(start);
567  timeFill = tab->getBeginTime();
568
569  LOG_INFO << "Using fill no. = " << fill << endm;
570  LOG_INFO << "Using fill date = " << date << endm;
571  LOG_INFO << "Using fill time = " << time << endm;
572  } else {
573  LOG_WARN << "Could not find beamInfo in database\n" <<
574  " Using event date/time." << endm;
575  UseEventDateTime();
576  }
577 }
578 //_____________________________________________________________________________
579 void StVertexSeedMaker::FitData() {
580  LOG_INFO << "Now fitting the data..." <<
581  "\n *****************************************************" << endm;
582  TVirtualFitter *minuit = TVirtualFitter::Fitter(0,26);
583  // above '26' must be >25, or ROOT dictates a max of 3 params?!?
584  minuit->SetFCN(fnch);
585
586 // Set starting values and step sizes for parameters
587  static Double_t pstart[4] = {0., 0., 0., 0.};
588  static Double_t pstep[4] = {0.0001, 0.0000001, 0.0001, 0.0000001};
589  static Double_t plow[4] = {-3., -0.05, -3., -0.05};
590  static Double_t phigh[4] = { 3., 0.05, 3., 0.05};
591  minuit->SetParameter(0, "x0" , pstart[0], pstep[0], plow[0], phigh[0]);
592  minuit->SetParameter(1, "dxdz", pstart[1], pstep[1], plow[1], phigh[1]);
593  minuit->SetParameter(2, "y0" , pstart[2], pstep[2], plow[2], phigh[2]);
594  minuit->SetParameter(3, "dydz", pstart[3], pstep[3], plow[3], phigh[3]);
595
596 // Now ready for minimization step
597  double arglist[10];
598  arglist[0] = 500;
599
600  beamWidth = 0;
601  do {
602
603  int status = minuit->ExecuteCommand("MIGRAD", arglist ,1);
604  if (status) {
605  LOG_ERROR << "StVertexMaker: error on migrad call, err = "
606  << status << endm;
607  return;
608  }
609
610  double amin,edm,errdef;
611  int nvpar,nparx;
612  minuit->GetStats(amin,edm,errdef,nvpar,nparx);
613  chi = amin/((double) (nverts-4));
614  LOG_INFO << "beamWidth = " << beamWidth << ", chisq = " << amin << ", chisq/dof = " << chi <<
615  "\n *****************************************************" << endm;
616
617  beamWidth += 0.005 * TMath::Min((int) (chi*0.1+1),10); // 50 micron steps
618  } while (chi>1.1 && beamWidth<=0.15);
619
620  char pname[10];
621  for (int i=0; i<4; i++)
622  minuit->GetParameter(i, pname, p[i], ep[i], plow[i], phigh[i]);
623 }
624 //_____________________________________________________________________________
625 TNtupleD* StVertexSeedMaker::newBLpars() {
626  // up to the user of this function to delete the ntuple
627  return new TNtupleD("BLpars","BeamLine parameters",
628  "days:x0:err_x0:y0:err_y0:dxdz:err_dxdz:dydz:err_dydz:stats:date:fill:zdc:chi:beamwidth");
629 }
630 //_____________________________________________________________________________
631 void StVertexSeedMaker::Packer(int firstbit, int nbits, int& var, unsigned short val) {
632  var = val;
633  // erase nbits of detmap, starting at firstbit,
634  // and fill the bits with val, capped at 2^nbits-1
635  int cap = ~((~0)<<nbits); // e.g. nbits=3 => cap=7
636  (detmap &= ~(cap<<firstbit)) |= (TMath::Min(var,cap)<<firstbit);
637 }
638 //_____________________________________________________________________________
639 int StVertexSeedMaker::Aggregate(char* dir, const char* cuts, const int offset) {
640  // Format of filenames for parsing must be:
641  // vertexseedhist.DDDDDDDD.TTTTTT.root
642  // where D and T are 8 and 6 digit representations of date and time
643
644  // offset will be applied as a time offset in the date.time use in writing new files
645  SetOffset(offset);
646
647  TFile parsOut("BLpars.root","RECREATE");
648  parsNtuple = newBLpars();
649
650  const char* defaultDir = "./";
651  TString dirStr = dir;
652  if (!dir) dirStr = defaultDir;
653  if (dirStr.EndsWith("/")) dirStr.Append("vertexseedhist.*.root");
654  StFile allFiles;
656  // allFiles.ls(); allFiles.Rewind();
657  TList fileList;
658  const char* fileName;
659  for (int fileb = 0; fileb < allFiles.GetNBundles(); fileb++) {
660  allFiles.GetNextBundle();
661  int filen=0;
662  while ((fileName = allFiles.GetFileName(filen++))) {
664  }
665  }
666  fileList.SetOwner();
667  fileList.Sort();
668
669  int aggMode = 0;
670  // 0 : by fill
671  // 1 : all data into one result
672  // 2 : by date (flawed for input files & runs that span dates)
673  // 3 : by file (retains the same number of files)
674
675  // Parse cuts for keywords
676  TString cutsStr = cuts;
677  if (cutsStr.Contains("ALL")) {
678  aggMode = 1;
679  cutsStr.ReplaceAll("ALL&&","");
680  cutsStr.ReplaceAll("&&ALL","");
681  cutsStr.ReplaceAll("ALL","");
682  } else if (cutsStr.Contains("DATE")) {
683  aggMode = 2;
684  cutsStr.ReplaceAll("DATE&&","");
685  cutsStr.ReplaceAll("&&DATE","");
686  cutsStr.ReplaceAll("DATE","");
687  } else if (cutsStr.Contains("FILE")) {
688  aggMode = 3;
689  cutsStr.ReplaceAll("FILE&&","");
690  cutsStr.ReplaceAll("&&FILE","");
691  cutsStr.ReplaceAll("FILE","");
692  }
693  const char* cleanedCuts = cutsStr.Data();
694
695  // Try to catch stuck values and ignore them
696  static float prevX = -987.0;
697
698  TFile* currentFile=0;
699  float* vals=0;
700  int nfiles = fileList.GetSize();
701  for (int filen = 0; filen < nfiles; filen++) {
702  int fillf = fill;
703  int datef = date;
704  int timef = time;
705  fileName = fileList.At(filen)->GetName();
706  if (aggMode != 1 || date==0) {
707  TString dateTime = fileName;
708  dateTime.Remove(0,dateTime.Last('/') + 1);
709  dateTime.Remove(0,dateTime.First('.') + 1).Remove(15);
710  TString dateStr = dateTime;
711  if (aggMode == 3) GetFillDateTime(); // fill date and time will be overwritten
712  date = atoi(dateStr.Remove(8).Data());
713  time = atoi(dateTime.Remove(0,9).Remove(6).Data());
714  if (aggMode != 3) GetFillDateTime(); // file date and time will be overwritten
715  if (aggMode == 2) time = 0;
716  }
717  if ((currentFile) && (
718  (aggMode == 0 && fill != fillf) ||
719  (aggMode == 2 && date != datef) ||
720  (aggMode == 3) )) {
721  LOG_INFO << "Aggregator has changed\n"
722  << " Processing data from previous aggregation before continuing" << endm;
723  int fillp = fill;
724  int datep = date;
725  int timep = time;
726  int timeFillp = timeFill;
727  fill = fillf;
728  date = datef;
729  time = timef;
730  currentFile->Close();
731  currentFile = 0;
732  FindResult(kFALSE);
733  Reset();
734  fill = fillp;
735  date = datep;
736  time = timep;
737  timeFill = timeFillp;
738  }
739
740  LOG_INFO << "Now opening file:\n " << fileName << endm;
741  if (currentFile) currentFile->Close();
742  currentFile = new TFile(fileName);
743  TNtuple* curNtuple = static_cast<TNtuple*>(currentFile->Get("resNtuple"));
744  if (!curNtuple) {
745  LOG_ERROR << "No resNtuple found in " << fileName << endm;
746  continue;
747  }
748  curNtuple->Draw(">>elistVtxSeed",cleanedCuts);
749  TEventList* elist = static_cast<TEventList*>(gDirectory->Get("elistVtxSeed"));
750  int nentries = (elist ? (int) elist->GetN() : 0);
751  int nvar = curNtuple->GetNvar();
752  int rvar = resNtuple->GetNvar();
753  for (int entryn = 0; entryn < nentries; entryn++) {
754  curNtuple->GetEntry(elist->GetEntry(entryn));
755  vals = curNtuple->GetArgs();
756  if (vals[1] == prevX) continue; // stuck value!
757  else prevX = vals[1];
758  unsigned int tid = (unsigned int) vals[5];
759  bool updateForTimeFill = (nvar > 24 && timeFill >=0 &&
760  vals[24] >= 0 && vals[25] < 0);
761  if (ValidTrigger(tid)) {
762  if (nvar < rvar || updateForTimeFill) {
763  float vals2[32];
764  memset(vals2,0,32*sizeof(float));
765  memcpy(vals2,vals,nvar*sizeof(float));
766  if (nvar < 20) {
767  // detmap should be converted...
768  detmap = (int) (vals[12]);
769  bmatch = (detmap)&7;
770  ematch = (detmap>>3)&7;
771  tmatch = (detmap>>6)&7;
772  cmatch = (detmap>>9)&3;
773  hmatch = (detmap>>11)&7;
774  vals2[16] = (float) bmatch;
775  vals2[17] = (float) ematch;
776  vals2[18] = (float) tmatch;
777  vals2[19] = (float) cmatch;
778  vals2[20] = (float) hmatch;
779  }
780  if (nvar < 22) {
781  vals2[21] = (float) pmatch;
782  vals2[22] = (float) pct;
783  vals2[23] = vpd_zvertex;
784  }
785  if (nvar < 25) {
786  vals2[24] = (float) (timeEvent % 86400);
787  vals2[25] = (float) (timeFill >= 0 ? timeEvent - timeFill : timeFill);
788  } else if (updateForTimeFill) {
789  // timeFill not obtained previously, but available now
790  int timeTemp = ((int) vals[24]) - (timeFill % 86400);
791  while (timeTemp < 0) timeTemp += 86400;
792  vals2[25] = (float) timeTemp;
793  }
794  resNtuple->Fill(vals2);
795  } else
796  resNtuple->Fill(vals);
797  if (nvar > 13) // errors are available
799  else
801  sumzdc += vals[8];
802  } else {
803  LOG_INFO << "Invalid trigger: " << tid << endm;
804  }
805  }
806  LOG_INFO << "Current statistics: " << nverts << endm;
807  }
808  if (currentFile) {
809  currentFile->Close();
810  delete currentFile;
811  currentFile = 0;
812  FindResult(kFALSE);
813  }
814  parsOut.Write();
815  parsOut.Close();
816  LOG_INFO << "Examined " << nfiles << " files" << endm;
817  return nfiles;
818 }
819 //_____________________________________________________________________________
820 // \$Id: StVertexSeedMaker.cxx,v 1.64 2017/08/08 03:58:20 genevb Exp \$
821 // \$Log: StVertexSeedMaker.cxx,v \$
822 // Revision 1.64 2017/08/08 03:58:20 genevb
823 // Add vertex-seed-finding with picoDsts
824 //
825 // Revision 1.63 2017/08/07 18:10:57 genevb
826 // Introduce modes for aggregation
827 //
828 // Revision 1.62 2016/08/02 21:17:17 genevb
829 // Added tDay,tFill to resNtuple, and improved C++11 compliance
830 //
831 // Revision 1.61 2015/08/31 19:17:00 genevb
832 // Add chi and beamwidth to BLpars ntuple
833 //
834 // Revision 1.60 2015/05/23 02:38:07 genevb
835 // Ntuples attached to files should not get deleted, more rapid chi2 convergence
836 //
837 // Revision 1.59 2015/05/19 19:36:09 genevb
838 // Code cleanup in preparation for C++11
839 //
840 // Revision 1.58 2015/05/19 17:58:15 genevb
841 // Better garbage collection
842 //
843 // Revision 1.57 2015/05/15 05:38:21 genevb
844 // Include prompt hits and post-crossing tracks, simplify detmap packing, update doxygen documentation
845 //
846 // Revision 1.56 2015/05/14 20:29:25 genevb
847 // Add z of VPD vertex
848 //
849 // Revision 1.55 2013/08/14 21:42:48 genevb
850 // Introduce time offsets, noclobber toggle, more matched-tracks controls
851 //
852 // Revision 1.54 2012/10/11 16:33:12 genevb
853 // Protect against zero entries, and a more unique entry list name
854 //
855 // Revision 1.53 2012/10/01 17:50:07 genevb
856 // Reduce some overhead DB queries by being more specific about needed tables
857 //
858 // Revision 1.52 2012/08/22 04:52:35 genevb
859 // Add BeamLine parameter ntuples to output
860 //
861 // Revision 1.51 2012/08/17 22:57:33 genevb
862 // Add index of vertex within event to ntuple
863 //
864 // Revision 1.50 2012/08/16 05:37:07 genevb
865 // Missing mean ZDC for aggregation
866 //
867 // Revision 1.49 2012/08/15 22:11:06 genevb
868 // Improved doxygen-ready documentation
869 //
870 // Revision 1.48 2012/08/14 23:56:06 genevb
871 // detmap now includes BEMC+EEMC+BTOF+CM, added mean zdc to log output
872 //
873 // Revision 1.47 2012/02/29 02:00:04 genevb
874 // Hack to get TVirtualFitter to allow more than 3 parameters
875 //
876 // Revision 1.46 2012/02/28 21:54:37 genevb
877 // Restore .root to .ROOT to avoid St_db_Maker read-in
878 //
879 // Revision 1.45 2010/07/19 20:48:32 genevb
880 // Forgot to shift ex,ey after detmap addition to ntuple
881 //
882 // Revision 1.44 2010/07/02 22:36:10 genevb
883 // Option for using all triggers
884 //
885 // Revision 1.43 2009/11/23 21:38:56 genevb
886 // Fix problems with memory-resident TNtuple by using a temporary disk file
887 //
888 // Revision 1.42 2009/11/16 22:31:11 genevb
889 // phase out usage of old tables
890 //
891 // Revision 1.41 2009/11/10 20:54:13 fisyak
892 // pams Cleanup
893 //
894 // Revision 1.40 2009/05/22 23:50:50 genevb
895 // Code mods for BEMC matches, BeamWidth
896 //
897 // Revision 1.39 2008/05/21 17:48:39 genevb
898 // Use vertex errors for weighting
899 //
900 // Revision 1.38 2008/04/29 23:30:33 genevb
901 // Added cuts capability to Aggregate
902 //
903 // Revision 1.37 2008/04/29 19:06:06 genevb
904 // handle no DB access
905 //
906 // Revision 1.36 2007/11/27 23:42:47 genevb
907 // Move valid triggers from code to DB
908 //
909 // Revision 1.35 2007/07/12 19:46:56 fisyak
910 // Add includes for ROOT 5.16
911 //
912 // Revision 1.34 2007/04/28 17:56:31 perev
913 // Redundant StChain.h removed
914 //
915 // Revision 1.33 2007/04/22 04:25:59 genevb
916 // printf, gMessMgr ==> STAR Logger
917 //
918 // Revision 1.32 2006/09/01 22:27:16 genevb
919 // More detailed info in ntuple
920 //
921 // Revision 1.31 2006/08/16 21:58:01 genevb
922 // Added 2006 pp62 triggers
923 //
924 // Revision 1.30 2006/05/27 19:54:03 genevb
925 // Yet more 2006 triggers
926 //
927 // Revision 1.29 2006/05/11 18:09:44 genevb
928 // More pp2006 triggers
929 //
930 // Revision 1.28 2006/05/10 03:57:08 genevb
931 // ppProductionTrans trigger for 2006
932 //
933 // Revision 1.27 2006/05/08 02:38:21 genevb
934 // Added 2006 triggers
935 //
936 // Revision 1.26 2005/07/16 21:24:03 genevb
937 // Fixed bug with pp400 data from 2005
938 //
939 // Revision 1.25 2005/07/15 18:38:47 genevb
940 // ppTrans triggers, and fix for too many open files
941 //
942 // Revision 1.24 2005/07/14 21:02:40 genevb
943 // Modified use of test triggers
944 //
945 // Revision 1.23 2005/07/01 21:46:01 genevb
946 // Specify output directory
947 //
948 // Revision 1.22 2005/06/14 18:51:31 genevb
949 // Updates to allow for pp2005 triggers, and inheritance
950 //
951 // Revision 1.21 2004/08/02 01:19:33 genevb
952 // minor fixes for getting directories correct
953 //
954 // Revision 1.20 2004/07/23 16:56:01 genevb
955 // 2004 pp triggers
956 //
957 // Revision 1.19 2003/09/02 17:58:45 perev
958 // gcc 3.2 updates + WarnOff
959 //
960 // Revision 1.18 2003/05/30 19:02:51 genevb
961 // Add ppFPDw-slow
962 //
963 // Revision 1.17 2003/05/19 14:34:54 genevb
964 // Add ppFPD slow triggers
965 //
966 // Revision 1.16 2003/04/22 19:20:52 genevb
967 // Add ppMinBias triggers
968 //
969 // Revision 1.15 2003/03/21 15:12:24 genevb
970 // Allow use of TOF triggers
971 //
972 // Revision 1.14 2003/02/25 03:49:16 genevb
973 // Choose only dAu minbias triggers
974 //
975 // Revision 1.13 2003/02/12 04:19:39 genevb
976 // Small improvements
977 //
978 // Revision 1.12 2003/02/11 22:24:20 genevb
979 // Update to use beamInfo table from database
980 //
981 // Revision 1.11 2003/01/22 23:32:23 genevb
982 // Type cast fix
983 //
984 // Revision 1.10 2002/04/19 22:24:16 perev
985 // fixes for ROOT/3.02.07
986 //
987 // Revision 1.9 2002/03/23 01:05:15 jeromel
988 // Create files with extension .ROOT instead of .root (db_Maker will read the .root
989 // and crash otherwise). Will fix this with a more elegant scheme later.
990 //
991 // Revision 1.8 2002/03/23 00:23:54 genevb
992 // switch from float arrays to TArrayF
993 //
994 // Revision 1.7 2002/03/22 20:29:09 jeromel
995 // Not .C but .root
996 //
997 // Revision 1.6 2002/03/20 21:09:10 genevb
998 // Oops...introduced a typo last time
999 //
1000 // Revision 1.5 2002/03/20 19:22:59 genevb
1001 // changed output directory of histogram files
1002 //
1003 // Revision 1.4 2002/03/20 00:40:42 genevb
1004 // Addition of Aggregate feature, minor updates
1005 //
1006 // Revision 1.3 2002/01/28 22:16:33 genevb
1007 // Some revisions to output date/time stamps
1008 //
1009 // Revision 1.2 2002/01/27 01:56:14 genevb
1010 // Write db files with date/time for fill
1011 //
1012 // Revision 1.1 2002/01/26 18:55:33 jeromel
1013 // StTpcT0Maker moved from directory of the same name. First version
1014 // of StVertexSeedMaker.
1015 //
1016 //
Definition: StTree.h:125
virtual void Remove(TDataSet *set)
Remiove the &quot;set&quot; from this TDataSet.
Definition: TDataSet.cxx:641
virtual void Clear(Option_t *option)
User defined functions.
virtual Int_t Make()
static StDbManager * Instance()
strdup(..) is not ANSI
Definition: StDbManager.cc:155
virtual Int_t Finish()
Definition: StMaker.cxx:776
BeamLine Constraint calibration base class.
Definition: Stypes.h:44
virtual Int_t Finish()
Definition: Stypes.h:41