CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiStripGainFromCalibTree.cc
Go to the documentation of this file.
1 // Original Author: Loic QUERTENMONT
2 // Created: Mon Nov 16 08:55:18 CET 2009
3 
4 #include <memory>
5 #include <iostream>
6 
16 
22 
29 
31 
34 
47 
50 
53 
58 
61 
64 
65 
66 #include "TFile.h"
67 #include "TObjString.h"
68 #include "TString.h"
69 #include "TH1F.h"
70 #include "TH2F.h"
71 #include "TProfile.h"
72 #include "TF1.h"
73 #include "TROOT.h"
74 #include "TTree.h"
75 #include "TChain.h"
76 
77 #include <ext/hash_map>
78 
79 using namespace edm;
80 using namespace reco;
81 using namespace std;
82 using namespace __gnu_cxx;
83 using __gnu_cxx::hash_map;
84 using __gnu_cxx::hash;
85 
86 struct stAPVGain{
87  unsigned int Index;
88  int Bin;
89  unsigned int DetId;
90  unsigned int APVId;
91  unsigned int SubDet;
92  float x;
93  float y;
94  float z;
95  float Eta;
96  float R;
97  float Phi;
98  float Thickness;
99  double FitMPV;
100  double FitMPVErr;
101  double FitWidth;
102  double FitWidthErr;
103  double FitChi2;
104  double FitNorm;
105  double Gain;
106  double CalibGain;
107  double PreviousGain;
109  double NEntries;
110  TH1F* HCharge;
111  TH1F* HChargeN;
112  bool isMasked;
113 };
114 
115 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
116 public:
119 
120 
121 private:
122 
123 
124  virtual void algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
125  virtual void algoEndRun (const edm::Run& run, const edm::EventSetup& iSetup) override;
126  virtual void algoBeginJob (const edm::EventSetup& iSetup) override;
127  virtual void algoEndJob () override;
128  virtual void algoAnalyze (const edm::Event &, const edm::EventSetup &) override;
129 
130  void merge(TH2* A, TH2* B); //needed to add histograms with different number of bins
131  void algoAnalyzeTheTree();
132  void algoComputeMPVandGain();
133  void processEvent(); //what really does the job
134 
135  void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
136  bool IsGoodLandauFit(double* FitResults);
137  void storeOnTree(TFileService* tfs);
138  void MakeCalibrationMap();
139  bool produceTagFilter();
140 
141  template<typename T>
142  inline edm::Handle<T> connect(const T* &ptr, edm::EDGetTokenT<T> token, const edm::Event &evt) {
144  evt.getByToken(token, handle);
145  ptr = handle.product();
146  return handle; //return handle to keep alive pointer (safety first)
147  }
148 
149  SiStripApvGain* getNewObject() override;
150 
154  double MinNrEntries;
155  double MaxMPVError;
159  double MinTrackEta;
160  double MaxTrackEta;
161  unsigned int MaxNrStrips;
162  unsigned int MinTrackHits;
169 
173 
176 
179  vector<string> VInputFiles;
180 
191 
192  unsigned int NEvent;
193  unsigned int NTrack;
194  unsigned int NClusterStrip;
195  unsigned int NClusterPixel;
198  unsigned int SRun;
199  unsigned int ERun;
200  unsigned int GOOD;
201  unsigned int BAD;
202  unsigned int MASKED;
203 
204  //data members for processing
205  unsigned int eventnumber = 0;
206  unsigned int runnumber = 0;
207  const std::vector<bool>* TrigTech = 0; edm::EDGetTokenT<std::vector<bool> > TrigTech_token_;
208 
209  const std::vector<double>* trackchi2ndof = 0; edm::EDGetTokenT<std::vector<double> > trackchi2ndof_token_;
210  const std::vector<float>* trackp = 0; edm::EDGetTokenT<std::vector<float> > trackp_token_;
211  const std::vector<float>* trackpt = 0; edm::EDGetTokenT<std::vector<float> > trackpt_token_;
212  const std::vector<double>* tracketa = 0; edm::EDGetTokenT<std::vector<double> > tracketa_token_;
213  const std::vector<double>* trackphi = 0; edm::EDGetTokenT<std::vector<double> > trackphi_token_;
214  const std::vector<unsigned int>* trackhitsvalid = 0; edm::EDGetTokenT<std::vector<unsigned int> > trackhitsvalid_token_;
215 
216  const std::vector<int>* trackindex = 0; edm::EDGetTokenT<std::vector<int> > trackindex_token_;
217  const std::vector<unsigned int>* rawid = 0; edm::EDGetTokenT<std::vector<unsigned int> > rawid_token_;
218  const std::vector<double>* localdirx = 0; edm::EDGetTokenT<std::vector<double> > localdirx_token_;
219  const std::vector<double>* localdiry = 0; edm::EDGetTokenT<std::vector<double> > localdiry_token_;
220  const std::vector<double>* localdirz = 0; edm::EDGetTokenT<std::vector<double> > localdirz_token_;
221  const std::vector<unsigned short>* firststrip = 0; edm::EDGetTokenT<std::vector<unsigned short> > firststrip_token_;
222  const std::vector<unsigned short>* nstrips = 0; edm::EDGetTokenT<std::vector<unsigned short> > nstrips_token_;
223  const std::vector<bool>* saturation = 0; edm::EDGetTokenT<std::vector<bool> > saturation_token_;
224  const std::vector<bool>* overlapping = 0; edm::EDGetTokenT<std::vector<bool> > overlapping_token_;
225  const std::vector<bool>* farfromedge = 0; edm::EDGetTokenT<std::vector<bool> > farfromedge_token_;
226  const std::vector<unsigned int>* charge = 0; edm::EDGetTokenT<std::vector<unsigned int> > charge_token_;
227  const std::vector<double>* path = 0; edm::EDGetTokenT<std::vector<double> > path_token_;
228  const std::vector<double>* chargeoverpath = 0; edm::EDGetTokenT<std::vector<double> > chargeoverpath_token_;
229  const std::vector<unsigned char>* amplitude = 0; edm::EDGetTokenT<std::vector<unsigned char> > amplitude_token_;
230  const std::vector<double>* gainused = 0; edm::EDGetTokenT<std::vector<double> > gainused_token_;
231 
232  string EventPrefix_; //("");
233  string EventSuffix_; //("");
234  string TrackPrefix_; //("track");
235  string TrackSuffix_; //("");
236  string CalibPrefix_; //("GainCalibration");
237  string CalibSuffix_; //("");
238  string tree_path_;
239 
240 private :
241  class isEqual{
242  public:
243  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
244  };
245 
246  std::vector<stAPVGain*> APVsCollOrdered;
247  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
248 };
249 
251  if(A->GetNbinsX() == B->GetNbinsX()){
252  A->Add(B);
253  }else{
254  for(int x=0;x<=B->GetNbinsX()+1; x++){
255  for(int y=0;y<=B->GetNbinsY()+1; y++){
256  A->SetBinContent(x,y,A->GetBinContent(x,y)+B->GetBinContent(x,y));
257  }}
258  }
259 }
260 
261 
262 
264 {
265  Charge_Vs_Index = NULL;//make sure this is initialized to NULL
266  OutputGains = iConfig.getParameter<std::string>("OutputGains");
267 
268  AlgoMode = iConfig.getUntrackedParameter<std::string>("AlgoMode", "CalibTree");
269  MinNrEntries = iConfig.getUntrackedParameter<double> ("minNrEntries" , 20);
270  MaxMPVError = iConfig.getUntrackedParameter<double> ("maxMPVError" , 500.0);
271  MaxChi2OverNDF = iConfig.getUntrackedParameter<double> ("maxChi2OverNDF" , 5.0);
272  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 3.0);
273  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
274  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
275  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
276  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 2);
277  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 8);
278  MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double> ("MaxTrackChiOverNdf" , 3);
279  AllowSaturation = iConfig.getUntrackedParameter<bool> ("AllowSaturation" , false);
280  FirstSetOfConstants = iConfig.getUntrackedParameter<bool> ("FirstSetOfConstants", true);
281  Validation = iConfig.getUntrackedParameter<bool> ("Validation" , false);
282  OldGainRemoving = iConfig.getUntrackedParameter<bool> ("OldGainRemoving" , false);
283 
284  CalibrationLevel = iConfig.getUntrackedParameter<int> ("CalibrationLevel" , 0);
285  VInputFiles = iConfig.getUntrackedParameter<vector<string> > ("InputFiles");
286 
287 
288  useCalibration = iConfig.getUntrackedParameter<bool> ("UseCalibration" , false);
289  m_calibrationPath = iConfig.getUntrackedParameter<string> ("calibrationPath");
290  harvestingMode = iConfig.getUntrackedParameter<bool> ("harvestingMode" , false);
291 
292  tagCondition_NClusters = iConfig.getUntrackedParameter<double> ("NClustersForTagProd" , 2E8);
293  tagCondition_GoodFrac = iConfig.getUntrackedParameter<double> ("GoodFracForTagProd" , 0.95);
294 
295  saveSummary = iConfig.getUntrackedParameter<bool> ("saveSummary" , false);
296  tree_path_ = iConfig.getUntrackedParameter<string> ("treePath");
298  dbe->setVerbose(10);
299 
300  edm::ParameterSet swhallowgain_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("gain");
301  string label = swhallowgain_pset.getUntrackedParameter<string>("label");
302  CalibPrefix_ = swhallowgain_pset.getUntrackedParameter<string>("prefix");
303  CalibSuffix_ = swhallowgain_pset.getUntrackedParameter<string>("suffix");
304 
305  trackindex_token_ = consumes<std::vector<int> >(edm::InputTag(label, CalibPrefix_ + "trackindex" + CalibSuffix_));
306  rawid_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, CalibPrefix_ + "rawid" + CalibSuffix_));
307  localdirx_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdirx" + CalibSuffix_));
308  localdiry_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdiry" + CalibSuffix_));
309  localdirz_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdirz" + CalibSuffix_));
310  firststrip_token_ = consumes<std::vector<unsigned short> >(edm::InputTag(label, CalibPrefix_ + "firststrip" + CalibSuffix_));
311  nstrips_token_ = consumes<std::vector<unsigned short> >(edm::InputTag(label, CalibPrefix_ + "nstrips" + CalibSuffix_));
312  saturation_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "saturation" + CalibSuffix_));
313  overlapping_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "overlapping" + CalibSuffix_));
314  farfromedge_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "farfromedge" + CalibSuffix_));
315  charge_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, CalibPrefix_ + "charge" + CalibSuffix_));
316  path_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "path" + CalibSuffix_));
317  chargeoverpath_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "chargeoverpath"+ CalibSuffix_));
318  amplitude_token_ = consumes<std::vector<unsigned char> >(edm::InputTag(label, CalibPrefix_ + "amplitude" + CalibSuffix_));
319  gainused_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "gainused" + CalibSuffix_));
320 
321  edm::ParameterSet evtinfo_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("evtinfo");
322  label = evtinfo_pset.getUntrackedParameter<string>("label");
323  EventPrefix_ = evtinfo_pset.getUntrackedParameter<string>("prefix");
324  EventSuffix_ = evtinfo_pset.getUntrackedParameter<string>("suffix");
325  TrigTech_token_ = consumes<std::vector<bool> >(edm::InputTag(label, EventPrefix_ + "TrigTech" + EventSuffix_));
326 
327  edm::ParameterSet track_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("tracks");
328  label = track_pset.getUntrackedParameter<string>("label");
329  TrackPrefix_ = track_pset.getUntrackedParameter<string>("prefix");
330  TrackSuffix_ = track_pset.getUntrackedParameter<string>("suffix");
331 
332  trackchi2ndof_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "chi2ndof" + TrackSuffix_));
333  trackp_token_ = consumes<std::vector<float> >(edm::InputTag(label, TrackPrefix_ + "momentum" + TrackSuffix_));
334  trackpt_token_ = consumes<std::vector<float> >(edm::InputTag(label, TrackPrefix_ + "pt" + TrackSuffix_));
335  tracketa_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "eta" + TrackSuffix_));
336  trackphi_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "phi" + TrackSuffix_));
337  trackhitsvalid_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, TrackPrefix_ + "hitsvalid" + TrackSuffix_));
338 }
339 
340 
342 {
343  dbe->setCurrentFolder("AlCaReco/SiStripGains/");
344  if(AlgoMode != "PCL" or harvestingMode)dbe->setCurrentFolder("AlCaReco/SiStripGainsHarvesting/");
345 
346  Charge_Vs_Index = dbe->book2D("Charge_Vs_Index" , "Charge_Vs_Index" , 88625, 0 , 88624,2000,0,4000);
347  Charge_Vs_Index_Absolute = dbe->book2D("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 88625, 0 , 88624,1000,0,4000);
348 // Charge_Vs_Index = dbe->book2D("Charge_Vs_Index" , "Charge_Vs_Index" , 72785, 0 , 72784,1000,0,2000);
349 // Charge_Vs_Index_Absolute = dbe->book2D("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 72785, 0 , 72784, 500,0,2000);
350  Charge_Vs_PathlengthTIB = dbe->book2D("Charge_Vs_PathlengthTIB" , "Charge_Vs_PathlengthTIB" , 20 , 0.3 , 1.3 , 250,0,2000);
351  Charge_Vs_PathlengthTOB = dbe->book2D("Charge_Vs_PathlengthTOB" , "Charge_Vs_PathlengthTOB" , 20 , 0.3 , 1.3 , 250,0,2000);
352  Charge_Vs_PathlengthTIDP = dbe->book2D("Charge_Vs_PathlengthTIDP" , "Charge_Vs_PathlengthTIDP" , 20 , 0.3 , 1.3 , 250,0,2000);
353  Charge_Vs_PathlengthTIDM = dbe->book2D("Charge_Vs_PathlengthTIDM" , "Charge_Vs_PathlengthTIDM" , 20 , 0.3 , 1.3 , 250,0,2000);
354  Charge_Vs_PathlengthTECP1 = dbe->book2D("Charge_Vs_PathlengthTECP1", "Charge_Vs_PathlengthTECP1", 20 , 0.3 , 1.3 , 250,0,2000);
355  Charge_Vs_PathlengthTECP2 = dbe->book2D("Charge_Vs_PathlengthTECP2", "Charge_Vs_PathlengthTECP2", 20 , 0.3 , 1.3 , 250,0,2000);
356  Charge_Vs_PathlengthTECM1 = dbe->book2D("Charge_Vs_PathlengthTECM1", "Charge_Vs_PathlengthTECM1", 20 , 0.3 , 1.3 , 250,0,2000);
357  Charge_Vs_PathlengthTECM2 = dbe->book2D("Charge_Vs_PathlengthTECM2", "Charge_Vs_PathlengthTECM2", 20 , 0.3 , 1.3 , 250,0,2000);
358 
360  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
361  auto const & Det = tkGeom->dets();
362 
363  NPixelDets = 0;
364  NStripAPVs = 0;
365  unsigned int Index=0;
366  for(unsigned int i=0;i<Det.size();i++){
367  DetId Detid = Det[i]->geographicalId();
368  int SubDet = Detid.subdetId();
369 
370  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
371  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
372 
373  auto DetUnit = dynamic_cast<const StripGeomDetUnit*> (Det[i]);
374  if(!DetUnit)continue;
375 
376  const StripTopology& Topo = DetUnit->specificTopology();
377  unsigned int NAPV = Topo.nstrips()/128;
378 
379  for(unsigned int j=0;j<NAPV;j++){
380  stAPVGain* APV = new stAPVGain;
381  APV->Index = Index;
382  APV->Bin = -1;
383  APV->DetId = Detid.rawId();
384  APV->APVId = j;
385  APV->SubDet = SubDet;
386  APV->FitMPV = -1;
387  APV->FitMPVErr = -1;
388  APV->FitWidth = -1;
389  APV->FitWidthErr = -1;
390  APV->FitChi2 = -1;
391  APV->FitNorm = -1;
392  APV->Gain = -1;
393  APV->PreviousGain = 1;
394  APV->PreviousGainTick = 1;
395  APV->x = DetUnit->position().basicVector().x();
396  APV->y = DetUnit->position().basicVector().y();
397  APV->z = DetUnit->position().basicVector().z();
398  APV->Eta = DetUnit->position().basicVector().eta();
399  APV->Phi = DetUnit->position().basicVector().phi();
400  APV->R = DetUnit->position().basicVector().transverse();
401  APV->Thickness = DetUnit->surface().bounds().thickness();
402  APV->NEntries = 0;
403  APV->isMasked = false;
404 
405  APVsCollOrdered.push_back(APV);
406  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
407  Index++;
408  NStripAPVs++;
409  }
410  }
411  }
412 
413  for(unsigned int i=0;i<Det.size();i++){ //Make two loop such that the Pixel information is added at the end --> make transition simpler
414  DetId Detid = Det[i]->geographicalId();
415  int SubDet = Detid.subdetId();
417  auto DetUnit = dynamic_cast<const PixelGeomDetUnit*> (Det[i]);
418  if(!DetUnit) continue;
419 
420  const PixelTopology& Topo = DetUnit->specificTopology();
421  unsigned int NROCRow = Topo.nrows()/(80.);
422  unsigned int NROCCol = Topo.ncolumns()/(52.);
423 
424  for(unsigned int j=0;j<NROCRow;j++){
425  for(unsigned int i=0;i<NROCCol;i++){
426 
427  stAPVGain* APV = new stAPVGain;
428  APV->Index = Index;
429  APV->Bin = -1;
430  APV->DetId = Detid.rawId();
431  APV->APVId = (j<<3 | i);
432  APV->SubDet = SubDet;
433  APV->FitMPV = -1;
434  APV->FitMPVErr = -1;
435  APV->FitWidth = -1;
436  APV->FitWidthErr = -1;
437  APV->FitChi2 = -1;
438  APV->Gain = -1;
439  APV->PreviousGain = 1;
440  APV->x = DetUnit->position().basicVector().x();
441  APV->y = DetUnit->position().basicVector().y();
442  APV->z = DetUnit->position().basicVector().z();
443  APV->Eta = DetUnit->position().basicVector().eta();
444  APV->Phi = DetUnit->position().basicVector().phi();
445  APV->R = DetUnit->position().basicVector().transverse();
446  APV->Thickness = DetUnit->surface().bounds().thickness();
447  APV->isMasked = false; //SiPixelQuality_->IsModuleBad(Detid.rawId());
448  APV->NEntries = 0;
449 
450  APVsCollOrdered.push_back(APV);
451  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
452  Index++;
453  NPixelDets++;
454  }}
455  }
456  }
457 
458 
460 
461  NEvent = 0;
462  NTrack = 0;
463  NClusterStrip = 0;
464  NClusterPixel = 0;
465  SRun = 1<<31;
466  ERun = 0;
467  GOOD = 0;
468  BAD = 0;
469  MASKED = 0;
470 }
471 
472 
473 
475 {
476  edm::ESHandle<SiStripGain> gainHandle;
477  iSetup.get<SiStripGainRcd>().get(gainHandle);
478  if(!gainHandle.isValid()){edm::LogError("SiStripGainFromCalibTree")<< "gainHandle is not valid\n"; exit(0);}
479 
480  edm::ESHandle<SiStripQuality> SiStripQuality_;
481  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
482 
483  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
485 
486  APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId,APV->APVId);
487 // if(!FirstSetOfConstants){
488  if(gainHandle->getNumberOfTags()!=2){edm::LogError("SiStripGainFromCalibTree")<< "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";fflush(stdout);exit(0);};
489  float newPreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
490  if(APV->PreviousGain!=1 and newPreviousGain!=APV->PreviousGain)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: ParticleGain in the global tag changed\n";
491  APV->PreviousGain = newPreviousGain;
492 
493  float newPreviousGainTick = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 0),1);
494  if(APV->PreviousGainTick!=1 and newPreviousGainTick!=APV->PreviousGainTick)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: TickMarkGain in the global tag changed\n";
495  APV->PreviousGainTick = newPreviousGainTick;
496 
497 
498  //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f (G1) x %8.4f (G2)\n",APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain);
499 
500 
501 
502 // }
503  }
504 }
505 
507  if(AlgoMode == "PCL" && !harvestingMode)return;//nothing to do in that case
508 
509  if(AlgoMode == "PCL" and harvestingMode){
510  // Load the 2D histograms from the DQM objects
511  // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
512  // DQMStore so that they can be used in the fit
513 
514 // edm::LogInfo("SiStripGainFromCalibTree")<< "Harvesting " << (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F()->GetEntries() << " more clusters\n";
515 
516  merge(Charge_Vs_Index ->getTH2F(), (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F() );
517 // Charge_Vs_Index ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F());
518  merge(Charge_Vs_Index_Absolute ->getTH2F(), (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index_Absolute"))->getTH2F() );
519 // Charge_Vs_Index_Absolute ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index_Absolute"))->getTH2F());
520  Charge_Vs_PathlengthTIB ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIB"))->getTH2F());
521  Charge_Vs_PathlengthTOB ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTOB"))->getTH2F());
522  Charge_Vs_PathlengthTIDP ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDP"))->getTH2F());
523  Charge_Vs_PathlengthTIDM ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDM"))->getTH2F());
524  Charge_Vs_PathlengthTECP1 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP1"))->getTH2F());
525  Charge_Vs_PathlengthTECP2 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP2"))->getTH2F());
526  Charge_Vs_PathlengthTECM1 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM1"))->getTH2F());
527  Charge_Vs_PathlengthTECM2 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM2"))->getTH2F());
528  }
529 }
530 void
532  if(AlgoMode == "PCL" && !harvestingMode)return;//nothing to do in that case
533 
534  if(AlgoMode == "CalibTree"){
535  // Loop on calibTrees to fill the 2D histograms
537  }else if(harvestingMode){
538  NClusterStrip = Charge_Vs_Index->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 );
539  NClusterPixel = Charge_Vs_Index->getTH2F()->Integral(NStripAPVs+2, NStripAPVs+NPixelDets+2, 0, 99999 );
540  }
541 
542  // Now that we have the full statistics we can extract the information of the 2D histograms
544 
545  if(AlgoMode != "PCL" or saveSummary){
546  //also save the 2D monitor elements to this file as TH2D tfs
548  tfs->make<TH2F> (*Charge_Vs_Index->getTH2F());
558 
559 
560  storeOnTree(tfs);
561  }
562 }
563 
564 
565 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
566 {
567  FitResults[0] = -0.5; //MPV
568  FitResults[1] = 0; //MPV error
569  FitResults[2] = -0.5; //Width
570  FitResults[3] = 0; //Width error
571  FitResults[4] = -0.5; //Fit Chi2/NDF
572  FitResults[5] = 0; //Normalization
573 
574  if( InputHisto->GetEntries() < MinNrEntries)return;
575 
576  // perform fit with standard landau
577  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
578  MyLandau->SetParameter(1,300);
579  InputHisto->Fit(MyLandau,"0QR WW");
580 
581  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
582  FitResults[0] = MyLandau->GetParameter(1); //MPV
583  FitResults[1] = MyLandau->GetParError(1); //MPV error
584  FitResults[2] = MyLandau->GetParameter(2); //Width
585  FitResults[3] = MyLandau->GetParError(2); //Width error
586  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
587  FitResults[5] = MyLandau->GetParameter(0);
588 
589  delete MyLandau;
590 }
591 
593  if(FitResults[0] <= 0 )return false;
594 // if(FitResults[1] > MaxMPVError )return false;
595 // if(FitResults[4] > MaxChi2OverNDF)return false;
596  return true;
597 }
598 
600 
603 
604  NEvent++;
605  NTrack+=(*trackp).size();
606 
607  unsigned int FirstAmplitude=0;
608  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
609  FirstAmplitude+=(*nstrips)[i];
610  int TI = (*trackindex)[i];
611  //printf("%i - %i - %i %i %i\n", (int)(*rawid)[i], (int)(*firststrip)[i]/128, (int)(*farfromedge)[i], (int)(*overlapping)[i], (int)(*saturation )[i] );
612  if((*tracketa )[TI] < MinTrackEta )continue;
613  if((*tracketa )[TI] > MaxTrackEta )continue;
614  if((*trackp )[TI] < MinTrackMomentum )continue;
615  if((*trackp )[TI] > MaxTrackMomentum )continue;
616  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
617  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
618 
619  stAPVGain* APV = APVsColl[((*rawid)[i]<<4) | ((*firststrip)[i]/128)]; //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
620 
621  if(APV->SubDet>2 && (*farfromedge)[i] == false )continue;
622  if(APV->SubDet>2 && (*overlapping)[i] == true )continue;
623  if(APV->SubDet>2 && (*saturation )[i] && !AllowSaturation)continue;
624  if(APV->SubDet>2 && (*nstrips )[i] > MaxNrStrips )continue;
625 
626 
627  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
628 
629  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
630  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
631  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
632 
633  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
634  //printf("Amplitudes: ");
635  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
636  //printf("\n");
637 
638  if(APV->SubDet>2){NClusterStrip++;}else{NClusterPixel++;}
639 
640  int Charge = 0;
641  if(APV->SubDet>2 && (useCalibration || !FirstSetOfConstants)){
642  bool Saturation = false;
643  for(unsigned int s=0;s<(*nstrips)[i];s++){
644  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
645  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
646  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
647  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
648  if(StripCharge>1024){
649  StripCharge = 255;
650  Saturation = true;
651  }else if(StripCharge>254){
652  StripCharge = 254;
653  Saturation = true;
654  }
655  Charge += StripCharge;
656  }
657  if(Saturation && !AllowSaturation)continue;
658  }else if(APV->SubDet>2){
659  Charge = (*charge)[i];
660  }else{
661  Charge = (*charge)[i]/265.0; //expected scale factor between pixel and strip charge
662  }
663 
664  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
665 
666  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
667  if(APV->SubDet>2){
668  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
669  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
670  }
671  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
672  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
673 
675  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill((*path)[i],Charge);
676  }else if(APV->SubDet==StripSubdetector::TID){
677  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
678  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
679  }
680  }else if(APV->SubDet==StripSubdetector::TEC){
681  if(APV->Eta<0){
682  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
683  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
684  }
685  }else if(APV->Eta>0){
686  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
687  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
688  }
689  }
690  }
691 
692  }// END OF ON-CLUSTER LOOP
693 }//END OF processEvent()
694 
696 {
697  for(unsigned int i=0;i<VInputFiles.size();i++){
698  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
699  TFile *tfile = TFile::Open(VInputFiles[i].c_str());
700  TTree* tree = dynamic_cast<TTree*> (tfile->Get(tree_path_.c_str()));
701 
702  tree->SetBranchAddress((EventPrefix_ + "event" + EventSuffix_).c_str(), &eventnumber , NULL);
703  tree->SetBranchAddress((EventPrefix_ + "run" + EventSuffix_).c_str(), &runnumber , NULL);
704  tree->SetBranchAddress((EventPrefix_ + "TrigTech" + EventSuffix_).c_str(), &TrigTech , NULL);
705 
706  tree->SetBranchAddress((TrackPrefix_ + "chi2ndof" + TrackSuffix_).c_str(), &trackchi2ndof , NULL);
707  tree->SetBranchAddress((TrackPrefix_ + "momentum" + TrackSuffix_).c_str(), &trackp , NULL);
708  tree->SetBranchAddress((TrackPrefix_ + "pt" + TrackSuffix_).c_str(), &trackpt , NULL);
709  tree->SetBranchAddress((TrackPrefix_ + "eta" + TrackSuffix_).c_str(), &tracketa , NULL);
710  tree->SetBranchAddress((TrackPrefix_ + "phi" + TrackSuffix_).c_str(), &trackphi , NULL);
711  tree->SetBranchAddress((TrackPrefix_ + "hitsvalid" + TrackSuffix_).c_str(), &trackhitsvalid, NULL);
712 
713  tree->SetBranchAddress((CalibPrefix_ + "trackindex" + CalibSuffix_).c_str(), &trackindex , NULL);
714  tree->SetBranchAddress((CalibPrefix_ + "rawid" + CalibSuffix_).c_str(), &rawid , NULL);
715  tree->SetBranchAddress((CalibPrefix_ + "localdirx" + CalibSuffix_).c_str(), &localdirx , NULL);
716  tree->SetBranchAddress((CalibPrefix_ + "localdiry" + CalibSuffix_).c_str(), &localdiry , NULL);
717  tree->SetBranchAddress((CalibPrefix_ + "localdirz" + CalibSuffix_).c_str(), &localdirz , NULL);
718  tree->SetBranchAddress((CalibPrefix_ + "firststrip" + CalibSuffix_).c_str(), &firststrip , NULL);
719  tree->SetBranchAddress((CalibPrefix_ + "nstrips" + CalibSuffix_).c_str(), &nstrips , NULL);
720  tree->SetBranchAddress((CalibPrefix_ + "saturation" + CalibSuffix_).c_str(), &saturation , NULL);
721  tree->SetBranchAddress((CalibPrefix_ + "overlapping" + CalibSuffix_).c_str(), &overlapping , NULL);
722  tree->SetBranchAddress((CalibPrefix_ + "farfromedge" + CalibSuffix_).c_str(), &farfromedge , NULL);
723  tree->SetBranchAddress((CalibPrefix_ + "charge" + CalibSuffix_).c_str(), &charge , NULL);
724  tree->SetBranchAddress((CalibPrefix_ + "path" + CalibSuffix_).c_str(), &path , NULL);
725  tree->SetBranchAddress((CalibPrefix_ + "chargeoverpath" + CalibSuffix_).c_str(), &chargeoverpath, NULL);
726  tree->SetBranchAddress((CalibPrefix_ + "amplitude" + CalibSuffix_).c_str(), &amplitude , NULL);
727  tree->SetBranchAddress((CalibPrefix_ + "gainused" + CalibSuffix_).c_str(), &gainused , NULL);
728 
729 
730  unsigned int nentries = tree->GetEntries();
731  printf("Number of Events = %i + %i = %i\n",NEvent,nentries,(NEvent+nentries));
732  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
733  printf("Looping on the Tree :");
734  int TreeStep = nentries/50;if(TreeStep<=1)TreeStep=1;
735  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
736  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
737  tree->GetEntry(ientry);
738  processEvent();
739  }printf("\n");// END OF EVENT LOOP
740  }
741 }
742 
743 
744 
746  unsigned int I=0;
747  TH1F* Proj = NULL;
748  double FitResults[6];
749  double MPVmean = 300;
750 
751  TH2F *chvsidx = Charge_Vs_Index->getTH2F();
752 
753 
754  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
755  printf("Fitting Charge Distribution :");
756  int TreeStep = APVsColl.size()/50;
757  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
758  if(I%TreeStep==0){printf(".");fflush(stdout);}
759  stAPVGain* APV = it->second;
760  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
761 
762  if(APV->isMasked){APV->Gain=APV->PreviousGain; MASKED++; continue;}
763 
764  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
765  if(!Proj)continue;
766 
767  if(CalibrationLevel==0){
768  }else if(CalibrationLevel==1){
769  int SecondAPVId = APV->APVId;
770  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
771  stAPVGain* APV2 = APVsColl[(APV->DetId<<4) | SecondAPVId];
772  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
773  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
774  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
775  }else if(CalibrationLevel==2){
776  for(unsigned int i=0;i<16;i++){ //loop up to 6APV for Strip and up to 16 for Pixels
777  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
778  tmpit = APVsColl.find((APV->DetId<<4) | i);
779  if(tmpit==APVsColl.end())continue;
780  stAPVGain* APV2 = tmpit->second;
781  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
782  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
783  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
784  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
785  }
786  }else{
787  CalibrationLevel = 0;
788  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
789  }
790 
791  getPeakOfLandau(Proj,FitResults);
792  APV->FitMPV = FitResults[0];
793  APV->FitMPVErr = FitResults[1];
794  APV->FitWidth = FitResults[2];
795  APV->FitWidthErr = FitResults[3];
796  APV->FitChi2 = FitResults[4];
797  APV->FitNorm = FitResults[5];
798  APV->NEntries = Proj->GetEntries();
799 
800  if(IsGoodLandauFit(FitResults)){
801  APV->Gain = APV->FitMPV / MPVmean;
802  if(APV->SubDet>2)GOOD++;
803  }else{
804  APV->Gain = APV->PreviousGain;
805  if(APV->SubDet>2)BAD++;
806  }
807  if(APV->Gain<=0) APV->Gain = 1;
808 
809  //printf("%5i/%5i: %6i - %1i %5E Entries --> MPV = %f +- %f\n",I,APVsColl.size(),APV->DetId, APV->APVId, Proj->GetEntries(), FitResults[0], FitResults[1]);fflush(stdout);
810  delete Proj;
811  }printf("\n");
812 }
813 
814 
816 {
817  unsigned int tree_Index;
818  unsigned int tree_Bin;
819  unsigned int tree_DetId;
820  unsigned char tree_APVId;
821  unsigned char tree_SubDet;
822  float tree_x;
823  float tree_y;
824  float tree_z;
825  float tree_Eta;
826  float tree_R;
827  float tree_Phi;
828  float tree_Thickness;
829  float tree_FitMPV;
830  float tree_FitMPVErr;
831  float tree_FitWidth;
832  float tree_FitWidthErr;
833  float tree_FitChi2NDF;
834  float tree_FitNorm;
835  double tree_Gain;
836  double tree_PrevGain;
837  double tree_PrevGainTick;
838  double tree_NEntries;
839  bool tree_isMasked;
840 
841  TTree* MyTree;
842  MyTree = tfs->make<TTree> ("APVGain","APVGain");
843  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
844  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
845  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
846  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
847  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
848  MyTree->Branch("x" ,&tree_x ,"x/F");
849  MyTree->Branch("y" ,&tree_y ,"y/F");
850  MyTree->Branch("z" ,&tree_z ,"z/F");
851  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
852  MyTree->Branch("R" ,&tree_R ,"R/F");
853  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
854  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
855  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
856  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
857  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
858  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
859  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
860  MyTree->Branch("FitNorm" ,&tree_FitNorm ,"FitNorm/F");
861  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
862  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
863  MyTree->Branch("PrevGainTick" ,&tree_PrevGainTick,"PrevGainTick/D");
864  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
865  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
866 
867 
868  FILE* Gains = stdout;
869  fprintf(Gains,"NEvents = %i\n",NEvent);
870  fprintf(Gains,"NTracks = %i\n",NTrack);
871  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
872  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
873  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
874  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
875  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
876 
877  Gains=fopen(OutputGains.c_str(),"w");
878  fprintf(Gains,"NEvents = %i\n",NEvent);
879  fprintf(Gains,"NTracks = %i\n",NTrack);
880  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
881  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
882  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
883  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
884  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
885 
886  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
888  if(APV==NULL)continue;
889 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
890  fprintf(Gains,"%i | %i | PreviousGain = %7.5f(tick) x %7.5f(particle) NewGain (particle) = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain,APV->Gain, APV->NEntries);
891 
892  tree_Index = APV->Index;
893  tree_Bin = Charge_Vs_Index->getTH2F()->GetXaxis()->FindBin(APV->Index);
894  tree_DetId = APV->DetId;
895  tree_APVId = APV->APVId;
896  tree_SubDet = APV->SubDet;
897  tree_x = APV->x;
898  tree_y = APV->y;
899  tree_z = APV->z;
900  tree_Eta = APV->Eta;
901  tree_R = APV->R;
902  tree_Phi = APV->Phi;
903  tree_Thickness = APV->Thickness;
904  tree_FitMPV = APV->FitMPV;
905  tree_FitMPVErr = APV->FitMPVErr;
906  tree_FitWidth = APV->FitWidth;
907  tree_FitWidthErr= APV->FitWidthErr;
908  tree_FitChi2NDF = APV->FitChi2;
909  tree_FitNorm = APV->FitNorm;
910  tree_Gain = APV->Gain;
911  tree_PrevGain = APV->PreviousGain;
912  tree_PrevGainTick = APV->PreviousGainTick;
913  tree_NEntries = APV->NEntries;
914  tree_isMasked = APV->isMasked;
915 
916 
917  if(tree_DetId==402673324){
918  printf("%i | %i : %f --> %f (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
919  }
920 
921 
922  MyTree->Fill();
923  }
924  if(Gains)fclose(Gains);
925 
926 
927 }
928 
930 
931  // The goal of this function is to check wether or not there is enough statistics to produce a meaningful tag for the DB or not
932  if(Charge_Vs_Index->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 ) < tagCondition_NClusters) {
933  edm::LogWarning("SiStripGainFromCalibTree")<< "produceTagFilter -> Return false: Statistics is too low : " << Charge_Vs_Index->getTH2F()->Integral() << endl;
934  return false;
935  }
936  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
937  edm::LogWarning("SiStripGainFromCalibTree")<< "produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
938  return false;
939  }
940  return true;
941 }
942 
944 {
946  if(!harvestingMode) return obj;
947 
948  if(!produceTagFilter()){
949  edm::LogWarning("SiStripGainFromCalibTree")<< "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
950  setDoStore(false);
951  return obj;
952  }
953 
954 
955  std::vector<float>* theSiStripVector = NULL;
956  unsigned int PreviousDetId = 0;
957  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
959  if(APV==NULL){ printf("Bug\n"); continue; }
960  if(APV->SubDet<=2)continue;
961  if(APV->DetId != PreviousDetId){
962  if(theSiStripVector!=NULL){
963  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
964  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
965  }
966  theSiStripVector = new std::vector<float>;
967  PreviousDetId = APV->DetId;
968  }
969  theSiStripVector->push_back(APV->Gain);
970  }
971  if(theSiStripVector!=NULL){
972  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
973  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
974  }
975 
976  return obj;
977 }
978 
979 
981 {
982 }
983 
985  if(!useCalibration)return;
986 
987  TChain* t1 = new TChain("SiStripCalib/APVGain");
988  t1->Add(m_calibrationPath.c_str());
989 
990  unsigned int tree_DetId;
991  unsigned char tree_APVId;
992  double tree_Gain;
993 
994  t1->SetBranchAddress("DetId" ,&tree_DetId );
995  t1->SetBranchAddress("APVId" ,&tree_APVId );
996  t1->SetBranchAddress("Gain" ,&tree_Gain );
997 
998  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
999  t1->GetEntry(ientry);
1000  stAPVGain* APV = APVsColl[(tree_DetId<<4) | (unsigned int)tree_APVId];
1001  APV->CalibGain = tree_Gain;
1002  }
1003 }
1004 
1005 void
1007 {
1008  // in AlCaHarvesting mode we just need to run the logic in the endJob step
1009  if(harvestingMode) return;
1010 
1011  if(AlgoMode=="CalibTree")return;
1012 
1013  eventnumber = iEvent.id().event();
1014  runnumber = iEvent.id().run();
1015  auto handle01 = connect(TrigTech , TrigTech_token_ , iEvent);
1016  auto handle02 = connect(trackchi2ndof , trackchi2ndof_token_ , iEvent);
1017  auto handle03 = connect(trackp , trackp_token_ , iEvent);
1018  auto handle04 = connect(trackpt , trackpt_token_ , iEvent);
1019  auto handle05 = connect(tracketa , tracketa_token_ , iEvent);
1020  auto handle06 = connect(trackphi , trackphi_token_ , iEvent);
1021  auto handle07 = connect(trackhitsvalid, trackhitsvalid_token_, iEvent);
1022  auto handle08 = connect(trackindex , trackindex_token_ , iEvent);
1023  auto handle09 = connect(rawid , rawid_token_ , iEvent);
1024  auto handle11 = connect(localdirx , localdirx_token_ , iEvent);
1025  auto handle12 = connect(localdiry , localdiry_token_ , iEvent);
1026  auto handle13 = connect(localdirz , localdirz_token_ , iEvent);
1027  auto handle14 = connect(firststrip , firststrip_token_ , iEvent);
1028  auto handle15 = connect(nstrips , nstrips_token_ , iEvent);
1029  auto handle16 = connect(saturation , saturation_token_ , iEvent);
1030  auto handle17 = connect(overlapping , overlapping_token_ , iEvent);
1031  auto handle18 = connect(farfromedge , farfromedge_token_ , iEvent);
1032  auto handle19 = connect(charge , charge_token_ , iEvent);
1033  auto handle21 = connect(path , path_token_ , iEvent);
1034  auto handle22 = connect(chargeoverpath, chargeoverpath_token_, iEvent);
1035  auto handle23 = connect(amplitude , amplitude_token_ , iEvent);
1036  auto handle24 = connect(gainused , gainused_token_ , iEvent);
1037 
1038  processEvent();
1039 }
1040 
RunNumber_t run() const
Definition: EventID.h:39
edm::EDGetTokenT< std::vector< int > > trackindex_token_
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
virtual int nstrips() const =0
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
const std::vector< float > * trackp
edm::EDGetTokenT< std::vector< double > > tracketa_token_
edm::EDGetTokenT< std::vector< unsigned char > > amplitude_token_
const std::vector< double > * chargeoverpath
edm::EDGetTokenT< std::vector< unsigned int > > rawid_token_
const std::vector< unsigned char > * amplitude
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< void, edm::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
SiStripApvGain * getNewObject() override
virtual int ncolumns() const =0
const std::vector< double > * localdirz
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
const std::vector< double > * trackchi2ndof
virtual void algoBeginRun(const edm::Run &run, const edm::EventSetup &iSetup) override
const std::vector< double > * localdirx
dispatcher processEvent(e, inputTag, standby)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual void algoEndRun(const edm::Run &run, const edm::EventSetup &iSetup) override
const std::vector< unsigned int > * rawid
edm::EDGetTokenT< std::vector< bool > > overlapping_token_
SiStripGainFromCalibTree(const edm::ParameterSet &)
edm::Handle< T > connect(const T *&ptr, edm::EDGetTokenT< T > token, const edm::Event &evt)
virtual int nrows() const =0
#define NULL
Definition: scimark2.h:8
void storeOnTree(TFileService *tfs)
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
const std::vector< float > * trackpt
const std::vector< unsigned int > * trackhitsvalid
void Fill(long long x)
edm::EDGetTokenT< std::vector< bool > > farfromedge_token_
bool put(const uint32_t &detID, Range input)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
edm::EDGetTokenT< std::vector< unsigned int > > trackhitsvalid_token_
T x() const
Cartesian x coordinate.
std::vector< stAPVGain * > APVsCollOrdered
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
edm::EDGetTokenT< std::vector< double > > path_token_
int iEvent
Definition: GenABIO.cc:230
edm::EDGetTokenT< std::vector< float > > trackp_token_
bool IsGoodLandauFit(double *FitResults)
tuple handle
Definition: patZpeak.py:22
virtual void algoBeginJob(const edm::EventSetup &iSetup) override
std::pair< ContainerIterator, ContainerIterator > Range
int j
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< std::vector< bool > > TrigTech_token_
edm::EDGetTokenT< std::vector< double > > chargeoverpath_token_
const std::complex< double > I
Definition: I.h:8
const std::vector< double > * gainused
edm::EDGetTokenT< std::vector< double > > localdirz_token_
const std::vector< unsigned short > * firststrip
const std::vector< double > * localdiry
edm::EDGetTokenT< std::vector< unsigned int > > charge_token_
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
virtual void algoEndJob() override
edm::EDGetTokenT< std::vector< unsigned short > > nstrips_token_
const std::vector< int > * trackindex
Definition: DetId.h:18
edm::EDGetTokenT< std::vector< unsigned short > > firststrip_token_
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
edm::EDGetTokenT< std::vector< float > > trackpt_token_
T const * product() const
Definition: Handle.h:81
const std::vector< bool > * overlapping
const std::vector< bool > * saturation
bool merge(LuminosityBlockRange &lh, LuminosityBlockRange &rh)
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
const T & get() const
Definition: EventSetup.h:56
const std::vector< bool > * farfromedge
const std::vector< bool > * TrigTech
void setDoStore(const bool doStore)
When set to false the payload will not be written to the db.
edm::EventID id() const
Definition: EventBase.h:59
double a
Definition: hdecay.h:121
TH2F * getTH2F(std::string name, std::string process, std::string rootfolder, DQMStore *dbe_, bool verb, bool clone)
edm::EDGetTokenT< std::vector< double > > localdirx_token_
edm::EDGetTokenT< std::vector< bool > > saturation_token_
edm::EDGetTokenT< std::vector< double > > trackchi2ndof_token_
edm::EDGetTokenT< std::vector< double > > localdiry_token_
volatile std::atomic< bool > shutdown_flag false
TH2F * getTH2F(void) const
const std::vector< unsigned int > * charge
bool isValid() const
Definition: ESHandle.h:47
const std::vector< double > * path
long double T
const std::vector< double > * trackphi
const std::vector< unsigned short > * nstrips
edm::EDGetTokenT< std::vector< double > > gainused_token_
edm::EDGetTokenT< std::vector< double > > trackphi_token_
const std::vector< double > * tracketa
Definition: Run.h:43