StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
calibFms.C
1 /*
2  FMS calibration macro: perform fit on each mass, update gainCorr, and print new gainCorr set
3 
4  Chong Kim
5  UC Riverside
6  ckim@bnl.gov
7 */
8 
9 #include "calibFmsShow.C" //calibFmsShow.C includes calibFmsTools.C dlready
10 
11 void calibFms(const int iIter = 0, const char* outPath = "./Iterations", bool PRINT = true)
12 {
13  TStopwatch SW;
14  SW.Start();
15 
16  //TFile* F = TFile::Open("./fmsCalib.root");
17  TFile* F = TFile::Open(Form("./fmsCalib_%i.root", iIter));
18  TFile* G = new TFile(Form("out_fmsCalib_%i.root", iIter), "RECREATE");
19 
20  TH2F* H2_mass[4];
21  for (int a=0; a<4; a++) H2_mass[a] = (TH2F*)F->Get(Form("mass_d%i", a+8));
22 
23  //Map with exclamation mark is essentially required to proceed
24  map<int, int> iToCh[4];
25  map<int, int> chToBS[4];
26  map<int, int> chToCellStat[4];
27  map<int, float> chToGainCorr[4];
28  map<int, st_pos> chToPos[4];
29  GetMapIndexToCh ("./FmsMapBase.txt", iToCh);
30  GetMapChToBS ("./FmsMapBitShift.txt", chToBS);
31  GetMapChToCellStat("./FmsCellStat.txt", chToCellStat);
32  GetMapChToGainCorr("./FmsGainCorr.txt", chToGainCorr);
33  GetMapChToPos ("./FmsMapBase.txt", chToPos);
34 
35  //Perform fit, Update gainCorr
36  //--------------------------------------------------------------------
37 
38  //Essential
39  map<int, st_fitR> chToFitR[4];
40  map<int, float> chToGainCorrNext[4]; //Results of "current" iteration: will be used in upcoming iteration
41 
42  //Watchdog + QA
43  map<int, int> chToSkew[4];
44  map<int, st_fitR> chToFitRPrev[4];
45  if (iIter > 0) GetMapChToSkew(Form("%s/out_fmsCalibSkew_%i.txt", outPath, iIter-1), chToSkew);
46  if (iIter > 0) GetMapChToFitR(Form("%s/out_fmsCalibFit_%i.txt", outPath, iIter-1), chToFitRPrev);
47 
48  //In-situ manual gainCorr update: enforce values in this list for next iteration
49  map<int, float> manGainCorr[4];
50  GetMapManualGainCorr(Form("FmsGainCorrManual_%i.txt", iIter), manGainCorr);
51 
52  //In-situ QA lists
53  vector<int> markBadF[4]; //Bad fit
54  vector<int> markCold[4]; //Not bad nor dead, but has too small # of entries (< 100)
55  vector<int> markLoHi[4]; //Too low (<0.5) or Too high (>3.0) gainCorr
56  vector<int> markSkew[4]; //Channels judged skewed
57 
58  #if 1
59  for (int a=0; a<4; a++)
60  {
61  //Perform fit, get position of each cell's mass
62  FitMass(iIter, H2_mass[a], iToCh[a], chToBS[a], chToCellStat[a], chToGainCorr[a], chToPos[a], chToFitR[a]);
63 
64  //Update gainCorr
65  for (unsigned int b=0; b<iToCh[a].size(); b++)
66  {
67  const int ch = iToCh[a][b];
68  const float gainCorr = chToGainCorr[a][ch]; //gainCorr applied for this iteration
69  const float massFitR = chToFitR[a][ch].mass;
70  const float massBook = 0.1349766;
71  int nSkewed = (iIter==0)?0:chToSkew[a][ch];
72 
73  bool doUpdate = false;
74  if (chToCellStat[a][ch]==GOOD)
75  {
76  //Return mass will be 0 if # of entries is too small
77  if (massFitR != 0.) doUpdate = true;
78  else markCold[a].push_back(ch);
79  }
80 
81  const char* infoQA = "";
82  float gainCorrNext = gainCorr;
83  if (doUpdate)
84  {
85  gainCorrNext = (massBook/massFitR) * gainCorr;
86 
87  //Check too low or too high gainCorr
88  if (gainCorr<0.3 || gainCorr>3.0) markLoHi[a].push_back(ch);
89 
90  //Check quality of the fit, limit ratio of allowed gainCorr update by 20 % if quality looks bad
91  if (fabs(massFitR - massBook) > 0.03 ||
92  chToFitR[a][ch].nRefit == 10 ||
93  chToFitR[a][ch].chi2 > 25. ||
94  chToFitR[a][ch].chi2 < 0.5)
95  {
96  const float upRatio = (gainCorrNext - gainCorr) / gainCorr;
97  if (fabs(upRatio)>0.2 && upRatio>0) gainCorrNext = gainCorr * 1.2;
98  if (fabs(upRatio)>0.2 && upRatio<0) gainCorrNext = gainCorr * 0.8;
99  infoQA = "BADFIT ";
100  markBadF[a].push_back(ch);
101  }
102 
103  //Check skewed mass: mass position don't move even if higher gainCorr applied, small cells
104  if (chToCellStat[a][ch]!=CONVERGED && iIter>0 && a>1 && massFitR<0.12)
105  {
106  const float dGainCorr = gainCorrNext - gainCorr;
107  const float dMassFitR = massFitR - chToFitRPrev[a][ch].mass;
108  if (dGainCorr>=0 && dMassFitR<0) //gainCorr increased but mass decreased: cancel update
109  {
110  gainCorrNext = gainCorr;
111  infoQA = Form("%sSKEWED ", infoQA);
112  markSkew[a].push_back(ch);
113 
114  nSkewed++;
115  if (nSkewed==3)
116  {
117  infoQA = Form("%sRESET! ", infoQA);
118  gainCorrNext = 0.750;
119  nSkewed = 0;
120  }
121  map<int, int>::iterator chNewSkew = chToSkew[a].find(ch);
122  if (chNewSkew != chToSkew[a].end()) chNewSkew->second = nSkewed;
123  }
124  }
125  }//doUpdate
126 
127  //In-situ manual gainCorr update by user, between iterations
128  if (manGainCorr[a].size()!=0 && manGainCorr[a][ch]!=0.)
129  {
130  infoQA = Form("%sMANUAL", infoQA);
131  gainCorrNext = manGainCorr[a][ch];
132  }
133 
134  //Final updated gainCorr list for next iteration
135  chToGainCorrNext[a].insert(pair<int, float>(ch, gainCorrNext));
136 
137  if (PRINT)
138  {
139  const char* infoUp = (gainCorrNext != gainCorr)?Form("-> %4.3f", gainCorrNext):Form("%8s", "");
140  cout <<Form("%2i %3i | %4.3f %s %s", a+8, ch, gainCorr, infoUp, infoQA) <<endl;
141  }
142  }//b, ch
143  }//a, detId
144  F->Close();
145  G->Close();
146  #endif
147 
148  //Print QA plots for in-situ monitoring
149  //-------------------------------------
150 
151  #if 1
152  cout <<"Generating in-situ QA plots..." <<endl;
153  DrawMap(iToCh, chToPos, Form("GainCorrL_i%i", iIter), false, false, true, chToGainCorr, chToCellStat);
154  DrawMap(iToCh, chToPos, Form("GainCorrS_i%i", iIter), true, false, true, chToGainCorr, chToCellStat);
155 
156  bool QaBadF = false;
157  bool QaCold = false;
158  bool QaLoHi = false;
159  bool QaSkew = false;
160  for (int a=0; a<4; a++) { if (markBadF[a].size() > 0) { QaBadF = true; break; } }
161  for (int a=0; a<4; a++) { if (markCold[a].size() > 0) { QaCold = true; break; } }
162  for (int a=0; a<4; a++) { if (markLoHi[a].size() > 0) { QaLoHi = true; break; } }
163  for (int a=0; a<4; a++) { if (markSkew[a].size() > 0) { QaSkew = true; break; } }
164  if (QaBadF) DrawMap(iToCh, chToPos, "QABadF", false, false, true, DumMapF, chToCellStat, markBadF);
165  if (QaCold) DrawMap(iToCh, chToPos, "QACold", false, false, true, DumMapF, chToCellStat, markCold);
166  if (QaLoHi) DrawMap(iToCh, chToPos, "QALoHi", false, false, true, DumMapF, chToCellStat, markLoHi);
167  if (QaSkew) DrawMap(iToCh, chToPos, "QASkew", false, false, true, DumMapF, chToCellStat, markSkew);
168  if (QaBadF || QaCold || QaLoHi || QaSkew)
169  {
170  const char* QAout = Form("QAMap_i%i", iIter);
171  gSystem->Exec(Form("mkdir %s", QAout));
172  gSystem->Exec(Form("mv *.png %s", QAout));
173  gSystem->Exec(Form("mv %s %s", QAout, outPath));
174  }
175  #endif
176 
177  //Print files for next iteration + Test convergence
178  //-------------------------------------------------
179 
180  #if 1
181  PrintFitResults(Form("out_fmsCalibFit_%i.txt", iIter), iToCh, chToFitR);
182  gSystem->Exec(Form("mv out_fmsCalibFit_%i.txt %s", iIter, outPath));
183 
184  PrintGainCorr(Form("FmsGainCorrNew_%i.txt", iIter), chToGainCorrNext);
185  gSystem->Exec(Form("mv FmsGainCorr.txt %s/FmsGainCorr_%i.txt", outPath, iIter));
186  gSystem->Exec(Form("mv FmsGainCorrNew_%i.txt FmsGainCorr.txt", iIter));
187 
188  PrintSkewMonitor(Form("out_fmsCalibSkew_%i.txt", iIter), iToCh, chToSkew);
189  gSystem->Exec(Form("mv out_fmsCalibSkew_%i.txt %s", iIter, outPath));
190 
191  TestConvergence(iIter, iToCh, chToCellStat, outPath); //Includes PrintCellStat
192  #endif
193 
194  SW.Stop();
195  SW.Print();
196  return;
197 }//Main
Generated by   doxygen 1.8.5