CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HIPplots.cc
Go to the documentation of this file.
1 #include "TROOT.h"
2 #include "TAttFill.h"
3 #include "TColor.h"
4 #include "HIPplots.h"
5 #include <string>
6 #include <sstream>
7 #include <cmath>
8 
9 #include "TProfile.h"
10 #include "TPaveStats.h"
11 #include "TList.h"
12 #include "TNtuple.h"
13 #include "TString.h"
14 #include "TTree.h"
15 
16 #include "TMatrixD.h"
17 #include "TVectorD.h"
18 #include "TLegend.h"
19 
20 HIPplots::HIPplots(int IOV, char* path, char* outFile ){
21 
22  sprintf( _path, path );
23  sprintf( _outFile, outFile );
24  sprintf( _inFile_params, "%s/IOAlignmentParameters.root", _path ); //Id (raw id), ObjId, Par (alignable movements), CovariantMatrix
25  sprintf( _inFile_uservars, "%s/IOUserVariables.root", _path ); //Id,ObjId,Nhit,Jtvj,Jtve,AlignableChi2,AlignableNdof
26  sprintf( _inFile_truepos, "%s/IOTruePositions.root", _path ); //redundant
27  sprintf( _inFile_alipos, "%s/IOAlignedPositions_%d.root", _path,IOV ); //aligned absolute positions
28 // sprintf( _inFile_alipos, "%s/IOAlignedPositions_262922.root", _path ); //aligned absolute positions
29 // sprintf( _inFile_alipos, "%s/IOAlignedPositions_272497.root", _path ); //aligned absolute positions
30  sprintf( _inFile_mispos, "%s/IOMisalignedPositions.root", _path );
31  sprintf( _inFile_HIPalign, "%s/HIPAlignmentAlignables.root", _path );
32  sprintf( _inFile_surveys, "%s/HIPSurveyResiduals.root", _path );
33  SetPeakThreshold(8.0);
34  plotbadchi2=true;
35 }
36 
37 
38 TLegend * HIPplots::MakeLegend (double x1,
39  double y1,
40  double x2,
41  double y2)
42 {
43  TLegend * legend = new TLegend (x1, y1, x2, y2, "", "NBNDC");
44  legend->SetNColumns(6);
45  legend->SetFillColor(0);
46  legend->SetBorderSize(0);
47  int COLOR_CODE[6]={28,2,3,4,6,7};
48 
49  TString detNames[6];
50  detNames[0] = TString("PXB");
51  detNames[1] = TString("PXF");
52  detNames[2] = TString("TIB");
53  detNames[3] = TString("TID");
54  detNames[4] = TString("TOB");
55  detNames[5] = TString("TEC");
56 
57  for (unsigned int isublevel = 0 ; isublevel < 6 ; isublevel++)
58  {
59  TGraph * g = new TGraph (0);
60  g->SetLineColor(COLOR_CODE[isublevel]);
61  g->SetMarkerColor(COLOR_CODE[isublevel]);
62  g->SetFillColor(COLOR_CODE[isublevel]);
63  g->SetFillColor(1);
64  g->SetMarkerStyle(8);
65  g->SetMarkerSize(1);
66  legend->AddEntry(g,detNames[isublevel], "lp");
67  }
68  return legend;
69 }
70 
71 
72 void HIPplots::extractAlignParams(int currentPar, int minHits, int subDet, int doubleSided ){
73 
74  cout<<"--- extracting AlignParams ; Par "<<currentPar<<" ---"<<endl<<endl; //0-5: u,v,w,alpha,beta,gamma
75  //void ExtractAlignPars(int pp, string SubDetString){
76  // const int par_now=pp;
77  //load first tree
78 
79 // TFile *fp = new TFile( _inFile_params,"READ");
80 // const TList* keysp = fp->GetListOfKeys();
81 // const unsigned int maxIteration = keysp->GetSize() - 1;
82  cout<<"Loaded par file -. OK"<<endl;
83  TFile *fv = new TFile( _inFile_uservars,"READ");
84  const TList* keysv = fv->GetListOfKeys();
85  const unsigned int maxIteration = keysv->GetSize() - 1;
86 
87  char fileaction[16];
88  if(CheckFileExistence(_outFile))sprintf(fileaction,"UPDATE");
89  else sprintf(fileaction,"NEW");
90 
91  TFile *fout = new TFile( _outFile, fileaction);
92 
93  TTree* tree0 = (TTree*) fv->Get(keysv->At(0)->GetName() );
94  //const char* tree0v_Name = keysv->At(0)->GetName();
95  //tree0->AddFriend( tree0v_Name, fv );
96 
97  int nHit;
98  unsigned int detId;
99  double par[6];
100  unsigned int detId0;
101 
102  tree0->SetBranchAddress("Id", &detId0 );
103 
104  const int ndets=tree0->GetEntries();
105  TH1D *hpar1[ndets];
106  char ppdirname[16];
107  sprintf(ppdirname,"ShiftsPar%d", currentPar);
108  fout->mkdir(ppdirname);
109  gDirectory->cd(ppdirname);
110  // delare histos
111 
112 
113  TH1D *hpariter[maxIteration];
114  for(unsigned int a = 0; a < maxIteration; a++){
115  char histoname[32],histotitle[32];
116  sprintf(histoname,"Par_%d_Iter_%d",currentPar,a);
117  sprintf(histotitle,"Par %d for iteration #%d",currentPar,a);
118  hpariter[a]=new TH1D(histoname,histoname,1000,-1.0,1.0);
119  }
120 
121  for(int i = 0; i < ndets; i++){
122  char histoname[32],histotitle[32];
123  sprintf(histoname,"Par_%d_SiDet_%d",currentPar,i);
124  sprintf(histotitle,"Par %d for detector #%d",currentPar,i);
125  hpar1[i]=new TH1D(histoname,histoname,maxIteration+1,-0.5,maxIteration+0.5);
126  }
127 
128  //file where to save aligned det id list
129  ofstream flist("./List_aligned_dets.txt",ios::out);
130 
131 
132  //Loop on trees
133  int modules_accepted = 0;
134  for(unsigned int iter = 0; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
135 
136 // TTree* tmpTree = (TTree*) fp->Get(keysp->At(iter)->GetName() );
137 // const char* tmpTreev = keysv->At(iter)->GetName();
138  TTree* tmpTree = (TTree*) fv->Get(keysv->At(iter)->GetName() );
139 // tmpTree->AddFriend( tmpTreev, fv ); // get access to both IOAlignmentParameters and IOUserVariables
140  tmpTree->SetBranchAddress("Id", &detId );
141  tmpTree->SetBranchAddress("Nhit", &nHit ); // taken from IOUserVariables
142  tmpTree->SetBranchAddress("Par", &par ); // taken from IOAlignmentParameters
143 
144  std::cout << "iteration: " << iter << "..." << std::endl;
145 
146  modules_accepted=0;
147  //loop on entries
148  for(int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
149 
150  tmpTree->GetEntry(j);
151  bool passSubdetCut = true;
152  if (subDet > 0){ if ( GetSubDet(detId) != subDet ) passSubdetCut = false; } // find the subDet from module detId
153  if (doubleSided > 0){
154  if ( GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3 ) passSubdetCut = false;
155  if ( GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2 ) passSubdetCut = false;
156  }
157 
158  if (( nHit >= minHits )&&( passSubdetCut )){ //cut on min nhits per module
159  hpar1[j]->SetBinContent(iter+1,par[currentPar]); //x-iteration, y-movement in this iteration
160  //std::cout << "iteration: " << iter << ",par"<<currentPar<<"= "<<par[currentPar] << std::endl;
161 
162  int COLOR_CODE[6]={28,2,3,4,6,7};
163  hpar1[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
164  hpar1[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
165 
166  // hpariter[iter]->Fill(par[currentPar]);
167  //if(currentPar<3&&par[currentPar]>0.5)cout<<"LARGE PARCHANGE from DET "<<detId<<" in subdet "<< subDet<<". Par#"<<currentPar<<" shifted by "<<par[currentPar]<<" at iter "<<iter<<endl;
168 
169  if((iter==maxIteration-1)&&(currentPar==0))flist<<detId<<endl;
170  modules_accepted++;
171 
172  }
173  }//end loop on entries (i.e. modules)
174  delete tmpTree;
175  }//end loop on iteration
176  std::cout << "Modules accepted: " << modules_accepted << std::endl;
177  std::cout << "Writing..." << std::endl;
178  //Save
179  fout->Write();
180  flist.close();
181  // delete hpar1;
182  std::cout << "Deleting..." << std::endl;
183  delete fout;
184  // std::cout << "Deleting..." << std::endl;
185 // delete fp;
186  //std::cout << "Deleting..." << std::endl;
187  delete fv;
188 }
189 
190 void HIPplots::extractAlignShifts( int currentPar, int minHits, int subDet ){
191 
192  cout<<"\n--- extracting AlignShifts ; Par "<<currentPar<<" ---"<<endl<<endl;
193  TFile *fa = new TFile( _inFile_alipos,"READ");
194  const TList* keysa = fa->GetListOfKeys();
195  const unsigned int maxIteration = keysa->GetSize();
196 
197  TFile *fv = new TFile( _inFile_uservars,"READ");
198  const TList* keysv = fv->GetListOfKeys();
199 
200 // TFile *fm = new TFile( _inFile_mispos,"READ");
201 // TFile *ft = new TFile( _inFile_truepos,"READ");
202 
203  char fileaction[16];
204  if(CheckFileExistence(_outFile))sprintf(fileaction,"UPDATE");
205  else sprintf(fileaction,"NEW");
206  TFile *fout = new TFile( _outFile, fileaction);
207 
208  TTree* tree0 = (TTree*) fa->Get(keysa->At(0)->GetName() );
209  const char* tree0v_Name = keysv->At(0)->GetName();
210  tree0->AddFriend( tree0v_Name, fv ); //get access to both IOAlignedPositions and IOUserVariables
211 
212  unsigned int detId0;
213 
214  tree0->SetBranchAddress("Id", &detId0 );
215 
216  const int ndets=tree0->GetEntries();
217  TH1D *hshift[ndets]; //hshift[i] absolute position change for det i, at all iterations
218 // TH1D *hiter_mis[maxIteration]; //hiter_mis[i] misalignment position change at iteration i (defualt case : misalignment=0)
219  TH1D *hiter_ali[maxIteration]; //hiter_ali[i] absolute position change at iteration i, for all detectors
220  char ppdirname[16];
221  sprintf(ppdirname,"Shifts%d", currentPar);
222  fout->mkdir(ppdirname);
223  gDirectory->cd(ppdirname);
224  // delare histos
225  for(unsigned int a = 0; a < maxIteration; a++){
226  char histoname[64],histotitle[64];
227 // sprintf(histoname,"MisalignedShift_%d_Iter_%d",currentPar,a);
228 // sprintf(histotitle,"Misaligned Shift %d for iteration #%d",currentPar,a);
229 // hiter_mis[a]=new TH1D(histoname,histoname,1000,-1.0,1.0);
230  sprintf(histoname,"AlignedShift_%d_Iter_%d",currentPar,a);
231  sprintf(histotitle,"Aligned Shift %d for iteration #%d",currentPar,a);
232  hiter_ali[a]=new TH1D(histoname,histoname,ndets,0,ndets);
233  }
234  for(int i = 0; i < ndets; i++){
235  char histoname[64],histotitle[64];
236  sprintf(histoname,"Shift_%d_SiDet_%d",currentPar,i);
237  sprintf(histotitle,"Shift %d for detector #%d",currentPar,i);
238  hshift[i]=new TH1D(histoname,histoname,maxIteration,-0.5,maxIteration-0.5);
239  }
240 
241  //Loop on trees
242  int modules_accepted = 0;
243  int nHit;
244  unsigned int detId;
245  double tpos[3];
246  double apos[3];
247 // double mpos[3];
248  double trot[9];
249  double arot[9];
250 // double mrot[9];
251  double aerr[6];
252 
253  TString detNames[6];
254  detNames[0] = TString("PXB");
255  detNames[1] = TString("PXF");
256  detNames[2] = TString("TIB");
257  detNames[3] = TString("TID");
258  detNames[4] = TString("TOB");
259  detNames[5] = TString("TEC");
260 
261 
262 // TTree* tmpTree_m = (TTree*)fm->Get("AlignablesAbsPos_1"); //misaligned position
263 // TTree* tmpTree_true = (TTree*)ft->Get("AlignablesOrgPos_1");
264  TTree* tmpTree_true = (TTree*)fa->Get("AlignablesAbsPos_0"); //starting position
265 
266  if (currentPar < 3){
267 // tmpTree_m->SetBranchAddress("Pos", &mpos );
268  tmpTree_true->SetBranchAddress("Pos", &tpos );
269  }
270  if (currentPar >= 3 ){
271 // tmpTree_m->SetBranchAddress("Rot", &mrot );
272  tmpTree_true->SetBranchAddress("Rot", &trot );
273  }
274 
275  for(unsigned int iter = 1; iter < maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
276 
277  TTree* tmpTree = (TTree*) fa->Get(keysa->At(iter)->GetName() );
278  const char* tmpTreev = keysv->At(iter)->GetName();
279  tmpTree->AddFriend( tmpTreev, fv );
280  tmpTree->SetBranchAddress("Id", &detId );
281  tmpTree->SetBranchAddress("Nhit", &nHit );
282  tmpTree->SetBranchAddress("ParError", &aerr );
283 
284  if (currentPar < 3){ tmpTree->SetBranchAddress("Pos", &apos ); }
285  tmpTree->SetBranchAddress("Rot", &arot );
286 
287  //std::cout << "iteration: " << iter << "..." << std::endl;
288  modules_accepted=0;
289  //loop on entries
290  for(int j = 0; j < tmpTree->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
291  tmpTree_true->GetEntry(j);
292  tmpTree->GetEntry(j);
293  //cout<<"det"<<j<<","<<detId<<","<<nHit<<endl;
294 // tmpTree_m->GetEntry(j);
295 // tmpTree_true->GetEntry(j);
296 
297 
298  bool passSubdetCut = true;
299  int mysubDet=GetSubDet(detId);
300  if (subDet > 0){ if ( GetSubDet(detId) != subDet ) passSubdetCut = false; }
301 
302  int COLOR_CODE[6]={28,2,3,4,6,7};
303  hshift[j]->SetLineColor(COLOR_CODE[GetSubDet(detId)-1]);
304  hshift[j]->SetMarkerColor(COLOR_CODE[GetSubDet(detId)-1]);
305 
306  if (( nHit >= minHits )&&( passSubdetCut )){
307  // maths
308  if (currentPar < 3){
309 // TVectorD dr_mis(3, mpos);
310  TVectorD dr_ali(3, apos);
311  TVectorD r_true(3, tpos);
312  TMatrixD R_true(3, 3, trot);
313  //std::cout << "dr_ali 00 : " << dr_ali[currentPar] << std::endl;
314 // dr_mis -= r_true;
315  dr_ali -= r_true;
316  //std::cout << "dr_ali 0 : " << dr_ali[currentPar] << std::endl;
317  // to local
318  //if (dr_mis != 0.) dr_mis = R_true * dr_mis;
319  //if (dr_ali != 0.) dr_ali = R_true * dr_ali;
320  //std::cout << "currentPar: " << currentPar << std::endl;
321  //std::cout << "dr_mis: " << dr_mis[currentPar] << std::endl;
322  //std::cout << "dr_ali: " << dr_ali[currentPar] << std::endl;
323  // if(currentPar<3&&dr_ali[currentPar]>1.0)cout<<"LARGE SHIFT for DET "<<detId<<" in subdet "<< mysubDet<<". Par#"<<currentPar<<" shifted by "<<dr_ali[currentPar]<<" at iter "<<iter<<endl;
324  hshift[j]->SetBinContent(iter+1,dr_ali[currentPar]); //aligned position - start position
325 // hiter_mis[iter]->SetBinContent(j+1,dr_mis[currentPar]);
326  //if(j=0&&currentPar==5) std::cout << "iter="<<iter<<",dr_ali: " << dr_ali[currentPar] << std::endl;
327 // cout << "iter="<<iter<<"dr_ali: " << dr_ali[currentPar] << endl;
328  hiter_ali[iter]->SetBinContent(j+1,dr_ali[currentPar]);
329 // cout << "bin: " <<hiter_ali[iter]->GetBinContent(j+1)<<endl;
330 // hiter_ali[iter]->SetBinError(j+1,5);
331  }
332  if (currentPar >= 3){
333 // TMatrixD dR_mis(3, 3, mrot);
334  TMatrixD dR_ali(3, 3, arot);
335  TMatrixD R_true(3, 3, trot);
336 // dR_mis = dR_mis * TMatrixD(TMatrixD::kTransposed, R_true);
337  dR_ali = dR_ali * TMatrixD(TMatrixD::kTransposed, R_true);
338  //-std::atan2(dR(2, 1), dR(2, 2)), std::asin(dR(2, 0)), -std::atan2(dR(1, 0), dR(0, 0)));
339 // double dR_mis_euler = 0;
340  double dR_ali_euler = 0;
341  if (currentPar == 3 ){
342 // dR_mis_euler = -std::atan2(dR_mis(2, 1), dR_mis(2, 2));
343  dR_ali_euler = -std::atan2(dR_ali(2, 1), dR_ali(2, 2));
344  }
345  if (currentPar == 4 ){
346 // dR_mis_euler = -std::asin(dR_mis(2, 0));
347  dR_ali_euler = -std::asin(dR_ali(2, 0));
348  }
349  if (currentPar == 5 ){
350 // dR_mis_euler = -std::atan2(dR_mis(1, 0), dR_mis(0, 0));
351  dR_ali_euler = -std::atan2(dR_ali(1, 0), dR_ali(0, 0));
352  }
353  hshift[j]->SetBinContent(iter+1,dR_ali_euler);
354 // hiter_mis[iter]->SetBinContent(j+1,dR_mis_euler);
355  hiter_ali[iter]->SetBinContent(j+1,dR_ali_euler);
356  }
357  modules_accepted++;
358  hiter_ali[iter]->GetXaxis()->SetBinLabel(j+1,detNames[GetSubDet(detId)-1]);
359  hiter_ali[iter]->SetBinError(j+1,aerr[currentPar]);
360  //hiter_ali[iter]->SetBinError(j+1,0.001);
361  hiter_ali[iter]->LabelsOption("u","x");
362  }
363  }
364 // delete tmpTree;
365  }
366  delete tmpTree_true;
367 
368  std::cout <<"Modules accepted: " << modules_accepted << std::endl;
369  std::cout << "Writing..." << std::endl;
370  //Save
371  fout->Write();
372  // delete hpar1;
373  std::cout << "Deleting..." << std::flush;
374  delete fout;
375  std::cout << "Deleting..." << std::flush;
376  delete fa;
377 // delete ft;
378 // delete fm;
379  std::cout << "Deleting..." << std::endl;
380  delete fv;
381 
382 
383 }
384 
385 void HIPplots::plotAlignParams(string ShiftsOrParams, char* plotName){
386 
387 
388  cout<<"_|_| plotting AlignParams |_|_"<<endl<<"---> "<<ShiftsOrParams <<endl;
389  bool bParams = false;
390  bool bShifts = false;
391  if (ShiftsOrParams == "PARAMS") bParams = true;
392  if (ShiftsOrParams == "SHIFTS") bShifts = true;
393 
394  int i = 0;
395 
396  TFile *f = new TFile( _outFile,"READ" );
397  // f->ls();
398 
399  TCanvas *c_params = new TCanvas("can_params","CAN_PARAMS",1200,900);
400  c_params->Divide(3,2);
401  // cout<<"(1) I am in "<<gDirectory->GetPath()<<endl;
402  TDirectory* d;
403  int ndets = 0;
404  if (bParams){
405  d = (TDirectory*)f->Get("ShiftsPar0");
406  ndets = GetNIterations(d,"Par_0_SiDet_");
407  // std::cout << "INHERE!" << std::endl;
408  }
409  if (bShifts){
410  d = (TDirectory*)f->Get("Shifts0");
411  ndets = GetNIterations(d,"Shift_0_SiDet_");
412  // std::cout << "INHERE!" << std::endl;
413  }
414 
415  for(int iPar = 0; iPar < 6; iPar++){
416 
417  c_params->cd( iPar+1 );
418  gPad->SetTopMargin(0.15);
419  gPad->SetBottomMargin(0.15);
420  char ppdirname[32];
421  if (bParams) sprintf( ppdirname,"ShiftsPar%d", iPar );
422  if (bShifts) sprintf( ppdirname,"Shifts%d", iPar );
423  if(iPar > 0)gDirectory->cd("../");
424  gDirectory->cd(ppdirname);
425  // cout<<"(2) I am in "<<gDirectory->GetPath()<<endl;
426 
427  TH1D *hpar1[ndets];
428  char histoname[16];
429  int sampling_ratio=1;
430  int ndets_plotted=(int)ndets/sampling_ratio;
431  cout<<"Plotting "<<ndets_plotted<<" detectors over a total of "<<ndets<<endl;
432  i=0;
433  double histomax, histomin;
434  if(iPar>=3){ // alpha, beta, gamma
435  histomax=0.5;// in mrad
436  histomin=-0.5;
437  }
438  else if (iPar>=2) { // w
439  if(bShifts){
440  histomax=200.0;//in microns
441  histomin=-200.0;
442  }
443  else{
444  histomax=100.0;//in microns
445  histomin=-100.0;
446  }
447  }
448  else { // u, v
449  if(bShifts){
450  histomax=200.0;//in microns
451  histomin=-200.0;
452  }
453  else{
454  histomax=100.0;//in microns
455  histomin=-100.0;
456  }
457  }
458  while((i<ndets_plotted) && (i*sampling_ratio<ndets) ){
459  if (bParams) sprintf(histoname,"Par_%d_SiDet_%d",iPar,i*sampling_ratio);
460  if (bShifts) sprintf(histoname,"Shift_%d_SiDet_%d",iPar,i*sampling_ratio);
461  hpar1[i]=(TH1D*)gDirectory->Get(histoname);
462 
463  if(iPar>=3)hpar1[i]->Scale(1000.0); //convert from rad to mrad
464  else hpar1[i]->Scale(10000.0); //convert from cm to um
465 
466  hpar1[i]->SetMarkerStyle(7);
467  hpar1[i]->SetStats(0);
468 
469  double tmpmax=hpar1[i]->GetBinContent(hpar1[i]->GetMaximumBin());
470  double tmpmin=hpar1[i]->GetBinContent(hpar1[i]->GetMinimumBin());
471 
472  //Auto-adjust axis range
473  if(tmpmax>histomax)histomax=tmpmax*1.2;
474  if(tmpmin<histomin)histomin=tmpmin*1.2;
475 
476  //if(i%1300==0)cout<<"Actual maximum is "<<histomax<<" Actual minimum is "<<histomin<<endl;
477  if(i==0){
478 
479  hpar1[i]->SetXTitle("Iteration");
480 
481  if (bParams){
482  if(iPar==0)hpar1[i]->SetYTitle("#delta u param (#mum)");
483  else if(iPar==1)hpar1[i]->SetYTitle("#delta v param (#mum)");
484  else if(iPar==2)hpar1[i]->SetYTitle("#delta w param (#mum)");
485  else if(iPar==3)hpar1[i]->SetYTitle("#delta #alpha param (mrad)");
486  else if(iPar==4)hpar1[i]->SetYTitle("#delta #beta param (mrad)");
487  else if(iPar==5)hpar1[i]->SetYTitle("#delta #gamma param (mrad)");
488  else hpar1[i]->SetYTitle("dunno");}
489  else if (bShifts){
490  if(iPar==0)hpar1[i]->SetYTitle("#delta u shift (#mum)");
491  else if(iPar==1)hpar1[i]->SetYTitle("#delta v shift (#mum)");
492  else if(iPar==2)hpar1[i]->SetYTitle("#delta w shift (#mum)");
493  else if(iPar==3)hpar1[i]->SetYTitle("#delta #alpha shift (mrad)");
494  else if(iPar==4)hpar1[i]->SetYTitle("#delta #beta shift (mrad)");
495  else if(iPar==5)hpar1[i]->SetYTitle("#delta #gamma shift (mrad)");
496  else hpar1[i]->SetYTitle("dunno");}
497 
498  hpar1[i]->GetYaxis()->SetTitleOffset(1.5);
499  // hpar1[i]->SetTitle("Par1: x shift");
500  hpar1[i]->SetTitle("");
501  hpar1[i]->SetMaximum(histomax);
502  hpar1[i]->SetMinimum(histomin);
503  hpar1[i]->Draw("PL");
504  }
505 
506  else hpar1[i]->Draw("PLsame");
507  i++;
508  }//end loop on NDets
509  hpar1[0]->SetMaximum(histomax);
510  hpar1[0]->SetMinimum(histomin);
511  //hpar1[0]->SetLineColor(1);
512  cout<<"Plotted "<<i<<" aligned detectors"<<endl;
513  }//end loop on pars
514  c_params->cd();
515  TLegend * legend = MakeLegend(.1,.93,.9,.98);
516  legend->Draw();
517  c_params->SaveAs(plotName);
518  std::cout << "Deleting..." << std::flush;
519  delete c_params;
520  std::cout << "Deleting..." << std::flush;
521  //delete f;
522  std::cout << "Deleting..." << std::endl;
523 
524 }//end PlotAlignPars
525 
526 
527 
528 void HIPplots::plotAlignParamsAtIter( int iter, string ShiftsOrParams, char* plotName ){
529 
530  cout<<"Welcome to HIPplots::plotAlignParamsAtIter "<<iter<<endl;
531  bool bParams = false;
532  bool bShifts = false;
533  if (ShiftsOrParams == "PARAMS") bParams = true;
534  else if (ShiftsOrParams == "SHIFTS") bShifts = true;
535  else {cout<<"ERROR in plotAliParamsAtIter!!! Wrong input argument: "<<ShiftsOrParams<<" . Exiting"<<endl; return;}
536 
537  int i = 0;
538  TFile *f = new TFile( _outFile,"READ" );
539  //f->ls();
540 
541  TCanvas *c_params = new TCanvas("can_params","CAN_PARAMS",1200,700);
542  c_params->Divide(3,2);
543  //cout<<"(1) I am in "<<gDirectory->GetPath()<<endl;
544 
545  //TDirectory *d = (TDirectory*)f->Get("ShiftsPar0");
546  //const int ndets = GetNIterations(d,"Par_0_SiDet_");
547 
548  for(int iPar = 0; iPar < 6; iPar++){
549 
550  c_params->cd( iPar+1 );
551 
552  gPad->SetTopMargin(0.15);
553  gPad->SetBottomMargin(0.15);
554  char ppdirname[32];
555  if (bParams) sprintf( ppdirname,"ShiftsPar%d", iPar );
556  if (bShifts) sprintf( ppdirname,"Shifts%d", iPar );
557  if(iPar > 0)gDirectory->cd("../");
558  gDirectory->cd(ppdirname);
559  //cout<<"(2) I am in "<<gDirectory->GetPath()<<endl;
560 
561  TH1D *hiter;
562  char histoname[16];
563  if (bParams) sprintf( histoname, "Par_%d_Iter_%d", iPar, iter );
564  if (bShifts) sprintf( histoname, "AlignedShift_%d_Iter_%d", iPar, iter );
565  hiter = (TH1D*) gDirectory->Get(histoname);
566 
567 
568  //hiter->SetMarkerStyle(1);
569  hiter->SetStats(1);
570 
571  hiter->SetXTitle("Sub-Det");
572  if(iPar==0)hiter->SetYTitle("#delta u (#mum)");
573  else if(iPar==1)hiter->SetYTitle("#delta v (#mum)");
574  else if(iPar==2)hiter->SetYTitle("#delta w (#mum)");
575  else if(iPar==3)hiter->SetYTitle("#delta #alpha (#murad)");
576  else if(iPar==4)hiter->SetYTitle("#delta #beta (#murad)");
577  else if(iPar==5)hiter->SetYTitle("#delta #gamma (#murad)");
578  else hiter->SetXTitle("dunno");
579  hiter->GetYaxis()->SetTitleOffset(1.5);
580 
581  double histomax,histomin;
582  if(iPar==2||iPar==5) {histomax=200; histomin=-200;}
583  else {histomax=100; histomin=-100;}
584 
585  if(iPar>=3)hiter->Scale(1000000.0); //convert from rad to micro rad
586  else hiter->Scale(10000.0); //convert from cm to um
587  double tmpmax=hiter->GetBinContent(hiter->GetMaximumBin());
588  double tmpmin=hiter->GetBinContent(hiter->GetMinimumBin());
589  if(tmpmax>histomax)histomax=tmpmax*1.2;
590  if(tmpmin<histomin)histomin=tmpmin*1.2;
591  hiter->SetMaximum(histomax);
592  hiter->SetMinimum(histomin);
593  //hiter->SetAxisRange(0, 6, "X");
594  hiter->SetLineColor(1);
595  TColor *col = gROOT->GetColor(kCyan+1);
596  col->SetAlpha(0.3);
597  hiter->SetFillColor(col->GetNumber());
598 // hiter->SetFillColor(kCyan);
599  hiter->SetMarkerStyle(7);
600  hiter->Draw();
601  hiter->Draw("PE1");
602 /*
603  if (bShifts) {
604  TH1D *hiter2;
605  char histoname2[16];
606  sprintf( histoname2, "MisalignedShift_%d_Iter_%d", iPar, iter );
607  hiter2 = (TH1D*) gDirectory->Get(histoname2);
608  hiter2->SetLineColor(1);
609  hiter2->SetFillColor(kRed-9);
610  hiter2->Draw("SAME");
611  }
612 
613 */
614 // cout<<"Plotted "<<i<<" aligned detectors"<<endl;
615  }//end loop on pars
616  c_params->SaveAs(plotName);
617  delete c_params;
618  std::cout << "Deleting..." << std::endl;
619  delete f;
620 
621 }//end PlotAlignParamsAtIter
622 
623 
624 
625 void HIPplots::extractAlignableChiSquare( int minHits, int subDet, int doubleSided){
626 
627  //TFile *fp = new TFile( _inFile_params,"READ");
628  // const TList* keysp = fp->GetListOfKeys();
629  // const unsigned int maxIteration = keysp->GetSize();
630  cout<<"\n--- Welcome to extractAlignableChiSquare ---"<<endl;
631  cout<<"\nInput parameters:\n\tMinimum number of hits per alignbale = "<<minHits<<"\n\tSubdetetctor selection"<<subDet<<endl;
632 
633  if(minHits<1){
634  cout<<"Warning ! Allowing to select modules with NO hits. Chi2 not defined for them. Setting automatically minNhits=1 !!!"<<endl;
635  minHits=1;
636  }
637 
638  TFile *fv = new TFile( _inFile_uservars,"READ");
639  const TList* keysv = fv->GetListOfKeys();
640  const unsigned int maxIteration = keysv->GetSize() - 1;
641  cout<<"MaxIteration is "<<maxIteration <<endl;
642 
643  char fileaction[16];
644  if(CheckFileExistence(_outFile))sprintf(fileaction,"UPDATE");
645  else sprintf(fileaction,"NEW");
646  TFile *fout = new TFile( _outFile, fileaction);
647 
648  TTree* tree0 = (TTree*) fv->Get(keysv->At(0)->GetName() );
649  unsigned int detId0;
650  tree0->SetBranchAddress("Id", &detId0 );
651  const int ndets=tree0->GetEntries();
652  TH1D *halichi2n[ndets];//norm chi2 for each module as a function of iteration
653  TH1D *htotchi2n[maxIteration];//distrib of norm chi2 for all modules at a given iteration
654  TH1D *hprobdist[maxIteration];
655  char ppdirname[16];
656  sprintf(ppdirname,"AlignablesChi2n");
657  fout->mkdir(ppdirname);
658  gDirectory->cd(ppdirname);
659 
660  for(int iter = 1; iter <=(int)maxIteration ; iter++){
661  char histoname[64],histotitle[128];
662  sprintf(histoname,"Chi2n_iter%d",iter);
663  sprintf(histotitle,"Distribution of Normalised #chi^{2} for All alignables at iter %d",iter);
664  htotchi2n[iter-1]=new TH1D(histoname,histotitle,1000,0.0,10.0);
665  sprintf(histoname,"ProbChi2_%d",iter);
666  sprintf(histotitle,"Distribution of Prob(#chi^{2},ndof) at iter %d",iter);
667  hprobdist[iter-1]=new TH1D(histoname,histotitle,100,0.0,1.0);
668  }
669  gDirectory->mkdir("AlignablewiseChi2n");
670  gDirectory->cd("AlignablewiseChi2n");
671  //declare histos
672  for(int i = 0; i < (int)ndets; i++){
673  tree0->GetEntry(i);
674  char histoname[64],histotitle[64];
675  sprintf(histoname,"Chi2n_%d",i);
676  sprintf(histotitle,"Normalised #chi^{2} for detector #%d",i);
677  halichi2n[i]=new TH1D(histoname,histotitle,maxIteration,0.5,maxIteration+0.5);
678  }
679 
680  //Loop on trees and fill histos
681  int modules_accepted = 0;
682  for(unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
683 
684  TTree* tmpTreeUV = (TTree*) fv->Get(keysv->At(iter)->GetName() ); //get the UserVariable tree at each iteration
685  cout<<"Taking tree "<<keysv->At(iter)->GetName()<<endl;
686  //tmpTreeUV->GetListOfLeaves()->ls();
687  tmpTreeUV->SetBranchStatus("*",0);
688  tmpTreeUV->SetBranchStatus("Id",1);
689  tmpTreeUV->SetBranchStatus("Nhit",1);
690  tmpTreeUV->SetBranchStatus("AlignableChi2",1);
691  tmpTreeUV->SetBranchStatus("AlignableNdof",1);
692  double alichi2=0.0;
693  unsigned int alindof=0;
694  int nHit=0;
695  unsigned int detId=0;
696  tmpTreeUV->SetBranchAddress("AlignableChi2",&alichi2);
697  tmpTreeUV->SetBranchAddress("AlignableNdof",&alindof);
698  tmpTreeUV->SetBranchAddress("Id",&detId);
699  tmpTreeUV->SetBranchAddress("Nhit",&nHit);
700 
701  modules_accepted=0;
702  //loop on entries
703  for(int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
704  tmpTreeUV->GetEntry(j);
705 
706  bool passSubdetCut = true;
707  if (subDet > 0){ if ( GetSubDet(detId) != subDet ) passSubdetCut = false; } //get subDet Id from module detId,and plot only the selected subDet
708  if (doubleSided > 0){
709  if ( GetSubDet(detId)%2 == 1 && doubleSided == 1 && GetBarrelLayer(detId) < 3 ) passSubdetCut = false;
710  if ( GetSubDet(detId)%2 == 1 && doubleSided == 2 && GetBarrelLayer(detId) > 2 ) passSubdetCut = false;
711  }
712 
713  if (( nHit >= minHits )&&( passSubdetCut ) ){ //cut on nhits per module
714  halichi2n[j]->SetBinContent(iter,double(alichi2 / alindof));
715 // halichi2n[j]->SetBinContent(iter,double(alichi2));
716  double prob=TMath::Prob(double(alichi2),int(alindof));
717 
718  htotchi2n[iter-1]->Fill(double(alichi2 / alindof));
719  hprobdist[iter-1]->Fill(prob);
720  modules_accepted++;
721  }
722 
723  }//end loop on j - alignables
724  cout<<"alignables accepted at iteration "<<iter<<" = "<<modules_accepted<<endl;
725  delete tmpTreeUV;
726  }//end loop on iterations
727 
728  // Ma che e'???
729  /* cout<<"Prob for chi2=0,02 ndof=40 -> "<<TMath::Prob(0.02,40)<<endl;
730  cout<<"Prob for chi2=20, ndof=40 -> "<<TMath::Prob(20.0,40)<<endl;
731  cout<<"Prob for chi2=40, ndof=40 -> "<<TMath::Prob(40.0,40)<<endl;
732  cout<<"Prob for chi2=60, ndof=40 -> "<<TMath::Prob(60.0,40)<<endl;
733  cout<<"Prob for chi2=2000, ndof=40 -> "<<TMath::Prob(2000.0,40)<<endl; */
734 
735  //Save
736  fout->Write();
737  delete fout;
738  delete fv;
739  cout<<"Finished extractAlignableChiSquare"<<endl;
740 }//end extractAlignableChiSquare
741 
742 
743 void HIPplots::extractSurveyResiduals( int currentPar,int subDet ){
744 
745  cout<<"\n--- extractSurveyResiduals has been called ---"<<endl;
746 
747  TFile *fv = new TFile( _inFile_surveys,"READ");
748  const TList* keysv = fv->GetListOfKeys();
749  const unsigned int maxIteration = keysv->GetSize() ;
750  cout<<"MaxIteration is "<<maxIteration <<endl;
751 
752  char fileaction[16];
753  if(CheckFileExistence(_outFile))sprintf(fileaction,"UPDATE");
754  else sprintf(fileaction,"NEW");
755  TFile *fout = new TFile( _outFile, fileaction);
756 
757  TTree* tree0 = (TTree*) fv->Get(keysv->At(0)->GetName() );
758  unsigned int detId0;
759  tree0->SetBranchAddress("Id", &detId0 );
760  const int ndets=tree0->GetEntries();
761  TH1D *hsurvey[ndets];//norm chi2 for each module as a function of iteration
762  TH1D *htotres[maxIteration];//distrib of norm chi2 for all modules at a given iteration
763 
764  char ppdirname[16];
765  sprintf(ppdirname,"SurveyResiduals");
766  fout->mkdir(ppdirname);
767  gDirectory->cd(ppdirname);
768 
769  for(int iter = 1; iter <=(int)maxIteration ; iter++){
770  char histoname[64],histotitle[128];
771  sprintf(histoname,"SurveyRes_Par%d_iter%d",currentPar,iter);
772  sprintf(histotitle,"Distribution of survey Residuals for All alignables ; Par=%d ; iter %d",currentPar,iter);
773  htotres[iter-1]=new TH1D(histoname,histotitle,1000,0.0,10.0);
774  }
775  gDirectory->mkdir("SurveyResiduals_alignables");
776  gDirectory->cd("SurveyResiduals_alignables");
777  //declare histos
778  for(int i = 0; i < (int)ndets; i++){
779  tree0->GetEntry(i);
780  char histoname[64],histotitle[64];
781  sprintf(histoname,"SurveyRes_Par%d_%d",currentPar,i);
782  sprintf(histotitle,"Survey residual for detector #%d - Par=%d",i,currentPar);
783  hsurvey[i]=new TH1D(histoname,histotitle,maxIteration,0.5,maxIteration+0.5);
784  }
785 
786  //Loop on trees and fill histos
787  int modules_accepted = 0;
788  for(unsigned int iter = 1; iter <= maxIteration; iter++){//loop on i (HIP iterations -> trees in file)
789 
790  TTree* tmpTreeUV = (TTree*) fv->Get(keysv->At(iter)->GetName() ); //get the UserVariable tree at each iteration
791  cout<<"Taking tree "<<keysv->At(iter)->GetName()<<endl;
792  //tmpTreeUV->GetListOfLeaves()->ls();
793  tmpTreeUV->SetBranchStatus("*",0);
794  tmpTreeUV->SetBranchStatus("Id",1);
795  tmpTreeUV->SetBranchStatus("Par",1);
796 
797  double par[6];
798  unsigned int detId=0;
799  tmpTreeUV->SetBranchAddress("Par",&par);
800  tmpTreeUV->SetBranchAddress("Id",&detId);
801 
802 
803  modules_accepted=0;
804  //loop on entries
805  for(int j = 0; j < tmpTreeUV->GetEntries(); j++){ //loop on j (modules -> entries in each tree)
806  tmpTreeUV->GetEntry(j);
807 
808  bool passSubdetCut = true;
809  if (subDet > 0){ if ( GetSubDet(detId) != subDet ) passSubdetCut = false; }
810  if ( passSubdetCut ){
811  hsurvey[j]->SetBinContent(iter,double(par[currentPar]));
812 
813  modules_accepted++;
814  }
815 
816  }//end loop on j - alignables
817  cout<<"alignables accepted at iteration "<<iter<<" = "<<modules_accepted<<endl;
818  delete tmpTreeUV;
819  }//end loop on iterations
820 
821  //Save
822  fout->Write();
823  delete fout;
824  delete fv;
825  cout<<"Finished extractAlignableChiSquare"<<endl;
826 }//end extractAlignableChiSquare
827 
828 
829 
830 
831 
832 
833 void HIPplots::plotAlignableChiSquare(char *plotName, float minChi2n){
834 
835  int i = 0;
836  TFile *f = new TFile( _outFile,"READ" );
837  TCanvas *c_alichi2n=new TCanvas("can_alichi2n","CAN_ALIGNABLECHI2N",900,900);
838  c_alichi2n->cd();
839  TDirectory *chi_d1=(TDirectory*)f->Get("AlignablesChi2n");
840  const int maxIteration= GetNIterations(chi_d1,"Chi2n_iter",1)-1;
841  cout<<"N iterations "<<maxIteration<<endl;
842  //take the histos prepared with extractAlignableChiSquare
843  gDirectory->cd("AlignablesChi2n");
844  TDirectory *chi_d=(TDirectory*)gDirectory->Get("AlignablewiseChi2n");
845  const int ndets= GetNIterations(chi_d,"Chi2n_");
846 
847  gDirectory->cd("AlignablewiseChi2n");
848  TH1D *hchi2n[ndets];
849  char histoname[64];
850  int sampling_ratio=1;
851  int ndets_plotted=(int)ndets/sampling_ratio;
852  cout<<"Sampling "<<ndets_plotted<<" detectors over a total of "<<ndets<<endl;
853  double histomax=0.1, histomin=-0.1;
854  bool firstplotted=false;
855  int firstplottedindex=0;
856  int totalplotted=0;
857  bool plothisto=true;
858  while((i<ndets_plotted) && (i*sampling_ratio<ndets) ){
859  sprintf(histoname,"Chi2n_%d",i);
860  hchi2n[i]=(TH1D*)gDirectory->Get(histoname);
861 
862  //if required, check that the chi2n exceeds at any point the threshold set as input
863  bool chi2ncut=false;
864  if(minChi2n>0.0){
865  for(int bin=1;bin<=hchi2n[i]->GetNbinsX();bin++){
866  if(hchi2n[i]->GetBinContent(bin)>minChi2n){chi2ncut=true;break;}
867  }
868  }
869  else chi2ncut=true;
870 
871  //if required check that the chi2 is increasing three iterations in a row and plot only those histos
872  if(plotbadchi2)plothisto=CheckHistoRising(hchi2n[i]);
873  else plothisto=true;
874 
875  if(chi2ncut&&plothisto){
876  double tmpmax=hchi2n[i]->GetBinContent(hchi2n[i]->GetMaximumBin());
877  double tmpmin=hchi2n[i]->GetBinContent(hchi2n[i]->GetMinimumBin());
878  histomax=2.0;
879  histomin=0.0;
880  if(tmpmax>histomax)histomax=tmpmax;
881  if(tmpmin<histomin)histomin=tmpmin;
882 
883  hchi2n[i]->SetMaximum(histomax);
884  hchi2n[i]->SetMinimum(histomin);
885  hchi2n[i]->SetStats(0);
886 
887 
888  if(!firstplotted){
889  hchi2n[i]->SetXTitle("Iteration");
890  hchi2n[i]->SetYTitle("#chi^{2} / # dof");
891 // hchi2n[i]->SetYTitle("#chi^{2}");
892  hchi2n[i]->SetTitle("Reduced #chi^{2} for alignables");
893  hchi2n[i]->Draw("PL");
894  firstplotted=true;
895  firstplottedindex=i;
896  }
897  else hchi2n[i]->Draw("PLsame");
898  totalplotted++;
899  }//END IF plot it because it passed the chi2n cut
900  i++;
901  }//end while loop on alignables
902  hchi2n[firstplottedindex]->SetMaximum(histomax*1.2);
903  hchi2n[firstplottedindex]->SetMinimum(histomin*0.8);
904  //override auto-resize of axis scale
905  hchi2n[firstplottedindex]->SetMaximum(histomax);
906  hchi2n[firstplottedindex]->SetMinimum(histomin);
907 
908  cout<<"Plotted "<<totalplotted<<" alignables over an initial sample of "<<ndets_plotted<<endl;
909  TText *txtchi2n_1=new TText();
910  txtchi2n_1->SetTextFont(63);
911  txtchi2n_1->SetTextSize(22);
912  char strchi2n_1[128];
913  sprintf(strchi2n_1,"Plotted alignables (Chi2n > %.3f): %d / %d",minChi2n,totalplotted,ndets_plotted);
914  txtchi2n_1->DrawText(1.2,0.0,strchi2n_1);
915  char finplotname[192];
916  sprintf(finplotname,"%s.png",plotName);
917  c_alichi2n->SaveAs(finplotname);
918  delete c_alichi2n;
919  //delete hchi2n;
920 
921 
922  cout<<"Doing distrib"<<endl;
923  gDirectory->cd("../");
924  TCanvas *c_chi2ndist=new TCanvas("can_chi2ndistr","CAN_CHI2N_DISTRIBUTION",900,900);
925  c_chi2ndist->cd();
926  TH1D *hiter[maxIteration];
927  TLegend *l=new TLegend(0.7,0.7,0.9,0.9);
928  l->SetFillColor(0);
929  l->SetBorderSize(0);
930  int colors[10]={1,2,8,4,6,7,94,52,41,45};
931  int taken=0;
932  float tmpmax=1.0;
933  float newmax=0.0;
934  for(i=0;i<maxIteration;i++){
935  sprintf(histoname,"Chi2n_iter%d",i+1);
936  hiter[i]=(TH1D*)gDirectory->Get(histoname);
937  hiter[i]->SetXTitle("#chi^{2} / dof");
938  hiter[i]->SetYTitle("Alignables");
939  hiter[i]->SetTitle("Normalised #chi^{2} of Alignables");
940  hiter[i]->Rebin(5);
941  hiter[i]->GetXaxis()->SetRangeUser(0.0,3.0);
942  hiter[i]->SetLineColor(i);
943 
944  char legend_entry[64];
945  float histmean=hiter[i]->GetMean();
946  if(i==0){
947  hiter[i]->SetLineColor(colors[taken]);
948  taken++;
949  sprintf(legend_entry,"Norm #chi^{2}; Iter %d; %.4f",i,histmean);
950  l->AddEntry(hiter[i],legend_entry,"L");
951  tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
952  newmax=tmpmax;
953  // hiter[i]->Draw("");
954  }
955  else if((i+1)%5==0){
956  hiter[i]->SetLineColor(colors[taken]);
957  taken++;
958  sprintf(legend_entry,"Norm #chi^{2}; Iter %d; %.4f",i+1,histmean);
959  l->AddEntry(hiter[i],legend_entry,"L");
960  tmpmax=hiter[i]->GetBinContent(hiter[i]->GetMaximumBin());
961  if(tmpmax>newmax)newmax=tmpmax;
962  // hiter[i]->Draw("same");
963  }
964  }
965  cout<<"NewMax after 1st loop -> "<<newmax<<endl;
966 
967  for(i=0;i<maxIteration;i++){
968  hiter[i]->SetMaximum(newmax*1.1);
969  if(i==1) hiter[i]->Draw("");
970  else if((i+1)%5==0) hiter[i]->Draw("same");
971  }
972  l->Draw();
973 
974  sprintf(finplotname,"%s_distr.png",plotName);
975  cout<<finplotname<<endl;
976  c_chi2ndist->SaveAs(finplotname);
977  c_chi2ndist->SetLogy();
978  sprintf(finplotname,"%s_distrlog.png",plotName);
979  cout<<finplotname<<endl;
980  c_chi2ndist->SaveAs(finplotname);
981 
982  delete c_chi2ndist;
983  delete f;
984 }//end plotAlignableChiSquare
985 
986 
987 
988 // -----------------------------------------------------------------
989 // private classes
990 int HIPplots::GetNIterations(TDirectory *f,char *tag,int startingcounter){
991  int fin_iter=0,i=startingcounter;
992  bool obj_exist=kTRUE;
993  while(obj_exist){
994  char objname[32];
995  sprintf(objname,"%s%d",tag,i);
996  if(!f->FindObjectAny(objname))obj_exist=kFALSE;
997  fin_iter=i;
998  i++;
999  }
1000  cout<<"Max Iterations is "<<fin_iter<<endl;
1001 
1002  return fin_iter;
1003 }
1004 
1005 int HIPplots::GetSubDet( unsigned int id ){
1006 
1007  const int reserved_subdetectorstartbit=25;
1008  const int reserved_subdetectorfinalbit=27;
1009 
1010  unsigned int detID = id;
1011 
1012  int shift = 31-reserved_subdetectorfinalbit;
1013  detID = detID<<(shift);
1014  shift = reserved_subdetectorstartbit + shift;
1015  detID = detID>>(shift);
1016 
1017  return detID;
1018 
1019 }
1020 
1021 int HIPplots::GetBarrelLayer( unsigned int id ){
1022 
1023  const int reserved_layerstartbit=14;
1024  const int reserved_layermask=0x7;
1025 
1026  return int((id>>reserved_layerstartbit) & reserved_layermask);
1027 
1028 }
1029 
1030 void HIPplots::SetMinMax( TH1* h ){
1031 
1032  Double_t xmin = 10000;
1033  Double_t xmax = -10000;
1034 
1035 
1036  for(int i = 1; i <= h->GetNbinsX(); ++i) {
1037  if((h->GetBinContent(i) > 0)&&(h->GetBinCenter(i)>xmax) ) xmax = h->GetBinCenter(i);
1038  }
1039  for(int i = h->GetNbinsX(); i >= 1; --i) {
1040  if((h->GetBinContent(i) > 0)&&(h->GetBinCenter(i)<xmin) ) xmin = h->GetBinCenter(i);
1041  }
1042 
1043  h->SetAxisRange((xmin-xmin*0.1), (xmax+xmax*0.1), "X");
1044  // std::cout << "xmin: " << xmin << ", xmax: " << xmax << std::endl;
1045 
1046 }
1047 
1048 
1049 
1050 void HIPplots::plotHitMap( char *outpath,int subDet,int minHits){
1051  cout<<"Starting plotHitMap"<<flush;
1052 
1053  TFile *falignable=new TFile(_inFile_HIPalign,"READ");
1054  cout<<"\tLoaded file"<<flush;
1055  //take the alignablewise tree and address useful branches
1056  TTree *talignable=(TTree*)falignable->Get("T2");
1057  cout<<"\t Loaded tree"<<endl;
1058 
1059  float eta=-999.0,phi=-55.0, xpos=-999.0,ypos=+999.0,zpos=-11111.0;
1060  int layer=-1,type=-1, nhit=-11111;
1061  int id=0;
1062  talignable->SetBranchAddress("Id",&id);
1063  talignable->SetBranchAddress("Layer",&layer);
1064  talignable->SetBranchAddress("Type",&type);
1065  talignable->SetBranchAddress("Nhit",&nhit);
1066  talignable->SetBranchAddress("Ypos",&ypos);
1067  talignable->SetBranchAddress("Eta",&eta);
1068  talignable->SetBranchAddress("Phi",&phi);
1069  talignable->SetBranchAddress("Ypos",&ypos);
1070  talignable->SetBranchAddress("Xpos",&xpos);
1071  talignable->SetBranchAddress("Zpos",&zpos);
1072 
1073  //loop on subdets
1074  char typetag[16];
1075  int maxLayers=0;
1076  if(subDet == TPBid){sprintf(typetag,"TPB");maxLayers=3;}
1077  else if(subDet == TPEid){sprintf(typetag,"TPE");maxLayers=2;}
1078  else if(subDet == TIBid){sprintf(typetag,"TIB");maxLayers=4;}
1079  else if(subDet == TIDid){sprintf(typetag,"TID");maxLayers=3;}
1080  else if(subDet == TOBid){sprintf(typetag,"TOB");maxLayers=6;}
1081  else if(subDet == TECid){sprintf(typetag,"TEC");maxLayers=9;}
1082  else {sprintf(typetag,"UNKNOWN");}
1083  cout<<"Starting to plot Hit Distributions for "<<typetag<<endl;
1084 
1085  bool printbinning=true;
1086  char psname[600];
1087  char ovtag[16];
1088  sprintf(psname,"%s/Hits_%s_Layers_Skimmed.eps",outpath,typetag);
1089 
1090  char binfilename[64];
1091  sprintf(binfilename,"./BinningHitMaps_%s.txt",typetag);
1092  ofstream binfile(binfilename,ios::out);
1093 
1094  if(printbinning){
1095  binfile<<"******** Binning for Subdet "<<typetag<<" *********"<<endl<<endl;
1096  }
1097 
1098  for(int layerindex=1;layerindex<=maxLayers;layerindex++){ //loop on layers
1099 
1100  cout<<"\n\n*** Layer # "<<layerindex<<" ***"<<endl;
1101  //activate only useful branches
1102  talignable->SetBranchStatus("*",0);
1103  talignable->SetBranchStatus("Id",1);
1104  talignable->SetBranchStatus("Type",1);
1105  talignable->SetBranchStatus("Layer",1);
1106  talignable->SetBranchStatus("Nhit",1);
1107  talignable->SetBranchStatus("Eta",1);
1108  talignable->SetBranchStatus("Phi",1);
1109  talignable->SetBranchStatus("Ypos",1);
1110  talignable->SetBranchStatus("Zpos",1);
1111  talignable->SetBranchStatus("Xpos",1);
1112 
1113  TCut selA,selECpos,selECneg;
1114  char selA_str[196],selECneg_str[196],selECpos_str[196];
1115  char varA_str[64],varB_str[64];
1116  char commonsense_str[128]="Xpos>-150.0&&Xpos<150.0&&Ypos>-150.0&&Ypos<150.0&&Zpos>-400.0&&Zpos<400.0";
1117  TCut commonsense=commonsense_str;
1118 
1119  sprintf(selECneg_str,"(Type==%d && Layer==%d && Zpos<0.0 && Nhit>=%d && sqrt(pow(Xpos,2)+pow(Ypos,2) )<125.0 )",subDet,layerindex,minHits);
1120  sprintf(selECpos_str,"(Type==%d && Layer==%d && Zpos>=0.0 && Nhit>=%d && sqrt(pow(Xpos,2)+pow(Ypos,2) )<125.0 )",subDet,layerindex,minHits);
1121  // sprintf(selB_str,"(Eta>0.0)&&(Eta<0.2)");
1122  sprintf(selA_str,"Type==%d && Layer==%d",subDet,layerindex);
1123 
1124  sprintf(varA_str,"Eta>>hvarx");
1125  sprintf(varB_str,"Phi>>hvary");
1126 
1127  selECneg=selECneg_str;
1128  selECpos=selECpos_str;
1129  selA=selA_str;
1130  cout<<"Cuts defined as "<<selA<<endl;
1134 
1135  //--------- (2) START bin definition -----
1136  //cout<<"Selection is "<<sel<<endl;
1137  // char sel[96];
1138 
1139  int nzentries= talignable->Draw("Zpos>>hZ(360,-270.0,270.0)",commonsense&&selA,"goff");
1140  TH1F *hZ=(TH1F*)gDirectory->Get("hZ");
1141  if(subDet == TOBid) SetPeakThreshold(8.0);
1142  else SetPeakThreshold(5.0);
1143  float Zpeaks[120];
1144  int bin=0;
1145  for(bin=0;bin<120;bin++){
1146  Zpeaks[bin]=-444.0;
1147  }
1148  const int nZpeaks=FindPeaks(hZ,Zpeaks,99);
1149  const int nZbinlims=nZpeaks+1;
1150  float Zwidth=(Zpeaks[nZpeaks-1]-Zpeaks[0])/ (nZpeaks-1) ;
1151  float Zmin=Zpeaks[0]- Zwidth/2.0 ;
1152  float Zmax=Zpeaks[nZpeaks-1] + Zwidth/2.0;//((Zpeaks[nZbinlims-1]-Zpeaks[nZbinlims-2])/2.0) ;
1153  cout<<"--> Zmin= "<<Zmin<<" - Zmax= "<<Zmax<<" Zwidth= "<<Zwidth<<" ; found "<<nZpeaks<<" Zpeaks"<<endl;
1154  cout<<"Zpeaks[0] is "<<Zpeaks[0]<<endl;
1155 
1156 
1157  float Phipeaks[120];
1158  if((subDet==TIBid||subDet==TOBid)&&layerindex<=2) sprintf(selA_str,"%s&&Zpos>%f&&Zpos<%f",selA_str,Zpeaks[0]-2.0,Zpeaks[0]+2.0);//DS modules
1159  else sprintf(selA_str,"%s&&Zpos>%f&&Zpos<%f",selA_str,Zpeaks[0]-2.0,Zpeaks[0]+2.0);
1160  int nphientries=talignable->Draw("Phi>>hPhi",selA_str,"goff");
1161  cout<<"N phi entries "<<nphientries<<" from sel "<<selA_str<< endl;
1162  TH1F *hPhi=(TH1F*)gDirectory->Get("hPhi");
1163  // nPhibins=FindPeaks(hPhi,Phipeaks,99);
1164  if(subDet==TPBid&&layerindex==1)nphientries=nphientries-1;
1165  const int nPhibins=nphientries;
1166  cout<<"+ It would have found "<<nPhibins<<" phi bins"<<endl;
1167 
1168  //fill Z array with binning. special case for DS layers of TIB with non equidistant z bins
1169  float phibin[nPhibins+1];
1170  float zbin[nZbinlims];
1171  float Phiwidth=6.28/nPhibins;
1172 
1173  if((subDet==TIBid||subDet==TOBid)&&layerindex<=2){//DS modules in TIB and TOB
1174 
1176  cout<<"Gonna loop over "<<nZpeaks<<" peaks / "<<nZbinlims<<" bin limits"<<endl;
1177  for(bin=0;bin<nZbinlims-1;bin++){
1178  float zup=0.0;
1179  float zdown=0.0;
1180 
1181  if(bin==0){ //don't go underflow!
1182  zup=(Zpeaks[bin+1]-Zpeaks[bin])/2.0;
1183  zdown=zup;
1184  }
1185  else if(bin==nZbinlims-2){//don't go overflow!
1186  cout<<"Don't go overflow !"<<endl;
1187  zdown=(Zpeaks[bin]-Zpeaks[bin-1])/2.0;
1188  if(layerindex==1) zup=(Zpeaks[bin-1]-Zpeaks[bin-2])/2.0;
1189  else zup=zdown;
1190  }
1191  else{
1192  zup=(Zpeaks[bin+1]-Zpeaks[bin])/2.0;
1193  zdown=(Zpeaks[bin]-Zpeaks[bin-1])/2.0;
1194  }
1195  zbin[bin] = Zpeaks[bin]-zdown;
1196  zbin[bin+1]= Zpeaks[bin]+zup;
1197 
1198  }//end loop on z bin
1200 
1201 
1202  for(int bin=0;bin<=nPhibins;++bin){
1203  if(bin==0)phibin[bin]=-3.14+bin*(Phiwidth)+Phiwidth/4;
1204  else if(bin==nPhibins-1)phibin[bin]=-3.14+bin*(Phiwidth)-Phiwidth/4;
1205  // else phibin[bin]=-3.14+bin*(Phiwidth);
1206  else phibin[bin]=phibin[bin-1]+ Phiwidth;
1207  }//end loop on phi bin
1208 
1209  }//end if DS layers of TIB/TOB
1210  else{
1211  for(int bin=0;bin<nZbinlims;++bin){
1212  zbin[bin]=Zmin+bin*(Zwidth);
1213  }//end loop on z bin
1214 
1215  for(int bin=0;bin<=nPhibins;++bin){
1216  phibin[bin]=-3.14+bin*(Phiwidth);
1217  }//end loop on phi bin
1218  }
1219 
1220 
1221  float Phimin=Phipeaks[0]- ((Phipeaks[1]-Phipeaks[0])/2.0) ;
1222  float Phimax=Phipeaks[nPhibins-1] + ((Phipeaks[nPhibins-1]-Phipeaks[nPhibins-2])/2.0) ;
1223 
1224  cout<<"N Z bin LIMITs = "<<nZbinlims<<" N Phi bin LIMITS = "<<nPhibins<<endl;
1225  //--------- END bin definition -----
1226 
1227 
1228 
1229  char histoname[64];
1230  sprintf(histoname,"%s_Layer%d",typetag,layerindex);
1231  TH2F *hetaphi=new TH2F(histoname,histoname,nZpeaks,zbin,nPhibins,phibin);
1232 
1233 
1234  cout<<"Starting to loop on entries"<<flush;
1235  int nmods=0;
1236  int nlowentrycells=0;
1237 
1238  for(int j=0;j<talignable->GetEntries();j++){ //loop on entries (-> alignables)
1239  if(j%1000==0)cout<<"."<<flush;
1240  talignable->GetEntry(j);
1241  if(type==subDet&&layer==layerindex){
1242  if(nhit>=minHits){
1243  //find the right eta bin
1244  hetaphi->Fill(zpos,phi,nhit);
1245  nmods++;
1246  }
1247  else{
1248  hetaphi->Fill(zpos,phi,-99);
1249  nlowentrycells++;
1250  }//end else if nhits < minhits
1251 
1252  }//end if the type and layer are the desired ones
1253  }//end loop on entries(->alignables)
1254 
1255  //Let's start with the funny things: fancy graphics!
1256  hetaphi->SetXTitle("Z [cm]");
1257  hetaphi->SetYTitle("#phi (rad)");
1258 
1259  int Nxbins=hetaphi->GetXaxis()->GetNbins();
1260  int Nybins=hetaphi->GetYaxis()->GetNbins();
1261  cout<<"On x-axis there are "<<Nxbins<<" bins "<<endl;
1262  cout<<"On y-axis there are "<<Nybins<<" bins "<<endl;
1263 
1264 
1265  bool smooth_etaphi=false;
1266  if(smooth_etaphi){
1267  for(int i=1;i<=Nxbins;i++){
1268  for(int j=1;j<=Nybins;j++){
1269  float bincont=hetaphi->GetBinContent(i,j);
1270  if(bincont==0){//average with the surrounding bins
1271  float binup1=hetaphi->GetBinContent(i,j+1);
1272  float bindown1=hetaphi->GetBinContent(i,j-1);
1273  float binlx1=hetaphi->GetBinContent(i-1,j);
1274  float binrx1=hetaphi->GetBinContent(i+1,j);
1275  if(i==1)binlx1=binrx1;//at the edges you don't have lx or rx bins; set a def value
1276  if(i==Nxbins)binrx1=binlx1;
1277  if(j==1)bindown1=binup1;
1278  if(j==Nybins)binup1=bindown1;
1279  int adiacentbins=0;
1280  if(binup1>0.0)adiacentbins++;
1281  if(bindown1>0.0)adiacentbins++;
1282  if(binlx1>0.0)adiacentbins++;
1283  if(binrx1>0.0)adiacentbins++;
1284  if(adiacentbins>=2){
1285  float avg=(binup1+bindown1+binlx1+binrx1)/adiacentbins;
1286  hetaphi->SetBinContent(i,j,avg);
1287  }
1288  }
1289  }
1290  }
1291  }//end if smooth_etaphi
1292 
1293 
1294  //for debugging purposes
1295  TGraph *gretaphi;
1296  bool plotAlignablePos=false;
1297  if(plotAlignablePos){
1298  const int ngrpoints=nmods;
1299  float etagr[ngrpoints],phigr[ngrpoints];
1300  nmods=0;
1301  for(int j=0;j<talignable->GetEntries();j++){ //loop on entries (-> alignables)
1302  if(j%1000==0)cout<<"."<<flush;
1303  talignable->GetEntry(j);
1304  if(type==subDet&&layer==layerindex){
1305  etagr[nmods]=zpos;
1306  phigr[nmods]=phi;
1307  nmods++;
1308  }
1309  }
1310 
1311  gretaphi=new TGraph(ngrpoints,etagr,phigr);
1312  gretaphi->SetMarkerStyle(20);
1313  gretaphi->SetMarkerColor(1);
1314  gretaphi->SetMarkerSize(0.75);
1315  }
1316 
1317 
1318 
1320  // cout<<"Layer #"<<i<<endl;
1321  float Zcellgr[512];
1322  float Phicellgr[512];
1323  int nemptycells=0;
1324  for(int zcells=1;zcells<=hetaphi->GetNbinsX();zcells++){
1325  for(int phicells=1;phicells<=hetaphi->GetNbinsY();phicells++){
1326  if( hetaphi->GetBinContent( zcells,phicells)==-99){
1327  hetaphi->SetBinContent( zcells,phicells,0);
1328  Zcellgr[nemptycells]= float(hetaphi->GetXaxis()->GetBinCenter(zcells));
1329  Phicellgr[nemptycells]= float(hetaphi->GetYaxis()->GetBinCenter(phicells));
1330  nemptycells++;
1331  }
1332  // if(a==2)cout<<"Finished Z value "<<hetaphi->GetXaxis()->GetBinCenter(zcells)<<" the emptycells are "<<nemptycells<<endl;
1333  }//end loop on Phi bins
1334  }//end loop on Z bins
1335  TGraph *gr_empty=new TGraph( nlowentrycells,Zcellgr,Phicellgr);
1336  sprintf(histoname,"gr_emptycells_subdet%d_layer%d",subDet,layerindex);
1337  // cout<<"Graph name: "<<histoname<<" with "<<gr_empty->GetN()<<endl;
1338  gr_empty->SetName(histoname);
1339  gr_empty->SetMarkerStyle(5);
1341 
1342 
1343  cout<<" Done! Used "<<nmods<<" alignables. Starting to plot !"<<endl;
1344 
1345  //plot them and save the canvas appending it to a ps file
1346  gStyle->SetPalette(1,0);
1347  TCanvas *c_barrels=new TCanvas("canvas_hits_barrels","CANVAS_HITS_BARRELS",1600,1600);
1348  TCanvas *c_endcaps=new TCanvas("canvas_hits_endcaps","CANVAS_HITS_ENDCAPS",3200,1600);
1349  TPad *pleft=new TPad("left_panel","Left Panel",0.0,0.0,0.49,0.99);
1350  TPad *pcent=new TPad("central_up_panel","Central Panel",0.01,0.00,0.99,0.99);
1351  TPad *pright=new TPad("right_panel","Right Panel",0.51,0.0,0.99,0.99);
1352 
1353 
1354  if(subDet==1 ||subDet==3 ||subDet==5 ){
1355  c_barrels->cd();
1356  pcent->Draw();
1357  pcent->cd();
1358  gPad->SetRightMargin(0.15);
1359 
1360  // hbarrel->SetMarkerSize(2.0);
1361  // hbarrel->SetMarkerSize(21);
1362  hetaphi->SetStats(0);
1363  if(subDet==TPBid||subDet==TIBid||subDet==TOBid)hetaphi->Draw("COLZtext");//COLZtext
1364  if(plotAlignablePos) gretaphi->Draw("P");
1365  gr_empty->Draw("P");
1366 
1367  if(printbinning){
1368  binfile<<"--> Layer #"<<layerindex<<endl;
1369  binfile<<"Eta Binning: "<<flush;
1370  for(int h=1;h<=hetaphi->GetNbinsX();h++){
1371  binfile<<hetaphi->GetXaxis()->GetBinLowEdge(h)<<"\t";
1372  if(h==hetaphi->GetNbinsX()) binfile<<hetaphi->GetXaxis()->GetBinLowEdge(h)+hetaphi->GetXaxis()->GetBinWidth(h);
1373  }
1374  binfile<<endl;
1375  binfile<<"Phi Binning: "<<flush;
1376  for(int h=1;h<=hetaphi->GetNbinsY();h++){
1377  binfile<<hetaphi->GetYaxis()->GetBinLowEdge(h)<<"\t";
1378  if(h==hetaphi->GetNbinsX()) binfile<<hetaphi->GetYaxis()->GetBinLowEdge(h)+hetaphi->GetYaxis()->GetBinWidth(h);
1379  }
1380  binfile<<endl;
1381  }
1382 
1383  }
1384  else{
1385  c_endcaps->cd();
1386  pleft->Draw();
1387  pright->Draw();
1388 
1389 
1390  pleft->cd();
1391  gPad->SetRightMargin(0.15);
1392  char selEC_str[192], varEC_str[128];
1393  float radlimit=0.0;
1394  if(subDet == TPBid){radlimit=45.0;}
1395  else if(subDet == TPEid){radlimit=45.0;}
1396  else if(subDet == TIBid){radlimit=70.0;}
1397  else if(subDet == TIDid){radlimit=70.0;}
1398  else if(subDet == TOBid){radlimit=130.0;}
1399  else if(subDet == TECid){radlimit=130.0;}
1400  else {radlimit=0.0;}
1401 
1402  sprintf(varEC_str,"Ypos:Xpos>>hxy_negz(30,%f,%f)",-radlimit,radlimit);
1403  sprintf(selEC_str,"Nhit*(%s&&%s)",selECneg_str,commonsense_str);
1404  cout<<selEC_str<<endl;
1405  int selentriesECneg=talignable->Draw(varEC_str,selEC_str,"goff"); //total number of alignables for thys type and layer
1406  if(selentriesECneg>0){
1407  TH2F *hxy_negz=(TH2F*)gDirectory->Get("hxy_negz");
1408  hxy_negz->GetXaxis()->SetRangeUser(-radlimit,radlimit);
1409  hxy_negz->GetYaxis()->SetRangeUser(-radlimit,radlimit);
1410  char histoname_negz[32];
1411  sprintf(histoname_negz,"%s (-Z)",histoname);
1412  hxy_negz->SetStats(0);
1413  hxy_negz->SetTitle( histoname_negz);
1414  hxy_negz->SetXTitle("X (cm)");
1415  hxy_negz->SetYTitle("Y (cm)");
1416  hxy_negz->Draw("COLZ");
1417  }
1418  else{
1419  cout<<"WARNING !!!! No hits on this layer ! not plotting (-Z) !"<<endl;
1420  }
1421 
1422 
1423 
1424  cout<<"PAD 3"<<endl;
1425  pright->cd();
1426  gPad->SetRightMargin(0.15);
1427  sprintf(selEC_str,"Nhit*(%s&&%s)",selECpos_str,commonsense_str);
1428  // selEC=selEC_str&&commonsense;
1429  cout<<"(2)"<<selEC_str<<endl;
1430  sprintf(varEC_str,"Ypos:Xpos>>hxy_posz(30,%f,%f)",-radlimit,radlimit);
1431  int selentriesECpos=talignable->Draw(varEC_str,selEC_str,"goff"); //total number of alignables for thys type and layer
1432  if(selentriesECpos>0){
1433  TH2F *hxy_posz=(TH2F*)gDirectory->Get("hxy_posz");
1434  char histoname_posz[32];
1435  hxy_posz->GetXaxis()->SetRangeUser(-radlimit,radlimit);
1436  hxy_posz->GetYaxis()->SetRangeUser(-radlimit,radlimit);
1437  sprintf(histoname_posz,"%s (+Z)",histoname);
1438  hxy_posz->SetStats(0);
1439  hxy_posz->SetTitle( histoname_posz);
1440  hxy_posz->SetXTitle("X (cm)");
1441  hxy_posz->SetYTitle("Y (cm)");
1442  hxy_posz->Draw("COLZ");
1443  }
1444  else{
1445  cout<<"WARNING !!!! No hits on this layer ! not plotting (+Z) !"<<endl;
1446  }
1447 
1448 
1449  }
1450 
1451  //save in a ps file
1452  cout << "******* " << typetag << endl;
1453  char psnamefinal[600];
1454 
1455  if(layerindex==1) sprintf(psnamefinal,"%s(",psname);
1456  else if(layerindex==maxLayers) sprintf(psnamefinal,"%s)",psname);
1457  else sprintf(psnamefinal,"%s",psname);
1458 
1459  cout<<"Saving in "<<psnamefinal<<endl;
1460  if(subDet==1 ||subDet==3 ||subDet==5 )c_barrels->SaveAs(psnamefinal);
1461  else c_endcaps->SaveAs(psnamefinal);
1462 
1463  if(subDet==1 ||subDet==3 ||subDet==5 )delete c_barrels;
1464  else delete c_endcaps;
1465 
1466 
1467  // ctmp->cd();
1468  // hbarrel->Draw("scat");
1469 
1470  }//end loop on layers
1471  cout<<"Finished "<<maxLayers<<" of the "<<typetag<<endl;
1472 
1473  delete talignable;
1474  delete falignable;
1475 }//end PlotHitDistribution
1476 
1477 
1478 void HIPplots::dumpAlignedModules(int nhits){
1479 
1481 
1482 }//end dumpAlignedModules
1483 
1484 
1485 
1486 int HIPplots::FindPeaks(TH1F *h1,float *peaklist,const int maxNpeaks,int startbin,int endbin){
1487 
1488  int Npeaks=0;
1489  if(startbin<0)startbin=1;
1490  if(endbin<0)endbin=h1->GetNbinsX();
1491  int bin=startbin;
1492  int prevbin=startbin;
1493  bool rising=true;
1494  while(bin<=endbin){
1495 
1496  float bincont=h1->GetBinContent(bin);
1497  float prevbincont=h1->GetBinContent(prevbin);
1498 
1499  if(prevbincont>peakthreshold||bincont>1.0){
1500  if(bincont>=prevbincont){//we are rising, keep going until we don't go down
1501  rising=true;
1502  }
1503  else{//ehi! we are decreasing. But check if we were decreasing already at the step before
1504  if(!rising){//we were already decreasing, this is not a maximum
1505  rising=true;
1506  }
1507  else{//we found a maximum
1508  rising=false;
1509  peaklist[Npeaks]=h1->GetBinCenter(prevbin);
1510  Npeaks++;
1511  }
1512  }
1513  }//end if bincont>1.0
1514  if(Npeaks>=maxNpeaks){//max number reached
1515  bin=endbin;
1516  }
1517  bin++;
1518  prevbin=bin-1;
1519  }//end while loop on bins
1520 
1521 
1522  return Npeaks;
1523 }//End FindPeaks
1524 
1525 
1526 void HIPplots::SetPeakThreshold(float newpeakthreshold){
1527  peakthreshold=newpeakthreshold;
1528 }
1529 
1530 bool HIPplots::CheckFileExistence(char *filename){
1531  bool flag = false;
1532  fstream fin;
1533  fin.open(filename,ios::in);
1534  if( fin.is_open() ) flag=true;
1535  fin.close();
1536  return flag;
1537 }
1538 
1539 void HIPplots::CheckFiles(int &ierr){
1540 
1541  if (! CheckFileExistence(_inFile_uservars)){
1542  cout<<"Missing file "<<_inFile_uservars<<endl;
1543  ierr++;
1544  }
1545 
1546  if(! CheckFileExistence(_inFile_HIPalign)){
1547  cout<<"Missing file "<<_inFile_HIPalign<<endl;
1548  ierr++;
1549  }
1550 
1551  if(! CheckFileExistence( _inFile_alipos)){
1552  cout<<"Missing file "<< _inFile_alipos<<endl;
1553  ierr++;
1554  }
1555 
1556  if( CheckFileExistence( _outFile)){
1557  cout<<"Output file already existing !"<<endl;
1558  ierr=-1*(ierr+1);
1559  }
1560 
1561 }
1562 
1563 
1564 bool HIPplots::CheckHistoRising(TH1D *h){
1565 
1566  bool rise1=false,rise2=false,rise3=false;
1567  int totbins=h->GetNbinsX();
1568  // for(int i=1;i<=totbins;i++){
1569  int i=1;
1570 
1571 
1572  for(i=1;i<=totbins-3;i++){
1573  double cont1 = h->GetBinContent(i);
1574  if( h->GetBinContent(i+1)>cont1)rise1=true;
1575  if(rise1){
1576  if( h->GetBinContent(i+2)>h->GetBinContent(i+1))rise2=true;
1577  if(rise2){
1578  if( h->GetBinContent(i+3)>h->GetBinContent(i+2)){
1579  rise3=true;
1580  break;
1581  }
1582  }
1583  }
1584 
1585  }//emnd loop on bins
1586 
1587  return rise3;
1588 
1589 }//end CheckHistoRising
1590 
type
type
Definition: HCALResponse.h:21
i
int i
Definition: DBlmapReader.cc:9
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
Gflash::Zmax
const double Zmax[kNumberCalorimeter]
Definition: GflashNameSpace.h:28
HIPplots::plotAlignParams
void plotAlignParams(string ShiftsOrParams, char *plotName="test.png")
Definition: HIPplots.cc:385
HIPplots::MakeLegend
TLegend * MakeLegend(double x1=0.1, double y1=0.1, double x2=0.1, double y2=0.1)
Definition: HIPplots.cc:38
HIPplots::CheckFileExistence
bool CheckFileExistence(char *filename)
Definition: HIPplots.cc:1530
HIPplots::GetSubDet
int GetSubDet(unsigned int id)
Definition: HIPplots.cc:1005
OccupancyClient_cfi.minHits
int minHits
Definition: OccupancyClient_cfi.py:5
HIPplots::plotAlignableChiSquare
void plotAlignableChiSquare(char *plotName="testchi2.png", float minChi2n=-1.0)
Definition: HIPplots.cc:833
HIPplots::SetMinMax
void SetMinMax(TH1 *h)
Definition: HIPplots.cc:1030
HIPplots::_path
char _path[256]
Definition: HIPplots.h:29
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
HIPplots::GetBarrelLayer
int GetBarrelLayer(unsigned int id)
Definition: HIPplots.cc:1021
HIPplots::CheckFiles
void CheckFiles(int &ierr)
Definition: HIPplots.cc:1539
ztail.d
tuple d
Definition: ztail.py:151
colors
vector< Color_t > colors
Definition: trackSplitPlot.h:37
GlobalPosition_Frontier_DevDB_cff.tag
tuple tag
Definition: GlobalPosition_Frontier_DevDB_cff.py:11
HIPplots::dumpAlignedModules
void dumpAlignedModules(int nhits=0)
Definition: HIPplots.cc:1478
GetRecoTauVFromDQM_MC_cff.outFile
tuple outFile
Definition: GetRecoTauVFromDQM_MC_cff.py:20
HIPplots::extractSurveyResiduals
void extractSurveyResiduals(int currentPar, int subDet=0)
Definition: HIPplots.cc:743
HIPplots::HIPplots
HIPplots(int IOV, char *path, char *outFile)
Definition: HIPplots.cc:20
j
int j
Definition: DBlmapReader.cc:9
HIPplots::_inFile_mispos
char _inFile_mispos[256]
Definition: HIPplots.h:35
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
HIPplots::CheckHistoRising
bool CheckHistoRising(TH1D *h)
Definition: HIPplots.cc:1564
HIPplots::extractAlignShifts
void extractAlignShifts(int i, int minHits=0, int subDet=0)
Definition: HIPplots.cc:190
HIPplots::plotbadchi2
bool plotbadchi2
Definition: HIPplots.h:50
patRefSel_outputModule_cff.outpath
tuple outpath
Definition: patRefSel_outputModule_cff.py:5
HIPplots::_inFile_surveys
char _inFile_surveys[256]
Definition: HIPplots.h:37
HIPplots::GetNIterations
int GetNIterations(TDirectory *f, char *tag, int startingcounter=0)
Definition: HIPplots.cc:990
HIPplots::_inFile_truepos
char _inFile_truepos[256]
Definition: HIPplots.h:33
HIPplots::extractAlignableChiSquare
void extractAlignableChiSquare(int minHits=0, int subDet=0, int doubleSided=0)
Definition: HIPplots.cc:625
HIPplots::_inFile_alipos
char _inFile_alipos[256]
Definition: HIPplots.h:34
Gflash::Zmin
const double Zmin[kNumberCalorimeter]
Definition: GflashNameSpace.h:27
hlt_dqm_clientPB-live_cfg.xmax
tuple xmax
Definition: hlt_dqm_clientPB-live_cfg.py:37
HIPplots::_inFile_params
char _inFile_params[256]
Definition: HIPplots.h:31
HIPplots::peakthreshold
float peakthreshold
Definition: HIPplots.h:49
fed_dqm_sourceclient-live_cfg.path
string path
Definition: fed_dqm_sourceclient-live_cfg.py:33
HIPplots::_inFile_HIPalign
char _inFile_HIPalign[256]
Definition: HIPplots.h:36
benchmark_cfg.fa
string fa
Definition: benchmark_cfg.py:13
HIPplots::extractAlignParams
void extractAlignParams(int i, int minHits=0, int subDet=0, int doubleSided=0)
Definition: HIPplots.cc:72
HIPplots::plotAlignParamsAtIter
void plotAlignParamsAtIter(int iter, string ShiftsOrParams, char *plotName="test.png")
Definition: HIPplots.cc:528
a
double a
Definition: hdecay.h:121
lut2db_cfg.filename
tuple filename
Definition: lut2db_cfg.py:20
edm::shift
static unsigned int const shift
Definition: LuminosityBlockID.cc:9
HIPplots::plotHitMap
void plotHitMap(char *outpath, int subDet, int minHits=0)
Definition: HIPplots.cc:1050
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:145
cuy.col
int col
Definition: cuy.py:1008
HIPplots::_inFile_uservars
char _inFile_uservars[256]
Definition: HIPplots.h:32
COLOR_CODE
#define COLOR_CODE(icolor)
HIPplots::_outFile
char _outFile[256]
Definition: HIPplots.h:30
hlt_dqm_clientPB-live_cfg.xmin
tuple xmin
Definition: hlt_dqm_clientPB-live_cfg.py:36
HIPplots::FindPeaks
int FindPeaks(TH1F *h1, float *peaklist, const int maxNpeaks, int startbin=-1, int endbin=-1)
Definition: HIPplots.cc:1486
HIPplots::SetPeakThreshold
void SetPeakThreshold(float newpeakthreshold)
Definition: HIPplots.cc:1526
Generated for CMSSW Reference Manual by doxygen 1.8.5