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 "TH2F.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  edm::LogInfo("SiStripGainFromCalibTree") << "Setting " << dqm_dir << "in DQM and booking histograms for tag "
414  << tag << std::endl;
415 
416  dbe->setCurrentFolder(dqm_dir);
417 
418  std::string stag(tag);
419  if(stag.size()!=0 && stag[0]!='_') stag.insert(0,1,'_');
420 
421  std::string cvi = std::string("Charge_Vs_Index") + stag;
422  //std::string cviA = std::string("Charge_Vs_Index_Absolute") + stag;
423  std::string cvpTIB = std::string("Charge_Vs_PathlengthTIB") + stag;
424  std::string cvpTOB = std::string("Charge_Vs_PathlengthTOB") + stag;
425  std::string cvpTIDP = std::string("Charge_Vs_PathlengthTIDP") + stag;
426  std::string cvpTIDM = std::string("Charge_Vs_PathlengthTIDM") + stag;
427  std::string cvpTECP1 = std::string("Charge_Vs_PathlengthTECP1") + stag;
428  std::string cvpTECP2 = std::string("Charge_Vs_PathlengthTECP2") + stag;
429  std::string cvpTECM1 = std::string("Charge_Vs_PathlengthTECM1") + stag;
430  std::string cvpTECM2 = std::string("Charge_Vs_PathlengthTECM2") + stag;
431 
432  int elepos = (m_harvestingMode && AlgoMode=="PCL")? Harvest : statCollectionFromMode(tag);
433 
434  Charge_Vs_Index[elepos] = dbe->book2D(cvi.c_str() , cvi.c_str() , 88625, 0 , 88624,2000,0,4000);
435  //Charge_Vs_Index_Absolute[elepos] = dbe->book2D(cviA.c_str() , cviA.c_str() , 88625, 0 , 88624,1000,0,4000);
436  Charge_Vs_PathlengthTIB[elepos] = dbe->book2D(cvpTIB.c_str() , cvpTIB.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
437  Charge_Vs_PathlengthTOB[elepos] = dbe->book2D(cvpTOB.c_str() , cvpTOB.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
438  Charge_Vs_PathlengthTIDP[elepos] = dbe->book2D(cvpTIDP.c_str() , cvpTIDP.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
439  Charge_Vs_PathlengthTIDM[elepos] = dbe->book2D(cvpTIDM.c_str() , cvpTIDM.c_str() , 20 , 0.3 , 1.3 , 250,0,2000);
440  Charge_Vs_PathlengthTECP1[elepos] = dbe->book2D(cvpTECP1.c_str(), cvpTECP1.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
441  Charge_Vs_PathlengthTECP2[elepos] = dbe->book2D(cvpTECP2.c_str(), cvpTECP2.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
442  Charge_Vs_PathlengthTECM1[elepos] = dbe->book2D(cvpTECM1.c_str(), cvpTECM1.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
443  Charge_Vs_PathlengthTECM2[elepos] = dbe->book2D(cvpTECM2.c_str(), cvpTECM2.c_str(), 20 , 0.3 , 1.3 , 250,0,2000);
444 }
445 
446 void SiStripGainFromCalibTree::algoBeginJob(const edm::EventSetup& iSetup)
447 {
448  edm::LogInfo("SiStripGainFromCalibTree") << "AlgoMode : " << AlgoMode << "\n"
449  << "CalibrationMode : " << m_calibrationMode << "\n"
450  << "HarvestingMode : " << m_harvestingMode << std::endl;
451  //Setup DQM histograms
452  if(AlgoMode != "PCL" or m_harvestingMode) {
453  const char * dqm_dir = "AlCaReco/SiStripGainsHarvesting/";
454  this->bookDQMHistos( dqm_dir, dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str() );
455  //if (m_harvestingMode) this->bookDQMHistos( dqm_dir, dqm_tag_[Harvest].c_str() );
456  //else this->bookDQMHistos( dqm_dir, dqm_tag_[statCollectionFromMode(m_calibrationMode.c_str())].c_str() );
457  } else {
458  std::string dqm_dir = m_DQMdir + ((m_splitDQMstat)? m_calibrationMode:"") + "/";
459  int elem = statCollectionFromMode(m_calibrationMode.c_str());
460  this->bookDQMHistos( dqm_dir.c_str(), dqm_tag_[elem].c_str() );
461  this->bookDQMHistos( dqm_dir.c_str(), dqm_tag_[( (elem%2)? elem-1: elem+1 )].c_str() );
462  }
463 
464 
465  edm::ESHandle<TrackerGeometry> tkGeom;
466  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
467  auto const & Det = tkGeom->dets();
468 
469  NPixelDets = 0;
470  NStripAPVs = 0;
471  unsigned int Index=0;
472  for(unsigned int i=0;i<Det.size();i++){
473  DetId Detid = Det[i]->geographicalId();
474  int SubDet = Detid.subdetId();
475 
476  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
477  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
478 
479  auto DetUnit = dynamic_cast<const StripGeomDetUnit*> (Det[i]);
480  if(!DetUnit)continue;
481 
482  const StripTopology& Topo = DetUnit->specificTopology();
483  unsigned int NAPV = Topo.nstrips()/128;
484 
485  for(unsigned int j=0;j<NAPV;j++){
486  stAPVGain* APV = new stAPVGain;
487  APV->Index = Index;
488  APV->Bin = -1;
489  APV->DetId = Detid.rawId();
490  APV->APVId = j;
491  APV->SubDet = SubDet;
492  APV->FitMPV = -1;
493  APV->FitMPVErr = -1;
494  APV->FitWidth = -1;
495  APV->FitWidthErr = -1;
496  APV->FitChi2 = -1;
497  APV->FitNorm = -1;
498  APV->Gain = -1;
499  APV->PreviousGain = 1;
500  APV->PreviousGainTick = 1;
501  APV->x = DetUnit->position().basicVector().x();
502  APV->y = DetUnit->position().basicVector().y();
503  APV->z = DetUnit->position().basicVector().z();
504  APV->Eta = DetUnit->position().basicVector().eta();
505  APV->Phi = DetUnit->position().basicVector().phi();
506  APV->R = DetUnit->position().basicVector().transverse();
507  APV->Thickness = DetUnit->surface().bounds().thickness();
508  APV->NEntries = 0;
509  APV->isMasked = false;
510 
511  APVsCollOrdered.push_back(APV);
512  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
513  Index++;
514  NStripAPVs++;
515  }
516  }
517  }
518 
519  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
520  DetId Detid = Det[i]->geographicalId();
521  int SubDet = Detid.subdetId();
522  if( SubDet == PixelSubdetector::PixelBarrel || PixelSubdetector::PixelEndcap ){
523  auto DetUnit = dynamic_cast<const PixelGeomDetUnit*> (Det[i]);
524  if(!DetUnit) continue;
525 
526  const PixelTopology& Topo = DetUnit->specificTopology();
527  unsigned int NROCRow = Topo.nrows()/(80.);
528  unsigned int NROCCol = Topo.ncolumns()/(52.);
529 
530  for(unsigned int j=0;j<NROCRow;j++){
531  for(unsigned int i=0;i<NROCCol;i++){
532 
533  stAPVGain* APV = new stAPVGain;
534  APV->Index = Index;
535  APV->Bin = -1;
536  APV->DetId = Detid.rawId();
537  APV->APVId = (j<<3 | i);
538  APV->SubDet = SubDet;
539  APV->FitMPV = -1;
540  APV->FitMPVErr = -1;
541  APV->FitWidth = -1;
542  APV->FitWidthErr = -1;
543  APV->FitChi2 = -1;
544  APV->Gain = -1;
545  APV->PreviousGain = 1;
546  APV->x = DetUnit->position().basicVector().x();
547  APV->y = DetUnit->position().basicVector().y();
548  APV->z = DetUnit->position().basicVector().z();
549  APV->Eta = DetUnit->position().basicVector().eta();
550  APV->Phi = DetUnit->position().basicVector().phi();
551  APV->R = DetUnit->position().basicVector().transverse();
552  APV->Thickness = DetUnit->surface().bounds().thickness();
553  APV->isMasked = false; //SiPixelQuality_->IsModuleBad(Detid.rawId());
554  APV->NEntries = 0;
555 
556  APVsCollOrdered.push_back(APV);
557  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
558  Index++;
559  NPixelDets++;
560  }}
561  }
562  }
563 
564 
565  MakeCalibrationMap();
566 
567  NEvent = 0;
568  NTrack = 0;
569  NClusterStrip = 0;
570  NClusterPixel = 0;
571  SRun = 1<<31;
572  ERun = 0;
573  GOOD = 0;
574  BAD = 0;
575  MASKED = 0;
576 }
577 
578 bool SiStripGainFromCalibTree::isBFieldConsistentWithMode( const edm::EventSetup& iSetup) const {
579  edm::ESHandle<RunInfo> runInfo;
580  iSetup.get<RunInfoRcd>().get(runInfo);
581 
582  double average_current = runInfo.product()->m_avg_current;
583  bool isOn = (average_current > MagFieldCurrentTh);
584  bool is0T = (m_calibrationMode.substr( m_calibrationMode.length()-2 )=="0T");
585 
586  return ( (isOn && !is0T) || (!isOn && is0T) );
587 }
588 
589 void SiStripGainFromCalibTree::swapBFieldMode() {
590  if ( m_calibrationMode.substr( m_calibrationMode.length()-2 )=="0T" ) {
591  m_calibrationMode.erase( m_calibrationMode.length()-2,2);
592  } else {
593  m_calibrationMode.append( "0T" );
594  }
595 }
596 
597 void SiStripGainFromCalibTree::algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup)
598 {
599  if( !m_harvestingMode && AlgoMode=="PCL") {
600  //Check consistency of calibration Mode and BField only for the ALCAPROMPT in the PCL workflow
601  if (!isBFieldConsistentWithMode(iSetup)) {
602  string prevMode = m_calibrationMode;
603  swapBFieldMode();
604  edm::LogInfo("SiStripGainFromCalibTree") << "Switching calibration mode for endorsing BField status: "
605  << prevMode << " ==> " << m_calibrationMode << std::endl;
606  }
607  }
608 
609  edm::ESHandle<SiStripGain> gainHandle;
610  iSetup.get<SiStripGainRcd>().get(gainHandle);
611  if(!gainHandle.isValid()){edm::LogError("SiStripGainFromCalibTree")<< "gainHandle is not valid\n"; exit(0);}
612 
613  edm::ESHandle<SiStripQuality> SiStripQuality_;
614  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
615 
616  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
617  stAPVGain* APV = APVsCollOrdered[a];
618 
619  APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId,APV->APVId);
620 // if(!FirstSetOfConstants){
621  if(gainHandle->getNumberOfTags()!=2){edm::LogError("SiStripGainFromCalibTree")<< "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";fflush(stdout);exit(0);};
622  float newPreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
623  if(APV->PreviousGain!=1 and newPreviousGain!=APV->PreviousGain)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: ParticleGain in the global tag changed\n";
624  APV->PreviousGain = newPreviousGain;
625 
626  float newPreviousGainTick = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 0),1);
627  if(APV->PreviousGainTick!=1 and newPreviousGainTick!=APV->PreviousGainTick)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: TickMarkGain in the global tag changed\n";
628  APV->PreviousGainTick = newPreviousGainTick;
629 
630 
631  //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f (G1) x %8.4f (G2)\n",APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain);
632 
633 
634 
635 // }
636  }
637 }
638 
639 void SiStripGainFromCalibTree::algoEndRun(const edm::Run& run, const edm::EventSetup& iSetup){
640  if(AlgoMode == "PCL" && !m_harvestingMode) return;//nothing to do in that case
641 
642  if(AlgoMode == "PCL" and m_harvestingMode){
643  // Load the 2D histograms from the DQM objects
644  // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
645  // DQMStore so that they can be used in the fit
646 
647  edm::LogInfo("SiStripGainFromCalibTree") << "Starting harvesting statistics" << std::endl;
648 
649  // check the required tag before adding histograms
650  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
651  if (elepos!=Harvest) {
652 
653  //collect statistics from DQM into the related monitored elements
654  std::string stag = m_calibrationMode;
655  if(stag.size()!=0 && stag[0]!='_') stag.insert(0,1,'_');
656 
657 
658  if (elepos==-1) {
659  //implememt backward compatibility
660  elepos=0;
661  stag="";
662  }
663 
664  std::string DQM_dir = m_DQMdir + ((m_splitDQMstat)? m_calibrationMode:"");
665 
666  std::string cvi = DQM_dir + std::string("/Charge_Vs_Index") + stag;
667  //std::string cviA = DQM_dir + std::string("/Charge_Vs_Index_Absolute") + stag;
668  std::string cvpTIB = DQM_dir + std::string("/Charge_Vs_PathlengthTIB") + stag;
669  std::string cvpTOB = DQM_dir + std::string("/Charge_Vs_PathlengthTOB") + stag;
670  std::string cvpTIDP = DQM_dir + std::string("/Charge_Vs_PathlengthTIDP") + stag;
671  std::string cvpTIDM = DQM_dir + std::string("/Charge_Vs_PathlengthTIDM") + stag;
672  std::string cvpTECP1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP1") + stag;
673  std::string cvpTECP2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECP2") + stag;
674  std::string cvpTECM1 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM1") + stag;
675  std::string cvpTECM2 = DQM_dir + std::string("/Charge_Vs_PathlengthTECM2") + stag;
676 
677  Charge_Vs_Index[elepos] = dbe->get(cvi.c_str());
678  //Charge_Vs_Index_Absolute[elepos] = dbe->get(cviA.c_str());
679  Charge_Vs_PathlengthTIB[elepos] = dbe->get(cvpTIB.c_str());
680  Charge_Vs_PathlengthTOB[elepos] = dbe->get(cvpTOB.c_str());
681  Charge_Vs_PathlengthTIDP[elepos] = dbe->get(cvpTIDP.c_str());
682  Charge_Vs_PathlengthTIDM[elepos] = dbe->get(cvpTIDM.c_str());
683  Charge_Vs_PathlengthTECP1[elepos] = dbe->get(cvpTECP1.c_str());
684  Charge_Vs_PathlengthTECP2[elepos] = dbe->get(cvpTECP2.c_str());
685  Charge_Vs_PathlengthTECM1[elepos] = dbe->get(cvpTECM1.c_str());
686  Charge_Vs_PathlengthTECM2[elepos] = dbe->get(cvpTECM2.c_str());
687 
688  if (Charge_Vs_Index[elepos]==0) {
689  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvi.c_str()
690  << ", statistics will not be summed!" << std::endl;
691  } else {
692  merge( (Charge_Vs_Index[Harvest])->getTH2F(), (Charge_Vs_Index[elepos])->getTH2F() );
693  edm::LogInfo("SiStripGainFromCalibTree") << "Harvesting "
694  << (Charge_Vs_Index[elepos])->getTH2F()->GetEntries() << " more clusters" << std::endl;
695  }
696 
697  //if (Charge_Vs_Index_Absolute[elepos]==0) {
698  // edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cviA.c_str()
699  // << ", statistics will not be summed!" << std::endl;
700  //} else merge( (Charge_Vs_Index_Absolute[Harvest])->getTH2F(), (Charge_Vs_Index_Absolute[elepos])->getTH2F() );
701 
702  if (Charge_Vs_PathlengthTIB[elepos]==0) {
703  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIB.c_str()
704  << ", statistics will not be summed!" << std::endl;
705  } else (Charge_Vs_PathlengthTIB[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTIB[elepos])->getTH2F());
706 
707  if (Charge_Vs_PathlengthTOB[elepos]==0) {
708  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTOB.c_str()
709  << ", statistics will not be summed!" << std::endl;
710  } else (Charge_Vs_PathlengthTOB[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTOB[elepos])->getTH2F());
711 
712  if (Charge_Vs_PathlengthTIDP[elepos]==0) {
713  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIDP.c_str()
714  << ", statistics will not be summed!" << std::endl;
715  } else (Charge_Vs_PathlengthTIDP[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTIDP[elepos])->getTH2F());
716 
717  if (Charge_Vs_PathlengthTIDM[elepos]==0) {
718  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTIDM.c_str()
719  << ", statistics will not be summed!" << std::endl;
720  } else (Charge_Vs_PathlengthTIDM[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTIDM[elepos])->getTH2F());
721 
722  if (Charge_Vs_PathlengthTECP1[elepos]==0) {
723  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECP1.c_str()
724  << ", statistics will not be summed!" << std::endl;
725  } else (Charge_Vs_PathlengthTECP1[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTECP1[elepos])->getTH2F());
726 
727  if (Charge_Vs_PathlengthTECP2[elepos]==0) {
728  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECP2.c_str()
729  << ", statistics will not be summed!" << std::endl;
730  } else (Charge_Vs_PathlengthTECP2[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTECP2[elepos])->getTH2F());
731 
732  if (Charge_Vs_PathlengthTECM1[elepos]==0) {
733  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECM1.c_str()
734  << ", statistics will not be summed!" << std::endl;
735  } else (Charge_Vs_PathlengthTECM1[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTECM1[elepos])->getTH2F());
736 
737  if (Charge_Vs_PathlengthTECM2[elepos]==0) {
738  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not retrieve " << cvpTECM2.c_str()
739  << ", statistics will not be summed!" << std::endl;
740  } else (Charge_Vs_PathlengthTECM2[Harvest])->getTH2F()->Add((Charge_Vs_PathlengthTECM2[elepos])->getTH2F());
741  }
742  }
743 }
744 
745 void
746 SiStripGainFromCalibTree::algoEndJob() {
747  if(AlgoMode == "PCL" && !m_harvestingMode) return;//nothing to do in that case
748 
749  if(AlgoMode == "CalibTree"){
750  edm::LogInfo("SiStripGainFromCalibTree") << "Analyzing calibration tree" << std::endl;
751  // Loop on calibTrees to fill the 2D histograms
752  algoAnalyzeTheTree();
753  }else if(m_harvestingMode){
754  NClusterStrip = (Charge_Vs_Index[Harvest])->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 );
755  NClusterPixel = (Charge_Vs_Index[Harvest])->getTH2F()->Integral(NStripAPVs+2, NStripAPVs+NPixelDets+2, 0, 99999 );
756  }
757 
758  // Now that we have the full statistics we can extract the information of the 2D histograms
759  algoComputeMPVandGain();
760 
761  if(AlgoMode != "PCL" or saveSummary){
762  edm::LogInfo("SiStripGainFromCalibTree") << "Saving summary into root file" << std::endl;
763 
764  //also save the 2D monitor elements to this file as TH2D tfs
765  tfs = edm::Service<TFileService>().operator->();
766 
767  //save only the statistics for the calibrationTag
768  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
769 
770  if( Charge_Vs_Index[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_Index[elepos])->getTH2F() );
771  //if( Charge_Vs_Index_Absolute[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_Index_Absolute[elepos])->getTH2F() );
772  if( Charge_Vs_PathlengthTIB[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTIB[elepos])->getTH2F() );
773  if( Charge_Vs_PathlengthTOB[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTOB[elepos])->getTH2F() );
774  if( Charge_Vs_PathlengthTIDP[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTIDP[elepos])->getTH2F() );
775  if( Charge_Vs_PathlengthTIDM[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTIDM[elepos])->getTH2F() );
776  if( Charge_Vs_PathlengthTECP1[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTECP1[elepos])->getTH2F() );
777  if( Charge_Vs_PathlengthTECP2[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTECP2[elepos])->getTH2F() );
778  if( Charge_Vs_PathlengthTECM1[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTECM1[elepos])->getTH2F() );
779  if( Charge_Vs_PathlengthTECM2[elepos]!=0 ) tfs->make<TH2F> ( *(Charge_Vs_PathlengthTECM2[elepos])->getTH2F() );
780 
781  storeOnTree(tfs);
782  }
783 }
784 
785 
786 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
787 {
788  FitResults[0] = -0.5; //MPV
789  FitResults[1] = 0; //MPV error
790  FitResults[2] = -0.5; //Width
791  FitResults[3] = 0; //Width error
792  FitResults[4] = -0.5; //Fit Chi2/NDF
793  FitResults[5] = 0; //Normalization
794 
795  if( InputHisto->GetEntries() < MinNrEntries)return;
796 
797  // perform fit with standard landau
798  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
799  MyLandau->SetParameter(1,300);
800  InputHisto->Fit(MyLandau,"0QR WW");
801 
802  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
803  FitResults[0] = MyLandau->GetParameter(1); //MPV
804  FitResults[1] = MyLandau->GetParError(1); //MPV error
805  FitResults[2] = MyLandau->GetParameter(2); //Width
806  FitResults[3] = MyLandau->GetParError(2); //Width error
807  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
808  FitResults[5] = MyLandau->GetParameter(0);
809 
810  delete MyLandau;
811 }
812 
813 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
814  if(FitResults[0] <= 0 )return false;
815 // if(FitResults[1] > MaxMPVError )return false;
816 // if(FitResults[4] > MaxChi2OverNDF)return false;
817  return true;
818 }
819 
820 void SiStripGainFromCalibTree::processEvent() {
821 
822  edm::LogInfo("SiStripGainFromCalibTree") << "Processing run " << runnumber
823  << " and event " << eventnumber
824  << " for " << m_calibrationMode << " calibration." << std::endl;
825 
826  if(runnumber<SRun)SRun=runnumber;
827  if(runnumber>ERun)ERun=runnumber;
828 
829  NEvent++;
830  NTrack+=(*trackp).size();
831 
832  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
833 
834  unsigned int FirstAmplitude=0;
835  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
836  FirstAmplitude+=(*nstrips)[i];
837  int TI = (*trackindex)[i];
838 
839  //printf("%i - %i - %i %i %i\n", (int)(*rawid)[i], (int)(*firststrip)[i]/128, (int)(*farfromedge)[i], (int)(*overlapping)[i], (int)(*saturation )[i] );
840  if((*tracketa )[TI] < MinTrackEta )continue;
841  if((*tracketa )[TI] > MaxTrackEta )continue;
842  if((*trackp )[TI] < MinTrackMomentum )continue;
843  if((*trackp )[TI] > MaxTrackMomentum )continue;
844  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
845  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
846  if((*trackalgo )[TI] > MaxTrackingIteration )continue;
847 
848  stAPVGain* APV = APVsColl[((*rawid)[i]<<4) | ((*firststrip)[i]/128)]; //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
849 
850  if(APV->SubDet>2 && (*farfromedge)[i] == false )continue;
851  if(APV->SubDet>2 && (*overlapping)[i] == true )continue;
852  if(APV->SubDet>2 && (*saturation )[i] && !AllowSaturation)continue;
853  if(APV->SubDet>2 && (*nstrips )[i] > MaxNrStrips )continue;
854 
855 
856  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
857 
858  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
859  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
860  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
861 
862  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
863  //printf("Amplitudes: ");
864  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
865  //printf("\n");
866 
867  if(APV->SubDet>2){NClusterStrip++;}else{NClusterPixel++;}
868 
869  int Charge = 0;
870  if(APV->SubDet>2 && (useCalibration || !FirstSetOfConstants)){
871  bool Saturation = false;
872  for(unsigned int s=0;s<(*nstrips)[i];s++){
873  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
874  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
875  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
876  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
877  if(StripCharge>1024){
878  StripCharge = 255;
879  Saturation = true;
880  }else if(StripCharge>254){
881  StripCharge = 254;
882  Saturation = true;
883  }
884  Charge += StripCharge;
885  }
886  if(Saturation && !AllowSaturation)continue;
887  }else if(APV->SubDet>2){
888  Charge = (*charge)[i];
889  }else{
890  Charge = (*charge)[i]/265.0; //expected scale factor between pixel and strip charge
891  }
892 
893  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
894 
895  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
896  if(APV->SubDet>2){
897  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
898  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
899  }
900  //(Charge_Vs_Index_Absolute[elepos])->Fill(APV->Index,Charge);
901  (Charge_Vs_Index[elepos]) ->Fill(APV->Index,ClusterChargeOverPath);
902 
903  if(APV->SubDet==StripSubdetector::TIB){ (Charge_Vs_PathlengthTIB[elepos]) ->Fill((*path)[i],Charge);
904  }else if(APV->SubDet==StripSubdetector::TOB){ (Charge_Vs_PathlengthTOB[elepos]) ->Fill((*path)[i],Charge);
905  }else if(APV->SubDet==StripSubdetector::TID){
906  if(APV->Eta<0){ (Charge_Vs_PathlengthTIDM[elepos]) ->Fill((*path)[i],Charge);
907  }else if(APV->Eta>0){ (Charge_Vs_PathlengthTIDP[elepos]) ->Fill((*path)[i],Charge);
908  }
909  }else if(APV->SubDet==StripSubdetector::TEC){
910  if(APV->Eta<0){
911  if(APV->Thickness<0.04){ (Charge_Vs_PathlengthTECM1[elepos])->Fill((*path)[i],Charge);
912  }else if(APV->Thickness>0.04){ (Charge_Vs_PathlengthTECM2[elepos])->Fill((*path)[i],Charge);
913  }
914  }else if(APV->Eta>0){
915  if(APV->Thickness<0.04){ (Charge_Vs_PathlengthTECP1[elepos])->Fill((*path)[i],Charge);
916  }else if(APV->Thickness>0.04){ (Charge_Vs_PathlengthTECP2[elepos])->Fill((*path)[i],Charge);
917  }
918  }
919  }
920 
921  }// END OF ON-CLUSTER LOOP
922 }//END OF processEvent()
923 
924 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
925 {
926  for(unsigned int i=0;i<VInputFiles.size();i++){
927  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
928  TFile *tfile = TFile::Open(VInputFiles[i].c_str());
929  TString tree_path = TString::Format("gainCalibrationTree%s/tree",m_calibrationMode.c_str());
930  TTree* tree = dynamic_cast<TTree*> (tfile->Get(tree_path.Data()));
931 
932  tree->SetBranchAddress((EventPrefix_ + "event" + EventSuffix_).c_str(), &eventnumber , NULL);
933  tree->SetBranchAddress((EventPrefix_ + "run" + EventSuffix_).c_str(), &runnumber , NULL);
934  tree->SetBranchAddress((EventPrefix_ + "TrigTech" + EventSuffix_).c_str(), &TrigTech , NULL);
935 
936  tree->SetBranchAddress((TrackPrefix_ + "chi2ndof" + TrackSuffix_).c_str(), &trackchi2ndof , NULL);
937  tree->SetBranchAddress((TrackPrefix_ + "momentum" + TrackSuffix_).c_str(), &trackp , NULL);
938  tree->SetBranchAddress((TrackPrefix_ + "pt" + TrackSuffix_).c_str(), &trackpt , NULL);
939  tree->SetBranchAddress((TrackPrefix_ + "eta" + TrackSuffix_).c_str(), &tracketa , NULL);
940  tree->SetBranchAddress((TrackPrefix_ + "phi" + TrackSuffix_).c_str(), &trackphi , NULL);
941  tree->SetBranchAddress((TrackPrefix_ + "hitsvalid" + TrackSuffix_).c_str(), &trackhitsvalid, NULL);
942  tree->SetBranchAddress((TrackPrefix_ + "algo" + TrackSuffix_).c_str(), &trackalgo , NULL);
943 
944  tree->SetBranchAddress((CalibPrefix_ + "trackindex" + CalibSuffix_).c_str(), &trackindex , NULL);
945  tree->SetBranchAddress((CalibPrefix_ + "rawid" + CalibSuffix_).c_str(), &rawid , NULL);
946  tree->SetBranchAddress((CalibPrefix_ + "localdirx" + CalibSuffix_).c_str(), &localdirx , NULL);
947  tree->SetBranchAddress((CalibPrefix_ + "localdiry" + CalibSuffix_).c_str(), &localdiry , NULL);
948  tree->SetBranchAddress((CalibPrefix_ + "localdirz" + CalibSuffix_).c_str(), &localdirz , NULL);
949  tree->SetBranchAddress((CalibPrefix_ + "firststrip" + CalibSuffix_).c_str(), &firststrip , NULL);
950  tree->SetBranchAddress((CalibPrefix_ + "nstrips" + CalibSuffix_).c_str(), &nstrips , NULL);
951  tree->SetBranchAddress((CalibPrefix_ + "saturation" + CalibSuffix_).c_str(), &saturation , NULL);
952  tree->SetBranchAddress((CalibPrefix_ + "overlapping" + CalibSuffix_).c_str(), &overlapping , NULL);
953  tree->SetBranchAddress((CalibPrefix_ + "farfromedge" + CalibSuffix_).c_str(), &farfromedge , NULL);
954  tree->SetBranchAddress((CalibPrefix_ + "charge" + CalibSuffix_).c_str(), &charge , NULL);
955  tree->SetBranchAddress((CalibPrefix_ + "path" + CalibSuffix_).c_str(), &path , NULL);
956  tree->SetBranchAddress((CalibPrefix_ + "chargeoverpath" + CalibSuffix_).c_str(), &chargeoverpath, NULL);
957  tree->SetBranchAddress((CalibPrefix_ + "amplitude" + CalibSuffix_).c_str(), &amplitude , NULL);
958  tree->SetBranchAddress((CalibPrefix_ + "gainused" + CalibSuffix_).c_str(), &gainused , NULL);
959 
960 
961  unsigned int nentries = tree->GetEntries();
962  printf("Number of Events = %i + %i = %i\n",NEvent,nentries,(NEvent+nentries));
963  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
964  printf("Looping on the Tree :");
965  int TreeStep = nentries/50;if(TreeStep<=1)TreeStep=1;
966  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
967  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
968  tree->GetEntry(ientry);
969  processEvent();
970  }printf("\n");// END OF EVENT LOOP
971  }
972 }
973 
974 
975 
976 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
977  unsigned int I=0;
978  TH1F* Proj = NULL;
979  double FitResults[6];
980  double MPVmean = 300;
981 
982  int elepos = (AlgoMode == "PCL")? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
983 
984  if ( Charge_Vs_Index[elepos]==0 ) {
985  edm::LogError("SiStripGainFromCalibTree") << "Harvesting: could not execute algoComputeMPVandGain method because "
986  << m_calibrationMode.c_str() << " statistics cannot be retrieved.\n"
987  << "Please check if input contains "
988  << m_calibrationMode.c_str() << " data." << std::endl;
989  return;
990  }
991 
992  TH2F *chvsidx = (Charge_Vs_Index[elepos])->getTH2F();
993 
994 
995  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
996  printf("Fitting Charge Distribution :");
997  int TreeStep = APVsColl.size()/50;
998  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
999  if(I%TreeStep==0){printf(".");fflush(stdout);}
1000  stAPVGain* APV = it->second;
1001  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
1002 
1003  if(APV->isMasked){APV->Gain=APV->PreviousGain; MASKED++; continue;}
1004 
1005  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
1006  if(!Proj)continue;
1007 
1008  if(CalibrationLevel==0){
1009  }else if(CalibrationLevel==1){
1010  int SecondAPVId = APV->APVId;
1011  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
1012  stAPVGain* APV2 = APVsColl[(APV->DetId<<4) | SecondAPVId];
1013  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1014  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
1015  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
1016  }else if(CalibrationLevel==2){
1017  for(unsigned int i=0;i<16;i++){ //loop up to 6APV for Strip and up to 16 for Pixels
1018  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
1019  tmpit = APVsColl.find((APV->DetId<<4) | i);
1020  if(tmpit==APVsColl.end())continue;
1021  stAPVGain* APV2 = tmpit->second;
1022  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
1023  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
1024  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
1025  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
1026  }
1027  }else{
1028  CalibrationLevel = 0;
1029  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
1030  }
1031 
1032  getPeakOfLandau(Proj,FitResults);
1033  APV->FitMPV = FitResults[0];
1034  APV->FitMPVErr = FitResults[1];
1035  APV->FitWidth = FitResults[2];
1036  APV->FitWidthErr = FitResults[3];
1037  APV->FitChi2 = FitResults[4];
1038  APV->FitNorm = FitResults[5];
1039  APV->NEntries = Proj->GetEntries();
1040 
1041  if(IsGoodLandauFit(FitResults)){
1042  APV->Gain = APV->FitMPV / MPVmean;
1043  if(APV->SubDet>2)GOOD++;
1044  }else{
1045  APV->Gain = APV->PreviousGain;
1046  if(APV->SubDet>2)BAD++;
1047  }
1048  if(APV->Gain<=0) APV->Gain = 1;
1049 
1050  //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);
1051  delete Proj;
1052  }printf("\n");
1053 }
1054 
1055 
1056 void SiStripGainFromCalibTree::storeOnTree(TFileService* tfs)
1057 {
1058  unsigned int tree_Index;
1059  unsigned int tree_Bin;
1060  unsigned int tree_DetId;
1061  unsigned char tree_APVId;
1062  unsigned char tree_SubDet;
1063  float tree_x;
1064  float tree_y;
1065  float tree_z;
1066  float tree_Eta;
1067  float tree_R;
1068  float tree_Phi;
1069  float tree_Thickness;
1070  float tree_FitMPV;
1071  float tree_FitMPVErr;
1072  float tree_FitWidth;
1073  float tree_FitWidthErr;
1074  float tree_FitChi2NDF;
1075  float tree_FitNorm;
1076  double tree_Gain;
1077  double tree_PrevGain;
1078  double tree_PrevGainTick;
1079  double tree_NEntries;
1080  bool tree_isMasked;
1081 
1082  TTree* MyTree;
1083  MyTree = tfs->make<TTree> ("APVGain","APVGain");
1084  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
1085  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
1086  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
1087  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
1088  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
1089  MyTree->Branch("x" ,&tree_x ,"x/F");
1090  MyTree->Branch("y" ,&tree_y ,"y/F");
1091  MyTree->Branch("z" ,&tree_z ,"z/F");
1092  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
1093  MyTree->Branch("R" ,&tree_R ,"R/F");
1094  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
1095  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
1096  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
1097  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
1098  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
1099  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
1100  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
1101  MyTree->Branch("FitNorm" ,&tree_FitNorm ,"FitNorm/F");
1102  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
1103  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
1104  MyTree->Branch("PrevGainTick" ,&tree_PrevGainTick,"PrevGainTick/D");
1105  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
1106  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
1107 
1108 
1109  FILE* Gains = stdout;
1110  fprintf(Gains,"NEvents = %i\n",NEvent);
1111  fprintf(Gains,"NTracks = %i\n",NTrack);
1112  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1113  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
1114  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1115  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
1116  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
1117 
1118  Gains=fopen(OutputGains.c_str(),"w");
1119  fprintf(Gains,"NEvents = %i\n",NEvent);
1120  fprintf(Gains,"NTracks = %i\n",NTrack);
1121  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
1122  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
1123  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
1124  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
1125  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
1126 
1127  int elepos = statCollectionFromMode(m_calibrationMode.c_str());
1128 
1129  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
1130  stAPVGain* APV = APVsCollOrdered[a];
1131  if(APV==NULL)continue;
1132 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
1133  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);
1134 
1135  tree_Index = APV->Index;
1136  tree_Bin = (Charge_Vs_Index[elepos])->getTH2F()->GetXaxis()->FindBin(APV->Index);
1137  tree_DetId = APV->DetId;
1138  tree_APVId = APV->APVId;
1139  tree_SubDet = APV->SubDet;
1140  tree_x = APV->x;
1141  tree_y = APV->y;
1142  tree_z = APV->z;
1143  tree_Eta = APV->Eta;
1144  tree_R = APV->R;
1145  tree_Phi = APV->Phi;
1146  tree_Thickness = APV->Thickness;
1147  tree_FitMPV = APV->FitMPV;
1148  tree_FitMPVErr = APV->FitMPVErr;
1149  tree_FitWidth = APV->FitWidth;
1150  tree_FitWidthErr= APV->FitWidthErr;
1151  tree_FitChi2NDF = APV->FitChi2;
1152  tree_FitNorm = APV->FitNorm;
1153  tree_Gain = APV->Gain;
1154  tree_PrevGain = APV->PreviousGain;
1155  tree_PrevGainTick = APV->PreviousGainTick;
1156  tree_NEntries = APV->NEntries;
1157  tree_isMasked = APV->isMasked;
1158 
1159 
1160  if(tree_DetId==402673324){
1161  printf("%i | %i : %f --> %f (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
1162  }
1163 
1164 
1165  MyTree->Fill();
1166  }
1167  if(Gains)fclose(Gains);
1168 
1169 
1170 }
1171 
1172 bool SiStripGainFromCalibTree::produceTagFilter(){
1173 
1174  // The goal of this function is to check wether or not there is enough statistics
1175  // to produce a meaningful tag for the DB
1176  int elepos = (AlgoMode == "PCL")? Harvest : statCollectionFromMode(m_calibrationMode.c_str());
1177  if( Charge_Vs_Index[elepos]==0 ) {
1178  edm::LogError("SiStripGainFromCalibTree") << "produceTagFilter -> Return false: could not retrieve the "
1179  << m_calibrationMode.c_str() << " statistics.\n"
1180  << "Please check if input contains "
1181  << m_calibrationMode.c_str() << " data." << std::endl;
1182  return false;
1183  }
1184 
1185 
1186  float integral = (Charge_Vs_Index[elepos])->getTH2F()->Integral();
1187  if( (Charge_Vs_Index[elepos])->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 ) < tagCondition_NClusters) {
1188  edm::LogWarning("SiStripGainFromCalibTree")
1189  << "calibrationMode -> " << m_calibrationMode << "\n"
1190  << "produceTagFilter -> Return false: Statistics is too low : " << integral << endl;
1191  return false;
1192  }
1193  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
1194  edm::LogWarning("SiStripGainFromCalibTree")
1195  << "calibrationMode -> " << m_calibrationMode << "\n"
1196  << "produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
1197  return false;
1198  }
1199  return true;
1200 }
1201 
1202 SiStripApvGain* SiStripGainFromCalibTree::getNewObject()
1203 {
1204  SiStripApvGain* obj = new SiStripApvGain();
1205  if(!m_harvestingMode) return obj;
1206 
1207  if(!produceTagFilter()){
1208  edm::LogWarning("SiStripGainFromCalibTree")<< "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
1209  setDoStore(false);
1210  return obj;
1211  }
1212 
1213 
1214  std::vector<float>* theSiStripVector = NULL;
1215  unsigned int PreviousDetId = 0;
1216  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
1217  stAPVGain* APV = APVsCollOrdered[a];
1218  if(APV==NULL){ printf("Bug\n"); continue; }
1219  if(APV->SubDet<=2)continue;
1220  if(APV->DetId != PreviousDetId){
1221  if(theSiStripVector!=NULL){
1222  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
1223  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
1224  }
1225  theSiStripVector = new std::vector<float>;
1226  PreviousDetId = APV->DetId;
1227  }
1228  theSiStripVector->push_back(APV->Gain);
1229  }
1230  if(theSiStripVector!=NULL){
1231  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
1232  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
1233  }
1234 
1235  return obj;
1236 }
1237 
1238 
1239 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
1240 {
1241 }
1242 
1243 void SiStripGainFromCalibTree::MakeCalibrationMap(){
1244  if(!useCalibration)return;
1245 
1246  TChain* t1 = new TChain("SiStripCalib/APVGain");
1247  t1->Add(m_calibrationPath.c_str());
1248 
1249  unsigned int tree_DetId;
1250  unsigned char tree_APVId;
1251  double tree_Gain;
1252 
1253  t1->SetBranchAddress("DetId" ,&tree_DetId );
1254  t1->SetBranchAddress("APVId" ,&tree_APVId );
1255  t1->SetBranchAddress("Gain" ,&tree_Gain );
1256 
1257  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
1258  t1->GetEntry(ientry);
1259  stAPVGain* APV = APVsColl[(tree_DetId<<4) | (unsigned int)tree_APVId];
1260  APV->CalibGain = tree_Gain;
1261  }
1262 }
1263 
1264 void
1265 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
1266 {
1267  // in AlCaHarvesting mode we just need to run the logic in the endJob step
1268  if(m_harvestingMode) return;
1269 
1270  if(AlgoMode=="CalibTree")return;
1271 
1272  eventnumber = iEvent.id().event();
1273  runnumber = iEvent.id().run();
1274  auto handle01 = connect(TrigTech , TrigTech_token_ , iEvent);
1275  auto handle02 = connect(trackchi2ndof , trackchi2ndof_token_ , iEvent);
1276  auto handle03 = connect(trackp , trackp_token_ , iEvent);
1277  auto handle04 = connect(trackpt , trackpt_token_ , iEvent);
1278  auto handle05 = connect(tracketa , tracketa_token_ , iEvent);
1279  auto handle06 = connect(trackphi , trackphi_token_ , iEvent);
1280  auto handle07 = connect(trackhitsvalid, trackhitsvalid_token_, iEvent);
1281  auto handle08 = connect(trackindex , trackindex_token_ , iEvent);
1282  auto handle09 = connect(rawid , rawid_token_ , iEvent);
1283  auto handle11 = connect(localdirx , localdirx_token_ , iEvent);
1284  auto handle12 = connect(localdiry , localdiry_token_ , iEvent);
1285  auto handle13 = connect(localdirz , localdirz_token_ , iEvent);
1286  auto handle14 = connect(firststrip , firststrip_token_ , iEvent);
1287  auto handle15 = connect(nstrips , nstrips_token_ , iEvent);
1288  auto handle16 = connect(saturation , saturation_token_ , iEvent);
1289  auto handle17 = connect(overlapping , overlapping_token_ , iEvent);
1290  auto handle18 = connect(farfromedge , farfromedge_token_ , iEvent);
1291  auto handle19 = connect(charge , charge_token_ , iEvent);
1292  auto handle21 = connect(path , path_token_ , iEvent);
1293  auto handle22 = connect(chargeoverpath, chargeoverpath_token_, iEvent);
1294  auto handle23 = connect(amplitude , amplitude_token_ , iEvent);
1295  auto handle24 = connect(gainused , gainused_token_ , iEvent);
1296 
1297  auto handle25 = connect(trackalgo , trackalgo_token_ , iEvent);
1298 
1299  processEvent();
1300 }
1301 
1302 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:578
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:1202
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:597
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:639
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:1056
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
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:813
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:446
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:746
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:1265
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:786
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
HLT_FULL_cff.InputTag
tuple InputTag
Definition: HLT_FULL_cff.py:72666
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:371
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
a
double a
Definition: hdecay.h:121
getTH2F
TH2F * getTH2F(std::string name, std::string process, std::string rootfolder, DQMStore *dbe_, bool verb, bool clone)
Definition: HcalHistoUtils.h:85
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