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 
7 #include "FWCore/Framework/interface/Frameworkfwd.h"
8 #include "FWCore/Framework/interface/EDAnalyzer.h"
9 #include "FWCore/Framework/interface/Event.h"
10 #include "FWCore/Framework/interface/MakerMacros.h"
11 #include "FWCore/Framework/interface/ESHandle.h"
12 #include "FWCore/ParameterSet/interface/ParameterSet.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 #include "FWCore/Utilities/interface/Exception.h"
15 #include "FWCore/Utilities/interface/EDGetToken.h"
16 
17 #include "CondFormats/RunInfo/interface/RunInfo.h"
18 
19 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
20 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
21 #include "Geometry/CommonTopologies/interface/StripTopology.h"
22 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
23 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
24 
25 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
26 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
27 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
28 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
29 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
30 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
31 
32 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
33 
34 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
35 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
36 
37 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
38 #include "DataFormats/TrackReco/interface/Track.h"
39 #include "DataFormats/TrackReco/interface/TrackFwd.h"
40 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
41 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
42 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
43 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
44 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
45 #include "DataFormats/DetId/interface/DetId.h"
46 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
47 #include "DataFormats/TrackReco/interface/DeDxHit.h"
48 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
49 
50 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
51 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
52 
53 #include "TrackingTools/PatternTools/interface/Trajectory.h"
54 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
55 
56 #include "DQMServices/Core/interface/DQMStore.h"
57 #include "DQMServices/Core/interface/MonitorElement.h"
58 #include "FWCore/ServiceRegistry/interface/Service.h"
59 #include "CommonTools/UtilAlgos/interface/TFileService.h"
60 
61 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
62 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
63 
64 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
65 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
66 
67 
68 #include "TFile.h"
69 #include "TObjString.h"
70 #include "TString.h"
71 #include "TH1F.h"
72 #include "TH2S.h"
73 #include "TProfile.h"
74 #include "TF1.h"
75 #include "TROOT.h"
76 #include "TTree.h"
77 #include "TChain.h"
78 
79 #include <ext/hash_map>
80 
81 
82 
83 using namespace edm;
84 using namespace reco;
85 using namespace std;
86 using namespace __gnu_cxx;
87 using __gnu_cxx::hash_map;
88 using __gnu_cxx::hash;
89 
90 struct stAPVGain{
91  unsigned int Index;
92  int Bin;
93  unsigned int DetId;
94  unsigned int APVId;
95  unsigned int SubDet;
96  float x;
97  float y;
98  float z;
99  float Eta;
100  float R;
101  float Phi;
102  float Thickness;
103  double FitMPV;
104  double FitMPVErr;
105  double FitWidth;
106  double FitWidthErr;
107  double FitChi2;
108  double FitNorm;
109  double Gain;
110  double CalibGain;
111  double PreviousGain;
112  double PreviousGainTick;
113  double NEntries;
114  TH1F* HCharge;
115  TH1F* HChargeN;
116  bool isMasked;
117 };
118 
119 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
120 public:
121  explicit SiStripGainFromCalibTree(const edm::ParameterSet&);
122  ~SiStripGainFromCalibTree();
123 
124 
125 private:
126 
127 
128  virtual void algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
129  virtual void algoEndRun (const edm::Run& run, const edm::EventSetup& iSetup) override;
130  virtual void algoBeginJob (const edm::EventSetup& iSetup) override;
131  virtual void algoEndJob () override;
132  virtual void algoAnalyze (const edm::Event &, const edm::EventSetup &) override;
133 
134  int statCollectionFromMode(const char* tag) const;
135  void bookDQMHistos(const char* dqm_dir, const char* tag);
136 
137  bool isBFieldConsistentWithMode( const edm::EventSetup& iSetup) const;
138  void swapBFieldMode(void);
139 
140  void merge(TH2* A, TH2* B); //needed to add histograms with different number of bins
141  void algoAnalyzeTheTree();
142  void algoComputeMPVandGain();
143  void processEvent(); //what really does the job
144 
145  void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
146  bool IsGoodLandauFit(double* FitResults);
147  void storeOnTree(TFileService* tfs);
148  void MakeCalibrationMap();
149  bool produceTagFilter();
150 
151  template<typename T>
152  inline edm::Handle<T> connect(const T* &ptr, edm::EDGetTokenT<T> token, const edm::Event &evt) {
153  edm::Handle<T> handle;
154  evt.getByToken(token, handle);
155  ptr = handle.product();
156  return handle; //return handle to keep alive pointer (safety first)
157  }
158 
159  SiStripApvGain* getNewObject() override;
160 
161  TFileService *tfs;
162  DQMStore* dbe;
163  double MagFieldCurrentTh;
164  double MinNrEntries;
165  double MaxMPVError;
166  double MaxChi2OverNDF;
167  double MinTrackMomentum;
168  double MaxTrackMomentum;
169  double MinTrackEta;
170  double MaxTrackEta;
171  unsigned int MaxNrStrips;
172  unsigned int MinTrackHits;
173  double MaxTrackChiOverNdf;
174  int MaxTrackingIteration;
175  bool AllowSaturation;
176  bool FirstSetOfConstants;
177  bool Validation;
178  bool OldGainRemoving;
179  int CalibrationLevel;
180 
181  bool saveSummary;
182  bool useCalibration;
183  bool m_harvestingMode;
184  bool m_splitDQMstat;
185  string m_calibrationMode;
186  string m_calibrationPath;
187  string m_DQMdir;
188 
189  double tagCondition_NClusters;
190  double tagCondition_GoodFrac;
191 
192  std::string AlgoMode;
193  std::string OutputGains;
194  vector<string> VInputFiles;
195 
196  enum statistic_type {None=-1, StdBunch, StdBunch0T, FaABunch, FaABunch0T, IsoBunch, IsoBunch0T, Harvest};
197 
198  std::vector<string> dqm_tag_;
199 
200  std::vector<MonitorElement*> Charge_Vs_Index;
201  //std::vector<MonitorElement*> Charge_Vs_Index_Absolute;
202  std::vector<MonitorElement*> Charge_Vs_PathlengthTIB;
203  std::vector<MonitorElement*> Charge_Vs_PathlengthTOB;
204  std::vector<MonitorElement*> Charge_Vs_PathlengthTIDP;
205  std::vector<MonitorElement*> Charge_Vs_PathlengthTIDM;
206  std::vector<MonitorElement*> Charge_Vs_PathlengthTECP1;
207  std::vector<MonitorElement*> Charge_Vs_PathlengthTECP2;
208  std::vector<MonitorElement*> Charge_Vs_PathlengthTECM1;
209  std::vector<MonitorElement*> Charge_Vs_PathlengthTECM2;
210 
211  unsigned int NEvent;
212  unsigned int NTrack;
213  unsigned int NClusterStrip;
214  unsigned int NClusterPixel;
215  int NStripAPVs;
216  int NPixelDets;
217  unsigned int SRun;
218  unsigned int ERun;
219  unsigned int GOOD;
220  unsigned int BAD;
221  unsigned int MASKED;
222 
223  //Data members for processing
224 
225  //Event data
226  unsigned int eventnumber =0;
227  unsigned int runnumber =0;
228  const std::vector<bool>* TrigTech =0; edm::EDGetTokenT<std::vector<bool> > TrigTech_token_;
229 
230  // Track data
231  const std::vector<double>* trackchi2ndof =0; edm::EDGetTokenT<std::vector<double> > trackchi2ndof_token_;
232  const std::vector<float>* trackp =0; edm::EDGetTokenT<std::vector<float> > trackp_token_;
233  const std::vector<float>* trackpt =0; edm::EDGetTokenT<std::vector<float> > trackpt_token_;
234  const std::vector<double>* tracketa =0; edm::EDGetTokenT<std::vector<double> > tracketa_token_;
235  const std::vector<double>* trackphi =0; edm::EDGetTokenT<std::vector<double> > trackphi_token_;
236  const std::vector<unsigned int>* trackhitsvalid =0; edm::EDGetTokenT<std::vector<unsigned int> > trackhitsvalid_token_;
237  const std::vector<int>* trackalgo =0; edm::EDGetTokenT<std::vector<int> > trackalgo_token_;
238 
239  // CalibTree data
240  const std::vector<int>* trackindex =0; edm::EDGetTokenT<std::vector<int> > trackindex_token_;
241  const std::vector<unsigned int>* rawid =0; edm::EDGetTokenT<std::vector<unsigned int> > rawid_token_;
242  const std::vector<double>* localdirx =0; edm::EDGetTokenT<std::vector<double> > localdirx_token_;
243  const std::vector<double>* localdiry =0; edm::EDGetTokenT<std::vector<double> > localdiry_token_;
244  const std::vector<double>* localdirz =0; edm::EDGetTokenT<std::vector<double> > localdirz_token_;
245  const std::vector<unsigned short>* firststrip =0; edm::EDGetTokenT<std::vector<unsigned short> > firststrip_token_;
246  const std::vector<unsigned short>* nstrips =0; edm::EDGetTokenT<std::vector<unsigned short> > nstrips_token_;
247  const std::vector<bool>* saturation =0; edm::EDGetTokenT<std::vector<bool> > saturation_token_;
248  const std::vector<bool>* overlapping =0; edm::EDGetTokenT<std::vector<bool> > overlapping_token_;
249  const std::vector<bool>* farfromedge =0; edm::EDGetTokenT<std::vector<bool> > farfromedge_token_;
250  const std::vector<unsigned int>* charge =0; edm::EDGetTokenT<std::vector<unsigned int> > charge_token_;
251  const std::vector<double>* path =0; edm::EDGetTokenT<std::vector<double> > path_token_;
252  const std::vector<double>* chargeoverpath =0; edm::EDGetTokenT<std::vector<double> > chargeoverpath_token_;
253  const std::vector<unsigned char>* amplitude =0; edm::EDGetTokenT<std::vector<unsigned char> > amplitude_token_;
254  const std::vector<double>* gainused =0; edm::EDGetTokenT<std::vector<double> > gainused_token_;
255 
256  string EventPrefix_; //("");
257  string EventSuffix_; //("");
258  string TrackPrefix_; //("track");
259  string TrackSuffix_; //("");
260  string CalibPrefix_; //("GainCalibration");
261  string CalibSuffix_; //("");
262 
263 private :
264  class isEqual{
265  public:
266  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
267  };
268 
269  std::vector<stAPVGain*> APVsCollOrdered;
270  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
271 };
272 
273 inline int
274 SiStripGainFromCalibTree::statCollectionFromMode(const char* tag) const
275 {
276  std::vector<string>::const_iterator it=dqm_tag_.begin();
277  while(it!=dqm_tag_.end()) {
278  if(*it==std::string(tag)) return it-dqm_tag_.begin();
279  it++;
280  }
281 
282  if (std::string(tag)=="") return 0; // return StdBunch calibration mode for backward compatibility
283 
284  return None;
285 }
286 
287 void SiStripGainFromCalibTree::merge(TH2* A, TH2* B){
288  if(A->GetNbinsX() == B->GetNbinsX()){
289  A->Add(B);
290  }else{
291  for(int x=0;x<=B->GetNbinsX()+1; x++){
292  for(int y=0;y<=B->GetNbinsY()+1; y++){
293  A->SetBinContent(x,y,A->GetBinContent(x,y)+B->GetBinContent(x,y));
294  }}
295  }
296 }
297 
298 
299 
300 SiStripGainFromCalibTree::SiStripGainFromCalibTree(const edm::ParameterSet& iConfig) : ConditionDBWriter<SiStripApvGain>(iConfig)
301 {
302  OutputGains = iConfig.getParameter<std::string>("OutputGains");
303 
304  AlgoMode = iConfig.getUntrackedParameter<std::string>("AlgoMode", "CalibTree");
305  MagFieldCurrentTh = iConfig.getUntrackedParameter<double> ("MagFieldCurrentTh" , 2000.);
306  MinNrEntries = iConfig.getUntrackedParameter<double> ("minNrEntries" , 20);
307  MaxMPVError = iConfig.getUntrackedParameter<double> ("maxMPVError" , 500.0);
308  MaxChi2OverNDF = iConfig.getUntrackedParameter<double> ("maxChi2OverNDF" , 5.0);
309  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 3.0);
310  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
311  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
312  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
313  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 2);
314  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 8);
315  MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double> ("MaxTrackChiOverNdf" , 3);
316  MaxTrackingIteration = iConfig.getUntrackedParameter<int> ("MaxTrackingIteration" , 7);
317  AllowSaturation = iConfig.getUntrackedParameter<bool> ("AllowSaturation" , false);
318  FirstSetOfConstants = iConfig.getUntrackedParameter<bool> ("FirstSetOfConstants" , true);
319  Validation = iConfig.getUntrackedParameter<bool> ("Validation" , false);
320  OldGainRemoving = iConfig.getUntrackedParameter<bool> ("OldGainRemoving" , false);
321 
322  CalibrationLevel = iConfig.getUntrackedParameter<int> ("CalibrationLevel" , 0);
323  VInputFiles = iConfig.getUntrackedParameter<vector<string> > ("InputFiles");
324 
325 
326  useCalibration = iConfig.getUntrackedParameter<bool> ("UseCalibration" , false);
327  m_harvestingMode = iConfig.getUntrackedParameter<bool> ("harvestingMode" , false);
328  m_splitDQMstat = iConfig.getUntrackedParameter<bool> ("splitDQMstat" , false);
329  m_calibrationMode = iConfig.getUntrackedParameter<string> ("calibrationMode" , "StdBunch");
330  m_calibrationPath = iConfig.getUntrackedParameter<string> ("calibrationPath");
331  m_DQMdir = iConfig.getUntrackedParameter<string> ("DQMdir" , "AlCaReco/SiStripGains");
332 
333  tagCondition_NClusters = iConfig.getUntrackedParameter<double> ("NClustersForTagProd" , 2E8);
334  tagCondition_GoodFrac = iConfig.getUntrackedParameter<double> ("GoodFracForTagProd" , 0.95);
335 
336  saveSummary = iConfig.getUntrackedParameter<bool> ("saveSummary" , false);
337 
338 
339  // Gather DQM Service
340  dbe = edm::Service<DQMStore>().operator->();
341  dbe->setVerbose(10);
342 
343  //Set the monitoring element tag and store
344  dqm_tag_.reserve(7);
345  dqm_tag_.clear();
346  dqm_tag_.push_back( "StdBunch" ); // statistic collection from Standard Collision Bunch @ 3.8 T
347  dqm_tag_.push_back( "StdBunch0T" ); // statistic collection from Standard Collision Bunch @ 0 T
348  dqm_tag_.push_back( "AagBunch" ); // statistic collection from First Collision After Abort Gap @ 3.8 T
349  dqm_tag_.push_back( "AagBunch0T" ); // statistic collection from First Collision After Abort Gap @ 0 T
350  dqm_tag_.push_back( "IsoMuon" ); // statistic collection from Isolated Muon @ 3.8 T
351  dqm_tag_.push_back( "IsoMuon0T" ); // statistic collection from Isolated Muon @ 0 T
352  dqm_tag_.push_back( "Harvest" ); // statistic collection: Harvest
353 
354  Charge_Vs_Index.insert( Charge_Vs_Index.begin(), dqm_tag_.size(), 0);
355  //Charge_Vs_Index_Absolute.insert( Charge_Vs_Index_Absolute.begin(), dqm_tag_.size(), 0);
356  Charge_Vs_PathlengthTIB.insert( Charge_Vs_PathlengthTIB.begin(), dqm_tag_.size(), 0);
357  Charge_Vs_PathlengthTOB.insert( Charge_Vs_PathlengthTOB.begin(), dqm_tag_.size(), 0);
358  Charge_Vs_PathlengthTIDP.insert( Charge_Vs_PathlengthTIDP.begin(), dqm_tag_.size(), 0);
359  Charge_Vs_PathlengthTIDM.insert( Charge_Vs_PathlengthTIDM.begin(), dqm_tag_.size(), 0);
360  Charge_Vs_PathlengthTECP1.insert( Charge_Vs_PathlengthTECP1.begin(), dqm_tag_.size(), 0);
361  Charge_Vs_PathlengthTECP2.insert( Charge_Vs_PathlengthTECP2.begin(), dqm_tag_.size(), 0);
362  Charge_Vs_PathlengthTECM1.insert( Charge_Vs_PathlengthTECM1.begin(), dqm_tag_.size(), 0);
363  Charge_Vs_PathlengthTECM2.insert( Charge_Vs_PathlengthTECM2.begin(), dqm_tag_.size(), 0);
364 
365 
366 
367 
368  // configure token for gathering the ntuple variables
369  edm::ParameterSet swhallowgain_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("gain");
370 
371  string label = swhallowgain_pset.getUntrackedParameter<string>("label");
372  CalibPrefix_ = swhallowgain_pset.getUntrackedParameter<string>("prefix");
373  CalibSuffix_ = swhallowgain_pset.getUntrackedParameter<string>("suffix");
374 
375  trackindex_token_ = consumes<std::vector<int> >(edm::InputTag(label, CalibPrefix_ + "trackindex" + CalibSuffix_));
376  rawid_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, CalibPrefix_ + "rawid" + CalibSuffix_));
377  localdirx_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdirx" + CalibSuffix_));
378  localdiry_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdiry" + CalibSuffix_));
379  localdirz_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "localdirz" + CalibSuffix_));
380  firststrip_token_ = consumes<std::vector<unsigned short> >(edm::InputTag(label, CalibPrefix_ + "firststrip" + CalibSuffix_));
381  nstrips_token_ = consumes<std::vector<unsigned short> >(edm::InputTag(label, CalibPrefix_ + "nstrips" + CalibSuffix_));
382  saturation_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "saturation" + CalibSuffix_));
383  overlapping_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "overlapping" + CalibSuffix_));
384  farfromedge_token_ = consumes<std::vector<bool> >(edm::InputTag(label, CalibPrefix_ + "farfromedge" + CalibSuffix_));
385  charge_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, CalibPrefix_ + "charge" + CalibSuffix_));
386  path_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "path" + CalibSuffix_));
387  chargeoverpath_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "chargeoverpath"+ CalibSuffix_));
388  amplitude_token_ = consumes<std::vector<unsigned char> >(edm::InputTag(label, CalibPrefix_ + "amplitude" + CalibSuffix_));
389  gainused_token_ = consumes<std::vector<double> >(edm::InputTag(label, CalibPrefix_ + "gainused" + CalibSuffix_));
390 
391  edm::ParameterSet evtinfo_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("evtinfo");
392  label = evtinfo_pset.getUntrackedParameter<string>("label");
393  EventPrefix_ = evtinfo_pset.getUntrackedParameter<string>("prefix");
394  EventSuffix_ = evtinfo_pset.getUntrackedParameter<string>("suffix");
395  TrigTech_token_ = consumes<std::vector<bool> >(edm::InputTag(label, EventPrefix_ + "TrigTech" + EventSuffix_));
396 
397  edm::ParameterSet track_pset = iConfig.getUntrackedParameter<edm::ParameterSet>("tracks");
398  label = track_pset.getUntrackedParameter<string>("label");
399  TrackPrefix_ = track_pset.getUntrackedParameter<string>("prefix");
400  TrackSuffix_ = track_pset.getUntrackedParameter<string>("suffix");
401 
402  trackchi2ndof_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "chi2ndof" + TrackSuffix_));
403  trackp_token_ = consumes<std::vector<float> >(edm::InputTag(label, TrackPrefix_ + "momentum" + TrackSuffix_));
404  trackpt_token_ = consumes<std::vector<float> >(edm::InputTag(label, TrackPrefix_ + "pt" + TrackSuffix_));
405  tracketa_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "eta" + TrackSuffix_));
406  trackphi_token_ = consumes<std::vector<double> >(edm::InputTag(label, TrackPrefix_ + "phi" + TrackSuffix_));
407  trackhitsvalid_token_ = consumes<std::vector<unsigned int> >(edm::InputTag(label, TrackPrefix_ + "hitsvalid" + TrackSuffix_));
408  trackalgo_token_ = consumes<std::vector<int> >(edm::InputTag(label, TrackPrefix_ + "algo" + TrackSuffix_));
409 }
410 
411 void SiStripGainFromCalibTree::bookDQMHistos(const char* dqm_dir, const char* tag)
412 {
413  static std::string booked_dir = "";
414 
415  edm::LogInfo("SiStripGainFromCalibTree") << "Setting " << dqm_dir << "in DQM and booking histograms for tag "
416  << tag << std::endl;
417 
418  if ( strcmp(booked_dir.c_str(),dqm_dir)!=0 ) {
419  booked_dir = dqm_dir;
420  dbe->setCurrentFolder(dqm_dir);
421  }
422 
423  std::string stag(tag);
424  if(stag.size()!=0 && stag[0]!='_') stag.insert(0,1,'_');
425 
426  std::string cvi = std::string("Charge_Vs_Index") + stag;
427  //std::string cviA = std::string("Charge_Vs_Index_Absolute") + stag;
428  std::string cvpTIB = std::string("Charge_Vs_PathlengthTIB") + stag;
429  std::string cvpTOB = std::string("Charge_Vs_PathlengthTOB") + stag;
430  std::string cvpTIDP = std::string("Charge_Vs_PathlengthTIDP") + stag;
431  std::string cvpTIDM = std::string("Charge_Vs_PathlengthTIDM") + stag;
432  std::string cvpTECP1 = std::string("Charge_Vs_PathlengthTECP1") + stag;
433  std::string cvpTECP2 = std::string("Charge_Vs_PathlengthTECP2") + stag;
434  std::string cvpTECM1 = std::string("Charge_Vs_PathlengthTECM1") + stag;
435  std::string cvpTECM2 = std::string("Charge_Vs_PathlengthTECM2") + stag;
436 
437  int elepos = (m_harvestingMode && AlgoMode=="PCL")? Harvest : statCollectionFromMode(tag);
438 
439  Charge_Vs_Index[elepos] = dbe->book2S(cvi.c_str() , cvi.c_str() , 88625, 0 , 88624,2000,0,4000);
440  //Charge_Vs_Index_Absolute[elepos] = dbe->book2S(cviA.c_str() , cviA.c_str() , 88625, 0 , 88624,1000,0,4000);
441  Charge_Vs_PathlengthTIB[elepos] = dbe->book2S(cvpTIB.c_str() , cvpTIB.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
442  Charge_Vs_PathlengthTOB[elepos] = dbe->book2S(cvpTOB.c_str() , cvpTOB.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
443  Charge_Vs_PathlengthTIDP[elepos] = dbe->book2S(cvpTIDP.c_str() , cvpTIDP.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
444  Charge_Vs_PathlengthTIDM[elepos] = dbe->book2S(cvpTIDM.c_str() , cvpTIDM.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
445  Charge_Vs_PathlengthTECP1[elepos] = dbe->book2S(cvpTECP1.c_str(), cvpTECP1.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
446  Charge_Vs_PathlengthTECP2[elepos] = dbe->book2S(cvpTECP2.c_str(), cvpTECP2.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
447  Charge_Vs_PathlengthTECM1[elepos] = dbe->book2S(cvpTECM1.c_str(), cvpTECM1.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
448  Charge_Vs_PathlengthTECM2[elepos] = dbe->book2S(cvpTECM2.c_str(), cvpTECM2.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
449 }
450 
451 void SiStripGainFromCalibTree::algoBeginJob(const edm::EventSetup& iSetup)
452 {
453  edm::LogInfo("SiStripGainFromCalibTree") << "AlgoMode : " << AlgoMode << "\n"
454  << "CalibrationMode : " << m_calibrationMode << "\n"
455  << "HarvestingMode : " << m_harvestingMode << std::endl;
456  //Setup DQM histograms
457  if(AlgoMode != "PCL" or m_harvestingMode) {
458  const char * dqm_dir = "AlCaReco/SiStripGainsHarvesting/";
459  this->bookDQMHistos( dqm_dir, dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str() );
460  } else {
461  //Check consistency of calibration Mode and BField only for the ALCAPROMPT in the PCL workflow
462  if (!isBFieldConsistentWithMode(iSetup)) {
463  string prevMode = m_calibrationMode;
464  swapBFieldMode();
465  edm::LogInfo("SiStripGainFromCalibTree") << "Switching calibration mode for endorsing BField status: "
466  << prevMode << " ==> " << m_calibrationMode << std::endl;
467  }
468  std::string dqm_dir = m_DQMdir + ((m_splitDQMstat)? m_calibrationMode:"") + "/";
469  int elem = statCollectionFromMode(m_calibrationMode.c_str());
470  this->bookDQMHistos( dqm_dir.c_str(), dqm_tag_[elem].c_str() );
471  }
472 
473 
474  edm::ESHandle<TrackerGeometry> tkGeom;
475  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
476  auto const & Det = tkGeom->dets();
477 
478  NPixelDets = 0;
479  NStripAPVs = 0;
480  unsigned int Index=0;
481  for(unsigned int i=0;i<Det.size();i++){
482  DetId Detid = Det[i]->geographicalId();
483  int SubDet = Detid.subdetId();
484 
485  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
486  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
487 
488  auto DetUnit = dynamic_cast<const StripGeomDetUnit*> (Det[i]);
489  if(!DetUnit)continue;
490 
491  const StripTopology& Topo = DetUnit->specificTopology();
492  unsigned int NAPV = Topo.nstrips()/128;
493 
494  for(unsigned int j=0;j<NAPV;j++){
495  stAPVGain* APV = new stAPVGain;
496  APV->Index = Index;
497  APV->Bin = -1;
498  APV->DetId = Detid.rawId();
499  APV->APVId = j;
500  APV->SubDet = SubDet;
501  APV->FitMPV = -1;
502  APV->FitMPVErr = -1;
503  APV->FitWidth = -1;
504  APV->FitWidthErr = -1;
505  APV->FitChi2 = -1;
506  APV->FitNorm = -1;
507  APV->Gain = -1;
508  APV->PreviousGain = 1;
509  APV->PreviousGainTick = 1;
510  APV->x = DetUnit->position().basicVector().x();
511  APV->y = DetUnit->position().basicVector().y();
512  APV->z = DetUnit->position().basicVector().z();
513  APV->Eta = DetUnit->position().basicVector().eta();
514  APV->Phi = DetUnit->position().basicVector().phi();
515  APV->R = DetUnit->position().basicVector().transverse();
516  APV->Thickness = DetUnit->surface().bounds().thickness();
517  APV->NEntries = 0;
518  APV->isMasked = false;
519 
520  APVsCollOrdered.push_back(APV);
521  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
522  Index++;
523  NStripAPVs++;
524  }
525  }
526  }
527 
528  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
529  DetId Detid = Det[i]->geographicalId();
530  int SubDet = Detid.subdetId();
531  if( SubDet == PixelSubdetector::PixelBarrel || PixelSubdetector::PixelEndcap ){
532  auto DetUnit = dynamic_cast<const PixelGeomDetUnit*> (Det[i]);
533  if(!DetUnit) continue;
534 
535  const PixelTopology& Topo = DetUnit->specificTopology();
536  unsigned int NROCRow = Topo.nrows()/(80.);
537  unsigned int NROCCol = Topo.ncolumns()/(52.);
538 
539  for(unsigned int j=0;j<NROCRow;j++){
540  for(unsigned int i=0;i<NROCCol;i++){
541 
542  stAPVGain* APV = new stAPVGain;
543  APV->Index = Index;
544  APV->Bin = -1;
545  APV->DetId = Detid.rawId();
546  APV->APVId = (j<<3 | i);
547  APV->SubDet = SubDet;
548  APV->FitMPV = -1;
549  APV->FitMPVErr = -1;
550  APV->FitWidth = -1;
551  APV->FitWidthErr = -1;
552  APV->FitChi2 = -1;
553  APV->Gain = -1;
554  APV->PreviousGain = 1;
555  APV->x = DetUnit->position().basicVector().x();
556  APV->y = DetUnit->position().basicVector().y();
557  APV->z = DetUnit->position().basicVector().z();
558  APV->Eta = DetUnit->position().basicVector().eta();
559  APV->Phi = DetUnit->position().basicVector().phi();
560  APV->R = DetUnit->position().basicVector().transverse();
561  APV->Thickness = DetUnit->surface().bounds().thickness();
562  APV->isMasked = false; //SiPixelQuality_->IsModuleBad(Detid.rawId());
563  APV->NEntries = 0;
564 
565  APVsCollOrdered.push_back(APV);
566  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
567  Index++;
568  NPixelDets++;
569  }}
570  }
571  }
572 
573 
574  MakeCalibrationMap();
575 
576  NEvent = 0;
577  NTrack = 0;
578  NClusterStrip = 0;
579  NClusterPixel = 0;
580  SRun = 1<<31;
581  ERun = 0;
582  GOOD = 0;
583  BAD = 0;
584  MASKED = 0;
585 }
586 
587 bool SiStripGainFromCalibTree::isBFieldConsistentWithMode( const edm::EventSetup& iSetup) const {
588  edm::ESHandle<RunInfo> runInfo;
589  iSetup.get<RunInfoRcd>().get(runInfo);
590 
591  double average_current = runInfo.product()->m_avg_current;
592  bool isOn = (average_current > MagFieldCurrentTh);
593  bool is0T = (m_calibrationMode.substr( m_calibrationMode.length()-2 )=="0T");
594 
595  return ( (isOn && !is0T) || (!isOn && is0T) );
596 }
597 
598 void SiStripGainFromCalibTree::swapBFieldMode() {
599  if ( m_calibrationMode.substr( m_calibrationMode.length()-2 )=="0T" ) {
600  m_calibrationMode.erase( m_calibrationMode.length()-2,2);
601  } else {
602  m_calibrationMode.append( "0T" );
603  }
604 }
605 
606 void SiStripGainFromCalibTree::algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup)
607 {
608  if( !m_harvestingMode && AlgoMode=="PCL") {
609  //Check consistency of calibration Mode and BField only for the ALCAPROMPT in the PCL workflow
610  if (!isBFieldConsistentWithMode(iSetup)) {
611  string prevMode = m_calibrationMode;
612  swapBFieldMode();
613  edm::LogInfo("SiStripGainFromCalibTree") << "Switching calibration mode for endorsing BField status: "
614  << prevMode << " ==> " << m_calibrationMode << std::endl;
615  }
616  }
617 
618  edm::ESHandle<SiStripGain> gainHandle;
619  iSetup.get<SiStripGainRcd>().get(gainHandle);
620  if(!gainHandle.isValid()){edm::LogError("SiStripGainFromCalibTree")<< "gainHandle is not valid\n"; exit(0);}
621 
622  edm::ESHandle<SiStripQuality> SiStripQuality_;
623  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
624 
625  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
626  stAPVGain* APV = APVsCollOrdered[a];
627 
628  APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId,APV->APVId);
629 // if(!FirstSetOfConstants){
630  if(gainHandle->getNumberOfTags()!=2){edm::LogError("SiStripGainFromCalibTree")<< "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";fflush(stdout);exit(0);};
631  float newPreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
632  if(APV->PreviousGain!=1 and newPreviousGain!=APV->PreviousGain)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: ParticleGain in the global tag changed\n";
633  APV->PreviousGain = newPreviousGain;
634 
635  float newPreviousGainTick = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 0),1);
636  if(APV->PreviousGainTick!=1 and newPreviousGainTick!=APV->PreviousGainTick)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: TickMarkGain in the global tag changed\n";
637  APV->PreviousGainTick = newPreviousGainTick;
638 
639 
640  //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f (G1) x %8.4f (G2)\n",APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain);
641 
642 
643 
644 // }
645  }
646 }
647 
648 void SiStripGainFromCalibTree::algoEndRun(const edm::Run& run, const edm::EventSetup& iSetup){
649  if(AlgoMode == "PCL" && !m_harvestingMode) return;//nothing to do in that case
650 
651  if(AlgoMode == "PCL" and m_harvestingMode){
652  // Load the 2D histograms from the DQM objects
653  // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
654  // DQMStore so that they can be used in the fit
655 
656  edm::LogInfo("SiStripGainFromCalibTree") << "Starting harvesting statistics" << std::endl;
657 
658  // check the required tag before adding histograms
659  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
660  if (elepos!=Harvest) {
661 
662  //collect statistics from DQM into the related monitored elements
663  std::string stag = m_calibrationMode;
664  if(stag.size()!=0 && stag[0]!='_') stag.insert(0,1,'_');
665 
666 
667  if (elepos==-1) {
668  //implememt backward compatibility
669  elepos=0;
670  stag="";
671  }
672 
673  std::string DQM_dir = m_DQMdir + ((m_splitDQMstat)? m_calibrationMode:"");
674 
675  std::string cvi = DQM_dir + std::string("/Charge_Vs_Index") + stag;
676  //std::string cviA = DQM_dir + std::string("/Charge_Vs_Index_Absolute") + stag;
677  std::string cvpTIB = DQM_dir + std::string("/Charge_Vs_PathlengthTIB") + stag;
678  std::string cvpTOB = DQM_dir + std::string("/Charge_Vs_PathlengthTOB") + stag;
679  std::string cvpTIDP = DQM_dir + std::string("/Charge_Vs_PathlengthTIDP") + stag;
680  std::string cvpTIDM = DQM_dir + std::string("/Charge_Vs_PathlengthTIDM") + stag;
681  std::string cvpTECP1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP1") + stag;
682  std::string cvpTECP2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP2") + stag;
683  std::string cvpTECM1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM1") + stag;
684  std::string cvpTECM2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM2") + stag;
685 
686  Charge_Vs_Index[elepos] = dbe->get(cvi.c_str());
687  //Charge_Vs_Index_Absolute[elepos] = dbe->get(cviA.c_str());
688  Charge_Vs_PathlengthTIB[elepos] = dbe->get(cvpTIB.c_str());
689  Charge_Vs_PathlengthTOB[elepos] = dbe->get(cvpTOB.c_str());
690  Charge_Vs_PathlengthTIDP[elepos] = dbe->get(cvpTIDP.c_str());
691  Charge_Vs_PathlengthTIDM[elepos] = dbe->get(cvpTIDM.c_str());
692  Charge_Vs_PathlengthTECP1[elepos] = dbe->get(cvpTECP1.c_str());
693  Charge_Vs_PathlengthTECP2[elepos] = dbe->get(cvpTECP2.c_str());
694  Charge_Vs_PathlengthTECM1[elepos] = dbe->get(cvpTECM1.c_str());
695  Charge_Vs_PathlengthTECM2[elepos] = dbe->get(cvpTECM2.c_str());
696 
697  if (Charge_Vs_Index[elepos]==0) {
698  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvi.c_str()
699  << ", statistics will not be summed!" << std::endl;
700  } else {
701  merge( (Charge_Vs_Index[Harvest])->getTH2S(), (Charge_Vs_Index[elepos])->getTH2S() );
702  edm::LogInfo("SiStripGainFromCalibTree") << "Harvesting "
703  << (Charge_Vs_Index[elepos])->getTH2S()->GetEntries() << " more clusters" << std::endl;
704  }
705 
706  //if (Charge_Vs_Index_Absolute[elepos]==0) {
707  // edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cviA.c_str()
708  // << ", statistics will not be summed!" << std::endl;
709  //} else merge( (Charge_Vs_Index_Absolute[Harvest])->getTH2S(), (Charge_Vs_Index_Absolute[elepos])->getTH2S() );
710 
711  if (Charge_Vs_PathlengthTIB[elepos]==0) {
712  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIB.c_str()
713  << ", statistics will not be summed!" << std::endl;
714  } else (Charge_Vs_PathlengthTIB[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIB[elepos])->getTH2S());
715 
716  if (Charge_Vs_PathlengthTOB[elepos]==0) {
717  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTOB.c_str()
718  << ", statistics will not be summed!" << std::endl;
719  } else (Charge_Vs_PathlengthTOB[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTOB[elepos])->getTH2S());
720 
721  if (Charge_Vs_PathlengthTIDP[elepos]==0) {
722  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIDP.c_str()
723  << ", statistics will not be summed!" << std::endl;
724  } else (Charge_Vs_PathlengthTIDP[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIDP[elepos])->getTH2S());
725 
726  if (Charge_Vs_PathlengthTIDM[elepos]==0) {
727  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIDM.c_str()
728  << ", statistics will not be summed!" << std::endl;
729  } else (Charge_Vs_PathlengthTIDM[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTIDM[elepos])->getTH2S());
730 
731  if (Charge_Vs_PathlengthTECP1[elepos]==0) {
732  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECP1.c_str()
733  << ", statistics will not be summed!" << std::endl;
734  } else (Charge_Vs_PathlengthTECP1[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECP1[elepos])->getTH2S());
735 
736  if (Charge_Vs_PathlengthTECP2[elepos]==0) {
737  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECP2.c_str()
738  << ", statistics will not be summed!" << std::endl;
739  } else (Charge_Vs_PathlengthTECP2[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECP2[elepos])->getTH2S());
740 
741  if (Charge_Vs_PathlengthTECM1[elepos]==0) {
742  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECM1.c_str()
743  << ", statistics will not be summed!" << std::endl;
744  } else (Charge_Vs_PathlengthTECM1[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECM1[elepos])->getTH2S());
745 
746  if (Charge_Vs_PathlengthTECM2[elepos]==0) {
747  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECM2.c_str()
748  << ", statistics will not be summed!" << std::endl;
749  } else (Charge_Vs_PathlengthTECM2[Harvest])->getTH2S()->Add((Charge_Vs_PathlengthTECM2[elepos])->getTH2S());
750  }
751  }
752 }
753 
754 void
755 SiStripGainFromCalibTree::algoEndJob() {
756  if(AlgoMode == "PCL" && !m_harvestingMode) return;//nothing to do in that case
757 
758  if(AlgoMode == "CalibTree"){
759  edm::LogInfo("SiStripGainFromCalibTree") << "Analyzing calibration tree" << std::endl;
760  // Loop on calibTrees to fill the 2D histograms
761  algoAnalyzeTheTree();
762  }else if(m_harvestingMode){
763  NClusterStrip = (Charge_Vs_Index[Harvest])->getTH2S()->Integral(0,NStripAPVs+1, 0, 99999 );
764  NClusterPixel = (Charge_Vs_Index[Harvest])->getTH2S()->Integral(NStripAPVs+2, NStripAPVs+NPixelDets+2, 0, 99999 );
765  }
766 
767  // Now that we have the full statistics we can extract the information of the 2D histograms
768  algoComputeMPVandGain();
769 
770  // Force the DB object writing,
771  // thus setting the IOV as the first processed run (if timeFromEndRun is set to false)
772  storeOnDbNow();
773 
774  if(AlgoMode != "PCL" or saveSummary){
775  edm::LogInfo("SiStripGainFromCalibTree") << "Saving summary into root file" << std::endl;
776 
777  //also save the 2D monitor elements to this file as TH2D tfs
778  tfs = edm::Service<TFileService>().operator->();
779 
780  //save only the statistics for the calibrationTag
781  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
782 
783  if( Charge_Vs_Index[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_Index[elepos])->getTH2S() );
784  //if( Charge_Vs_Index_Absolute[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_Index_Absolute[elepos])->getTH2S() );
785  if( Charge_Vs_PathlengthTIB[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTIB[elepos])->getTH2S() );
786  if( Charge_Vs_PathlengthTOB[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTOB[elepos])->getTH2S() );
787  if( Charge_Vs_PathlengthTIDP[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTIDP[elepos])->getTH2S() );
788  if( Charge_Vs_PathlengthTIDM[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTIDM[elepos])->getTH2S() );
789  if( Charge_Vs_PathlengthTECP1[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTECP1[elepos])->getTH2S() );
790  if( Charge_Vs_PathlengthTECP2[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTECP2[elepos])->getTH2S() );
791  if( Charge_Vs_PathlengthTECM1[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTECM1[elepos])->getTH2S() );
792  if( Charge_Vs_PathlengthTECM2[elepos]!=0 ) tfs->make<TH2S> ( *(Charge_Vs_PathlengthTECM2[elepos])->getTH2S() );
793 
794  storeOnTree(tfs);
795  }
796 }
797 
798 
799 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
800 {
801  FitResults[0] = -0.5; //MPV
802  FitResults[1] = 0; //MPV error
803  FitResults[2] = -0.5; //Width
804  FitResults[3] = 0; //Width error
805  FitResults[4] = -0.5; //Fit Chi2/NDF
806  FitResults[5] = 0; //Normalization
807 
808  if( InputHisto->GetEntries() < MinNrEntries)return;
809 
810  // perform fit with standard landau
811  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
812  MyLandau->SetParameter(1,300);
813  InputHisto->Fit(MyLandau,"0QR WW");
814 
815  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
816  FitResults[0] = MyLandau->GetParameter(1); //MPV
817  FitResults[1] = MyLandau->GetParError(1); //MPV error
818  FitResults[2] = MyLandau->GetParameter(2); //Width
819  FitResults[3] = MyLandau->GetParError(2); //Width error
820  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
821  FitResults[5] = MyLandau->GetParameter(0);
822 
823  delete MyLandau;
824 }
825 
826 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
827  if(FitResults[0] <= 0 )return false;
828 // if(FitResults[1] > MaxMPVError )return false;
829 // if(FitResults[4] > MaxChi2OverNDF)return false;
830  return true;
831 }
832 
833 void SiStripGainFromCalibTree::processEvent() {
834 
835  edm::LogInfo("SiStripGainFromCalibTree") << "Processing run " << runnumber
836  << " and event " << eventnumber
837  << " for " << m_calibrationMode << " calibration." << std::endl;
838 
839  if(runnumber<SRun)SRun=runnumber;
840  if(runnumber>ERun)ERun=runnumber;
841 
842  NEvent++;
843  NTrack+=(*trackp).size();
844 
845  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
846 
847  unsigned int FirstAmplitude=0;
848  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
849  FirstAmplitude+=(*nstrips)[i];
850  int TI = (*trackindex)[i];
851 
852  //printf("%i - %i - %i %i %i\n", (int)(*rawid)[i], (int)(*firststrip)[i]/128, (int)(*farfromedge)[i], (int)(*overlapping)[i], (int)(*saturation )[i] );
853  if((*tracketa )[TI] < MinTrackEta )continue;
854  if((*tracketa )[TI] > MaxTrackEta )continue;
855  if((*trackp )[TI] < MinTrackMomentum )continue;
856  if((*trackp )[TI] > MaxTrackMomentum )continue;
857  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
858  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
859  if((*trackalgo )[TI] > MaxTrackingIteration )continue;
860 
861  stAPVGain* APV = APVsColl[((*rawid)[i]<<4) | ((*firststrip)[i]/128)]; //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
862 
863  if(APV->SubDet>2 && (*farfromedge)[i] == false )continue;
864  if(APV->SubDet>2 && (*overlapping)[i] == true )continue;
865  if(APV->SubDet>2 && (*saturation )[i] && !AllowSaturation)continue;
866  if(APV->SubDet>2 && (*nstrips )[i] > MaxNrStrips )continue;
867 
868 
869  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
870 
871  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
872  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
873  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
874 
875  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
876  //printf("Amplitudes: ");
877  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
878  //printf("\n");
879 
880  if(APV->SubDet>2){NClusterStrip++;}else{NClusterPixel++;}
881 
882  int Charge = 0;
883  if(APV->SubDet>2 && (useCalibration || !FirstSetOfConstants)){
884  bool Saturation = false;
885  for(unsigned int s=0;s<(*nstrips)[i];s++){
886  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
887  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
888  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
889  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
890  if(StripCharge>1024){
891  StripCharge = 255;
892  Saturation = true;
893  }else if(StripCharge>254){
894  StripCharge = 254;
895  Saturation = true;
896  }
897  Charge += StripCharge;
898  }
899  if(Saturation && !AllowSaturation)continue;
900  }else if(APV->SubDet>2){
901  Charge = (*charge)[i];
902  }else{
903  Charge = (*charge)[i]/265.0; //expected scale factor between pixel and strip charge
904  }
905 
906  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
907 
908  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
909  if(APV->SubDet>2){
910  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
911  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
912  }
913  //(Charge_Vs_Index_Absolute[elepos])->Fill(APV->Index,Charge);
914  (Charge_Vs_Index[elepos]) ->Fill(APV->Index,ClusterChargeOverPath);
915 
916  if(APV->SubDet==StripSubdetector::TIB){ (Charge_Vs_PathlengthTIB[elepos]) ->Fill((*path)[i],Charge);
917  }else if(APV->SubDet==StripSubdetector::TOB){ (Charge_Vs_PathlengthTOB[elepos]) ->Fill((*path)[i],Charge);
918  }else if(APV->SubDet==StripSubdetector::TID){
919  if(APV->Eta<0){ (Charge_Vs_PathlengthTIDM[elepos]) ->Fill((*path)[i],Charge);
920  }else if(APV->Eta>0){ (Charge_Vs_PathlengthTIDP[elepos]) ->Fill((*path)[i],Charge);
921  }
922  }else if(APV->SubDet==StripSubdetector::TEC){
923  if(APV->Eta<0){
924  if(APV->Thickness<0.04){ (Charge_Vs_PathlengthTECM1[elepos])->Fill((*path)[i],Charge);
925  }else if(APV->Thickness>0.04){ (Charge_Vs_PathlengthTECM2[elepos])->Fill((*path)[i],Charge);
926  }
927  }else if(APV->Eta>0){
928  if(APV->Thickness<0.04){ (Charge_Vs_PathlengthTECP1[elepos])->Fill((*path)[i],Charge);
929  }else if(APV->Thickness>0.04){ (Charge_Vs_PathlengthTECP2[elepos])->Fill((*path)[i],Charge);
930  }
931  }
932  }
933 
934  }// END OF ON-CLUSTER LOOP
935 }//END OF processEvent()
936 
937 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
938 {
939  for(unsigned int i=0;i<VInputFiles.size();i++){
940  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
941  TFile *tfile = TFile::Open(VInputFiles[i].c_str());
942  TString tree_path = TString::Format("gainCalibrationTree%s/tree",m_calibrationMode.c_str());
943  TTree* tree = dynamic_cast<TTree*> (tfile->Get(tree_path.Data()));
944 
945  tree->SetBranchAddress((EventPrefix_ + "event" + EventSuffix_).c_str(), &eventnumber , NULL);
946  tree->SetBranchAddress((EventPrefix_ + "run" + EventSuffix_).c_str(), &runnumber , NULL);
947  tree->SetBranchAddress((EventPrefix_ + "TrigTech" + EventSuffix_).c_str(), &TrigTech , NULL);
948 
949  tree->SetBranchAddress((TrackPrefix_ + "chi2ndof" + TrackSuffix_).c_str(), &trackchi2ndof , NULL);
950  tree->SetBranchAddress((TrackPrefix_ + "momentum" + TrackSuffix_).c_str(), &trackp , NULL);
951  tree->SetBranchAddress((TrackPrefix_ + "pt" + TrackSuffix_).c_str(), &trackpt , NULL);
952  tree->SetBranchAddress((TrackPrefix_ + "eta" + TrackSuffix_).c_str(), &tracketa , NULL);
953  tree->SetBranchAddress((TrackPrefix_ + "phi" + TrackSuffix_).c_str(), &trackphi , NULL);
954  tree->SetBranchAddress((TrackPrefix_ + "hitsvalid" + TrackSuffix_).c_str(), &trackhitsvalid, NULL);
955  tree->SetBranchAddress((TrackPrefix_ + "algo" + TrackSuffix_).c_str(), &trackalgo , NULL);
956 
957  tree->SetBranchAddress((CalibPrefix_ + "trackindex" + CalibSuffix_).c_str(), &trackindex , NULL);
958  tree->SetBranchAddress((CalibPrefix_ + "rawid" + CalibSuffix_).c_str(), &rawid , NULL);
959  tree->SetBranchAddress((CalibPrefix_ + "localdirx" + CalibSuffix_).c_str(), &localdirx , NULL);
960  tree->SetBranchAddress((CalibPrefix_ + "localdiry" + CalibSuffix_).c_str(), &localdiry , NULL);
961  tree->SetBranchAddress((CalibPrefix_ + "localdirz" + CalibSuffix_).c_str(), &localdirz , NULL);
962  tree->SetBranchAddress((CalibPrefix_ + "firststrip" + CalibSuffix_).c_str(), &firststrip , NULL);
963  tree->SetBranchAddress((CalibPrefix_ + "nstrips" + CalibSuffix_).c_str(), &nstrips , NULL);
964  tree->SetBranchAddress((CalibPrefix_ + "saturation" + CalibSuffix_).c_str(), &saturation , NULL);
965  tree->SetBranchAddress((CalibPrefix_ + "overlapping" + CalibSuffix_).c_str(), &overlapping , NULL);
966  tree->SetBranchAddress((CalibPrefix_ + "farfromedge" + CalibSuffix_).c_str(), &farfromedge , NULL);
967  tree->SetBranchAddress((CalibPrefix_ + "charge" + CalibSuffix_).c_str(), &charge , NULL);
968  tree->SetBranchAddress((CalibPrefix_ + "path" + CalibSuffix_).c_str(), &path , NULL);
969  tree->SetBranchAddress((CalibPrefix_ + "chargeoverpath" + CalibSuffix_).c_str(), &chargeoverpath, NULL);
970  tree->SetBranchAddress((CalibPrefix_ + "amplitude" + CalibSuffix_).c_str(), &amplitude , NULL);
971  tree->SetBranchAddress((CalibPrefix_ + "gainused" + CalibSuffix_).c_str(), &gainused , NULL);
972 
973 
974  unsigned int nentries = tree->GetEntries();
975  printf("Number of Events = %i + %i = %i\n",NEvent,nentries,(NEvent+nentries));
976  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
977  printf("Looping on the Tree :");
978  int TreeStep = nentries/50;if(TreeStep<=1)TreeStep=1;
979  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
980  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
981  tree->GetEntry(ientry);
982  processEvent();
983  }printf("\n");// END OF EVENT LOOP
984  }
985 }
986 
987 
988 
989 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
990  unsigned int I=0;
991  TH1F* Proj = NULL;
992  double FitResults[6];
993  double MPVmean = 300;
994 
995  int elepos = (AlgoMode == "PCL")? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
996 
997  if ( Charge_Vs_Index[elepos]==0 ) {
998  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not execute algoComputeMPVandGain method because "
999  << m_calibrationMode.c_str() << " statistics cannot be retrieved.\n"
1000  << "Please check if input contains "
1001  << m_calibrationMode.c_str() << " data." << std::endl;
1002  return;
1003  }
1004 
1005  TH2S *chvsidx = (Charge_Vs_Index[elepos])->getTH2S();
1006 
1007 
1008  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
1009  printf("Fitting Charge Distribution :");
1010  int TreeStep = APVsColl.size()/50;
1011  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
1012  if(I%TreeStep==0){printf(".");fflush(stdout);}
1013  stAPVGain* APV = it->second;
1014  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
1015 
1016  if(APV->isMasked){APV->Gain=APV->PreviousGain; MASKED++; continue;}
1017 
1018  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
1019  if(!Proj)continue;
1020 
1021  if(CalibrationLevel==0){
1022  }else if(CalibrationLevel==1){
1023  int SecondAPVId = APV->APVId;
1024  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
1025  stAPVGain* APV2 = APVsColl[(APV->DetId<<4) | SecondAPVId];
1026  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1027  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
1028  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
1029  }else if(CalibrationLevel==2){
1030  for(unsigned int i=0;i<16;i++){ //loop up to 6APV for Strip and up to 16 for Pixels
1031  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
1032  tmpit = APVsColl.find((APV->DetId<<4) | i);
1033  if(tmpit==APVsColl.end())continue;
1034  stAPVGain* APV2 = tmpit->second;
1035  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
1036  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1037  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
1038  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
1039  }
1040  }else{
1041  CalibrationLevel = 0;
1042  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
1043  }
1044 
1045  getPeakOfLandau(Proj,FitResults);
1046  APV->FitMPV = FitResults[0];
1047  APV->FitMPVErr = FitResults[1];
1048  APV->FitWidth = FitResults[2];
1049  APV->FitWidthErr = FitResults[3];
1050  APV->FitChi2 = FitResults[4];
1051  APV->FitNorm = FitResults[5];
1052  APV->NEntries = Proj->GetEntries();
1053 
1054  if(IsGoodLandauFit(FitResults)){
1055  APV->Gain = APV->FitMPV / MPVmean;
1056  if(APV->SubDet>2)GOOD++;
1057  }else{
1058  APV->Gain = APV->PreviousGain;
1059  if(APV->SubDet>2)BAD++;
1060  }
1061  if(APV->Gain<=0) APV->Gain = 1;
1062 
1063  //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);
1064  delete Proj;
1065  }printf("\n");
1066 }
1067 
1068 
1069 void SiStripGainFromCalibTree::storeOnTree(TFileService* tfs)
1070 {
1071  unsigned int tree_Index;
1072  unsigned int tree_Bin;
1073  unsigned int tree_DetId;
1074  unsigned char tree_APVId;
1075  unsigned char tree_SubDet;
1076  float tree_x;
1077  float tree_y;
1078  float tree_z;
1079  float tree_Eta;
1080  float tree_R;
1081  float tree_Phi;
1082  float tree_Thickness;
1083  float tree_FitMPV;
1084  float tree_FitMPVErr;
1085  float tree_FitWidth;
1086  float tree_FitWidthErr;
1087  float tree_FitChi2NDF;
1088  float tree_FitNorm;
1089  double tree_Gain;
1090  double tree_PrevGain;
1091  double tree_PrevGainTick;
1092  double tree_NEntries;
1093  bool tree_isMasked;
1094 
1095  TTree* MyTree;
1096  MyTree = tfs->make<TTree> ("APVGain","APVGain");
1097  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
1098  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
1099  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
1100  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
1101  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
1102  MyTree->Branch("x" ,&tree_x ,"x/F");
1103  MyTree->Branch("y" ,&tree_y ,"y/F");
1104  MyTree->Branch("z" ,&tree_z ,"z/F");
1105  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
1106  MyTree->Branch("R" ,&tree_R ,"R/F");
1107  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
1108  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
1109  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
1110  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
1111  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
1112  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
1113  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
1114  MyTree->Branch("FitNorm" ,&tree_FitNorm ,"FitNorm/F");
1115  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
1116  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
1117  MyTree->Branch("PrevGainTick" ,&tree_PrevGainTick,"PrevGainTick/D");
1118  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
1119  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
1120 
1121 
1122  FILE* Gains = stdout;
1123  fprintf(Gains,"NEvents = %i\n",NEvent);
1124  fprintf(Gains,"NTracks = %i\n",NTrack);
1125  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1126  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
1127  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1128  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
1129  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
1130 
1131  Gains=fopen(OutputGains.c_str(),"w");
1132  fprintf(Gains,"NEvents = %i\n",NEvent);
1133  fprintf(Gains,"NTracks = %i\n",NTrack);
1134  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1135  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
1136  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
1137  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1138  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
1139 
1140  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
1141 
1142  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
1143  stAPVGain* APV = APVsCollOrdered[a];
1144  if(APV==NULL)continue;
1145 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
1146  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);
1147 
1148  tree_Index = APV->Index;
1149  tree_Bin = (Charge_Vs_Index[elepos])->getTH2S()->GetXaxis()->FindBin(APV->Index);
1150  tree_DetId = APV->DetId;
1151  tree_APVId = APV->APVId;
1152  tree_SubDet = APV->SubDet;
1153  tree_x = APV->x;
1154  tree_y = APV->y;
1155  tree_z = APV->z;
1156  tree_Eta = APV->Eta;
1157  tree_R = APV->R;
1158  tree_Phi = APV->Phi;
1159  tree_Thickness = APV->Thickness;
1160  tree_FitMPV = APV->FitMPV;
1161  tree_FitMPVErr = APV->FitMPVErr;
1162  tree_FitWidth = APV->FitWidth;
1163  tree_FitWidthErr= APV->FitWidthErr;
1164  tree_FitChi2NDF = APV->FitChi2;
1165  tree_FitNorm = APV->FitNorm;
1166  tree_Gain = APV->Gain;
1167  tree_PrevGain = APV->PreviousGain;
1168  tree_PrevGainTick = APV->PreviousGainTick;
1169  tree_NEntries = APV->NEntries;
1170  tree_isMasked = APV->isMasked;
1171 
1172 
1173  if(tree_DetId==402673324){
1174  printf("%i | %i : %f --> %f (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
1175  }
1176 
1177 
1178  MyTree->Fill();
1179  }
1180  if(Gains)fclose(Gains);
1181 
1182 
1183 }
1184 
1185 bool SiStripGainFromCalibTree::produceTagFilter(){
1186 
1187  // The goal of this function is to check wether or not there is enough statistics
1188  // to produce a meaningful tag for the DB
1189  int elepos = (AlgoMode == "PCL")? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
1190  if( Charge_Vs_Index[elepos]==0 ) {
1191  edm::LogError("SiStripGainFromCalibTree") << "produceTagFilter -> Return false: could not retrieve the "
1192  << m_calibrationMode.c_str() << " statistics.\n"
1193  << "Please check if input contains "
1194  << m_calibrationMode.c_str() << " data." << std::endl;
1195  return false;
1196  }
1197 
1198 
1199  float integral = (Charge_Vs_Index[elepos])->getTH2S()->Integral();
1200  if( (Charge_Vs_Index[elepos])->getTH2S()->Integral(0,NStripAPVs+1, 0, 99999 ) < tagCondition_NClusters) {
1201  edm::LogWarning("SiStripGainFromCalibTree")
1202  << "calibrationMode -> " << m_calibrationMode << "\n"
1203  << "produceTagFilter -> Return false: Statistics is too low : " << integral << endl;
1204  return false;
1205  }
1206  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
1207  edm::LogWarning("SiStripGainFromCalibTree")
1208  << "calibrationMode -> " << m_calibrationMode << "\n"
1209  << "produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
1210  return false;
1211  }
1212  return true;
1213 }
1214 
1215 SiStripApvGain* SiStripGainFromCalibTree::getNewObject()
1216 {
1217  SiStripApvGain* obj = new SiStripApvGain();
1218  if(!m_harvestingMode) return obj;
1219 
1220  if(!produceTagFilter()){
1221  edm::LogWarning("SiStripGainFromCalibTree")<< "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
1222  setDoStore(false);
1223  return obj;
1224  }
1225 
1226 
1227  std::vector<float>* theSiStripVector = NULL;
1228  unsigned int PreviousDetId = 0;
1229  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
1230  stAPVGain* APV = APVsCollOrdered[a];
1231  if(APV==NULL){ printf("Bug\n"); continue; }
1232  if(APV->SubDet<=2)continue;
1233  if(APV->DetId != PreviousDetId){
1234  if(theSiStripVector!=NULL){
1235  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
1236  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
1237  }
1238  theSiStripVector = new std::vector<float>;
1239  PreviousDetId = APV->DetId;
1240  }
1241  theSiStripVector->push_back(APV->Gain);
1242  }
1243  if(theSiStripVector!=NULL){
1244  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
1245  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
1246  }
1247 
1248  return obj;
1249 }
1250 
1251 
1252 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
1253 {
1254 }
1255 
1256 void SiStripGainFromCalibTree::MakeCalibrationMap(){
1257  if(!useCalibration)return;
1258 
1259  TChain* t1 = new TChain("SiStripCalib/APVGain");
1260  t1->Add(m_calibrationPath.c_str());
1261 
1262  unsigned int tree_DetId;
1263  unsigned char tree_APVId;
1264  double tree_Gain;
1265 
1266  t1->SetBranchAddress("DetId" ,&tree_DetId );
1267  t1->SetBranchAddress("APVId" ,&tree_APVId );
1268  t1->SetBranchAddress("Gain" ,&tree_Gain );
1269 
1270  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
1271  t1->GetEntry(ientry);
1272  stAPVGain* APV = APVsColl[(tree_DetId<<4) | (unsigned int)tree_APVId];
1273  APV->CalibGain = tree_Gain;
1274  }
1275 }
1276 
1277 void
1278 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
1279 {
1280  // in AlCaHarvesting mode we just need to run the logic in the endJob step
1281  if(m_harvestingMode) return;
1282 
1283  if(AlgoMode=="CalibTree")return;
1284 
1285  eventnumber = iEvent.id().event();
1286  runnumber = iEvent.id().run();
1287  auto handle01 = connect(TrigTech , TrigTech_token_ , iEvent);
1288  auto handle02 = connect(trackchi2ndof , trackchi2ndof_token_ , iEvent);
1289  auto handle03 = connect(trackp , trackp_token_ , iEvent);
1290  auto handle04 = connect(trackpt , trackpt_token_ , iEvent);
1291  auto handle05 = connect(tracketa , tracketa_token_ , iEvent);
1292  auto handle06 = connect(trackphi , trackphi_token_ , iEvent);
1293  auto handle07 = connect(trackhitsvalid, trackhitsvalid_token_, iEvent);
1294  auto handle08 = connect(trackindex , trackindex_token_ , iEvent);
1295  auto handle09 = connect(rawid , rawid_token_ , iEvent);
1296  auto handle11 = connect(localdirx , localdirx_token_ , iEvent);
1297  auto handle12 = connect(localdiry , localdiry_token_ , iEvent);
1298  auto handle13 = connect(localdirz , localdirz_token_ , iEvent);
1299  auto handle14 = connect(firststrip , firststrip_token_ , iEvent);
1300  auto handle15 = connect(nstrips , nstrips_token_ , iEvent);
1301  auto handle16 = connect(saturation , saturation_token_ , iEvent);
1302  auto handle17 = connect(overlapping , overlapping_token_ , iEvent);
1303  auto handle18 = connect(farfromedge , farfromedge_token_ , iEvent);
1304  auto handle19 = connect(charge , charge_token_ , iEvent);
1305  auto handle21 = connect(path , path_token_ , iEvent);
1306  auto handle22 = connect(chargeoverpath, chargeoverpath_token_, iEvent);
1307  auto handle23 = connect(amplitude , amplitude_token_ , iEvent);
1308  auto handle24 = connect(gainused , gainused_token_ , iEvent);
1309 
1310  auto handle25 = connect(trackalgo , trackalgo_token_ , iEvent);
1311 
1312  processEvent();
1313 }
1314 
1315 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
SiStripGainFromCalibTree::trackindex_token_
edm::EDGetTokenT< std::vector< int > > trackindex_token_
Definition: SiStripGainFromCalibTree.cc:240
SiStripGainFromCalibTree::Charge_Vs_Index
std::vector< MonitorElement * > Charge_Vs_Index
Definition: SiStripGainFromCalibTree.cc:200
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
StripTopology::nstrips
virtual int nstrips() const =0
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
SiStripGainFromCalibTree::trackp
const std::vector< float > * trackp
Definition: SiStripGainFromCalibTree.cc:232
SiStripGainFromCalibTree::tracketa_token_
edm::EDGetTokenT< std::vector< double > > tracketa_token_
Definition: SiStripGainFromCalibTree.cc:234
SiStripGainFromCalibTree::amplitude_token_
edm::EDGetTokenT< std::vector< unsigned char > > amplitude_token_
Definition: SiStripGainFromCalibTree.cc:253
SiStripGainFromCalibTree::chargeoverpath
const std::vector< double > * chargeoverpath
Definition: SiStripGainFromCalibTree.cc:252
SiStripGainFromCalibTree::rawid_token_
edm::EDGetTokenT< std::vector< unsigned int > > rawid_token_
Definition: SiStripGainFromCalibTree.cc:241
SiStripGainFromCalibTree::amplitude
const std::vector< unsigned char > * amplitude
Definition: SiStripGainFromCalibTree.cc:253
SiStripGainFromCalibTree::isBFieldConsistentWithMode
bool isBFieldConsistentWithMode(const edm::EventSetup &iSetup) const
Definition: SiStripGainFromCalibTree.cc:587
or
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
SiStripGainFromCalibTree::getNewObject
SiStripApvGain * getNewObject() override
Definition: SiStripGainFromCalibTree.cc:1215
PixelTopology::ncolumns
virtual int ncolumns() const =0
SiStripGainFromCalibTree::localdirz
const std::vector< double > * localdirz
Definition: SiStripGainFromCalibTree.cc:244
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
SiStripGainFromCalibTree::trackchi2ndof
const std::vector< double > * trackchi2ndof
Definition: SiStripGainFromCalibTree.cc:231
SiStripGainFromCalibTree::algoBeginRun
virtual void algoBeginRun(const edm::Run &run, const edm::EventSetup &iSetup) override
Definition: SiStripGainFromCalibTree.cc:606
SiStripGainFromCalibTree::localdirx
const std::vector< double > * localdirx
Definition: SiStripGainFromCalibTree.cc:242
AlcaSiStripGainsAfterAbortGapHarvester_cff.DQMStore
tuple DQMStore
Definition: AlcaSiStripGainsAfterAbortGapHarvester_cff.py:11
processEvent
dispatcher processEvent(e, inputTag, standby)
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
SiStripGainFromCalibTree::algoEndRun
virtual void algoEndRun(const edm::Run &run, const edm::EventSetup &iSetup) override
Definition: SiStripGainFromCalibTree.cc:648
SiStripGainFromCalibTree::rawid
const std::vector< unsigned int > * rawid
Definition: SiStripGainFromCalibTree.cc:241
SiStripGainFromCalibTree::overlapping_token_
edm::EDGetTokenT< std::vector< bool > > overlapping_token_
Definition: SiStripGainFromCalibTree.cc:248
SiStripGainFromCalibTree::SiStripGainFromCalibTree
SiStripGainFromCalibTree(const edm::ParameterSet &)
Definition: SiStripGainFromCalibTree.cc:300
SiStripGainFromCalibTree::connect
edm::Handle< T > connect(const T *&ptr, edm::EDGetTokenT< T > token, const edm::Event &evt)
Definition: SiStripGainFromCalibTree.cc:152
PixelTopology::nrows
virtual int nrows() const =0
SiStripGainFromCalibTree::trackalgo
const std::vector< int > * trackalgo
Definition: SiStripGainFromCalibTree.cc:237
SiStripGainFromCalibTree::Charge_Vs_PathlengthTOB
std::vector< MonitorElement * > Charge_Vs_PathlengthTOB
Definition: SiStripGainFromCalibTree.cc:203
NULL
#define NULL
Definition: scimark2.h:8
SiStripGainFromCalibTree::storeOnTree
void storeOnTree(TFileService *tfs)
Definition: SiStripGainFromCalibTree.cc:1069
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
SiStripGainFromCalibTree::trackpt
const std::vector< float > * trackpt
Definition: SiStripGainFromCalibTree.cc:233
SiStripGainFromCalibTree::Charge_Vs_PathlengthTECM1
std::vector< MonitorElement * > Charge_Vs_PathlengthTECM1
Definition: SiStripGainFromCalibTree.cc:208
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
SiStripGainFromCalibTree::trackhitsvalid
const std::vector< unsigned int > * trackhitsvalid
Definition: SiStripGainFromCalibTree.cc:236
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
SiStripGainFromCalibTree::Charge_Vs_PathlengthTECM2
std::vector< MonitorElement * > Charge_Vs_PathlengthTECM2
Definition: SiStripGainFromCalibTree.cc:209
SiStripGainFromCalibTree::farfromedge_token_
edm::EDGetTokenT< std::vector< bool > > farfromedge_token_
Definition: SiStripGainFromCalibTree.cc:249
SiStripApvGain::put
bool put(const uint32_t &detID, Range input)
Definition: SiStripApvGain.cc:8
double_binary
Definition: DDStreamer.cc:227
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
SiStripGainFromCalibTree::trackhitsvalid_token_
edm::EDGetTokenT< std::vector< unsigned int > > trackhitsvalid_token_
Definition: SiStripGainFromCalibTree.cc:236
HLT_25ns10e33_v2_cff.InputTag
tuple InputTag
Definition: HLT_25ns10e33_v2_cff.py:58467
x
T x() const
Cartesian x coordinate.
Definition: Basic3DVectorLD.h:127
SiStripGainFromCalibTree::APVsCollOrdered
std::vector< stAPVGain * > APVsCollOrdered
Definition: SiStripGainFromCalibTree.cc:269
SiStripGainFromCalibTree::APVsColl
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
Definition: SiStripGainFromCalibTree.cc:270
SiStripGainFromCalibTree::path_token_
edm::EDGetTokenT< std::vector< double > > path_token_
Definition: SiStripGainFromCalibTree.cc:251
iEvent
int iEvent
Definition: GenABIO.cc:230
SiStripGainFromCalibTree::trackalgo_token_
edm::EDGetTokenT< std::vector< int > > trackalgo_token_
Definition: SiStripGainFromCalibTree.cc:237
GlobalPosition_Frontier_DevDB_cff.tag
tuple tag
Definition: GlobalPosition_Frontier_DevDB_cff.py:11
SiStripGainFromCalibTree::trackp_token_
edm::EDGetTokenT< std::vector< float > > trackp_token_
Definition: SiStripGainFromCalibTree.cc:232
SiStripGainFromCalibTree::IsGoodLandauFit
bool IsGoodLandauFit(double *FitResults)
Definition: SiStripGainFromCalibTree.cc:826
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
patZpeak.handle
tuple handle
Definition: patZpeak.py:22
SiStripGainFromCalibTree::algoBeginJob
virtual void algoBeginJob(const edm::EventSetup &iSetup) override
Definition: SiStripGainFromCalibTree.cc:451
SiStripApvGain::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripApvGain.h:27
j
int j
Definition: DBlmapReader.cc:9
SiStripGainFromCalibTree::TrigTech_token_
edm::EDGetTokenT< std::vector< bool > > TrigTech_token_
Definition: SiStripGainFromCalibTree.cc:228
SiStripGainFromCalibTree::chargeoverpath_token_
edm::EDGetTokenT< std::vector< double > > chargeoverpath_token_
Definition: SiStripGainFromCalibTree.cc:252
Exhume::I
const std::complex< double > I
Definition: I.h:8
SiStripGainFromCalibTree::gainused
const std::vector< double > * gainused
Definition: SiStripGainFromCalibTree.cc:254
SiStripGainFromCalibTree::Charge_Vs_PathlengthTECP1
std::vector< MonitorElement * > Charge_Vs_PathlengthTECP1
Definition: SiStripGainFromCalibTree.cc:206
SiStripGainFromCalibTree::localdirz_token_
edm::EDGetTokenT< std::vector< double > > localdirz_token_
Definition: SiStripGainFromCalibTree.cc:244
SiStripGainFromCalibTree::firststrip
const std::vector< unsigned short > * firststrip
Definition: SiStripGainFromCalibTree.cc:245
SiStripGainFromCalibTree::Charge_Vs_PathlengthTIDM
std::vector< MonitorElement * > Charge_Vs_PathlengthTIDM
Definition: SiStripGainFromCalibTree.cc:205
SiStripGainFromCalibTree::localdiry
const std::vector< double > * localdiry
Definition: SiStripGainFromCalibTree.cc:243
SiStripGainFromCalibTree::charge_token_
edm::EDGetTokenT< std::vector< unsigned int > > charge_token_
Definition: SiStripGainFromCalibTree.cc:250
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
funct::integral
Integral< F, X >::type integral(const F &f)
Definition: Integral.h:69
SiStripGainFromCalibTree::algoEndJob
virtual void algoEndJob() override
Definition: SiStripGainFromCalibTree.cc:755
SiStripGainFromCalibTree::nstrips_token_
edm::EDGetTokenT< std::vector< unsigned short > > nstrips_token_
Definition: SiStripGainFromCalibTree.cc:246
SiStripGainFromCalibTree::Charge_Vs_PathlengthTECP2
std::vector< MonitorElement * > Charge_Vs_PathlengthTECP2
Definition: SiStripGainFromCalibTree.cc:207
SiStripGainFromCalibTree::trackindex
const std::vector< int > * trackindex
Definition: SiStripGainFromCalibTree.cc:240
DetId
Definition: DetId.h:18
SiStripGainFromCalibTree::statCollectionFromMode
int statCollectionFromMode(const char *tag) const
Definition: SiStripGainFromCalibTree.cc:274
SiStripGainFromCalibTree::Charge_Vs_PathlengthTIB
std::vector< MonitorElement * > Charge_Vs_PathlengthTIB
Definition: SiStripGainFromCalibTree.cc:202
SiStripGainFromCalibTree::firststrip_token_
edm::EDGetTokenT< std::vector< unsigned short > > firststrip_token_
Definition: SiStripGainFromCalibTree.cc:245
SiStripGainFromCalibTree::algoAnalyze
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
Definition: SiStripGainFromCalibTree.cc:1278
SiStripGainFromCalibTree::trackpt_token_
edm::EDGetTokenT< std::vector< float > > trackpt_token_
Definition: SiStripGainFromCalibTree.cc:233
edm::Handle::product
T const * product() const
Definition: Handle.h:81
SiStripGainFromCalibTree::overlapping
const std::vector< bool > * overlapping
Definition: SiStripGainFromCalibTree.cc:248
SiStripGainFromCalibTree::saturation
const std::vector< bool > * saturation
Definition: SiStripGainFromCalibTree.cc:247
fed_dqm_sourceclient-live_cfg.path
string path
Definition: fed_dqm_sourceclient-live_cfg.py:33
edm::merge
bool merge(LuminosityBlockRange &lh, LuminosityBlockRange &rh)
Definition: LuminosityBlockRange.cc:78
SiStripGainFromCalibTree::getPeakOfLandau
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
Definition: SiStripGainFromCalibTree.cc:799
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
SiStripGainFromCalibTree::farfromedge
const std::vector< bool > * farfromedge
Definition: SiStripGainFromCalibTree.cc:249
SiStripGainFromCalibTree::TrigTech
const std::vector< bool > * TrigTech
Definition: SiStripGainFromCalibTree.cc:228
ConditionDBWriter< SiStripApvGain >::setDoStore
void setDoStore(const bool doStore)
When set to false the payload will not be written to the db.
Definition: ConditionDBWriter.h:379
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
a
double a
Definition: hdecay.h:121
SiStripGainFromCalibTree::localdirx_token_
edm::EDGetTokenT< std::vector< double > > localdirx_token_
Definition: SiStripGainFromCalibTree.cc:242
SiStripGainFromCalibTree::saturation_token_
edm::EDGetTokenT< std::vector< bool > > saturation_token_
Definition: SiStripGainFromCalibTree.cc:247
SiStripGainFromCalibTree::trackchi2ndof_token_
edm::EDGetTokenT< std::vector< double > > trackchi2ndof_token_
Definition: SiStripGainFromCalibTree.cc:231
SiStripGainFromCalibTree::Charge_Vs_PathlengthTIDP
std::vector< MonitorElement * > Charge_Vs_PathlengthTIDP
Definition: SiStripGainFromCalibTree.cc:204
SiStripGainFromCalibTree::localdiry_token_
edm::EDGetTokenT< std::vector< double > > localdiry_token_
Definition: SiStripGainFromCalibTree.cc:243
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
SiStripGainFromCalibTree::charge
const std::vector< unsigned int > * charge
Definition: SiStripGainFromCalibTree.cc:250
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:47
SiStripGainFromCalibTree::path
const std::vector< double > * path
Definition: SiStripGainFromCalibTree.cc:251
T
long double T
Definition: Basic3DVectorLD.h:57
SiStripGainFromCalibTree::trackphi
const std::vector< double > * trackphi
Definition: SiStripGainFromCalibTree.cc:235
SiStripGainFromCalibTree::nstrips
const std::vector< unsigned short > * nstrips
Definition: SiStripGainFromCalibTree.cc:246
SiStripGainFromCalibTree::gainused_token_
edm::EDGetTokenT< std::vector< double > > gainused_token_
Definition: SiStripGainFromCalibTree.cc:254
SiStripGainFromCalibTree::bookDQMHistos
void bookDQMHistos(const char *dqm_dir, const char *tag)
Definition: SiStripGainFromCalibTree.cc:411
SiStripGainFromCalibTree::trackphi_token_
edm::EDGetTokenT< std::vector< double > > trackphi_token_
Definition: SiStripGainFromCalibTree.cc:235
SiStripGainFromCalibTree::tracketa
const std::vector< double > * tracketa
Definition: SiStripGainFromCalibTree.cc:234
edm::Run
Definition: Run.h:43
Generated for CMSSW Reference Manual by doxygen 1.8.5