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 
6 #include "FWCore/Framework/interface/Frameworkfwd.h"
7 #include "FWCore/Framework/interface/EDAnalyzer.h"
8 #include "FWCore/Framework/interface/Event.h"
9 #include "FWCore/Framework/interface/MakerMacros.h"
10 #include "FWCore/Framework/interface/ESHandle.h"
11 #include "FWCore/ParameterSet/interface/ParameterSet.h"
12 #include "FWCore/MessageLogger/interface/MessageLogger.h"
13 #include "FWCore/Utilities/interface/Exception.h"
14 
15 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
16 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
17 #include "Geometry/CommonTopologies/interface/StripTopology.h"
18 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
19 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
20 
21 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
22 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
23 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
24 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
25 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
26 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
27 
28 #include "DataFormats/SiStripCluster/interface/SiStripClusterCollection.h"
29 
30 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
31 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
32 
33 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
34 #include "DataFormats/TrackReco/interface/Track.h"
35 #include "DataFormats/TrackReco/interface/TrackFwd.h"
36 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
37 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
38 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
39 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
40 #include "DataFormats/DetId/interface/DetId.h"
41 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
42 #include "DataFormats/TrackReco/interface/DeDxHit.h"
43 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
44 
45 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
46 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
47 
48 #include "TrackingTools/PatternTools/interface/Trajectory.h"
49 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
50 
51 #include "DQMServices/Core/interface/DQMStore.h"
52 #include "DQMServices/Core/interface/MonitorElement.h"
53 #include "FWCore/ServiceRegistry/interface/Service.h"
54 #include "CommonTools/UtilAlgos/interface/TFileService.h"
55 
56 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
57 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
58 
59 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
60 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
61 
62 
63 #include "TFile.h"
64 #include "TObjString.h"
65 #include "TString.h"
66 #include "TH1F.h"
67 #include "TH2F.h"
68 #include "TProfile.h"
69 #include "TF1.h"
70 #include "TROOT.h"
71 #include "TTree.h"
72 #include "TChain.h"
73 
74 #include <ext/hash_map>
75 
76 
77 
78 using namespace edm;
79 using namespace reco;
80 using namespace std;
81 using namespace __gnu_cxx;
82 using __gnu_cxx::hash_map;
83 using __gnu_cxx::hash;
84 
85 struct stAPVGain{
86  unsigned int Index;
87  int Bin;
88  unsigned int DetId;
89  unsigned int APVId;
90  unsigned int SubDet;
91  float x;
92  float y;
93  float z;
94  float Eta;
95  float R;
96  float Phi;
97  float Thickness;
98  double FitMPV;
99  double FitMPVErr;
100  double FitWidth;
101  double FitWidthErr;
102  double FitChi2;
103  double Gain;
104  double CalibGain;
105  double PreviousGain;
106  double NEntries;
107  TH1F* HCharge;
108  TH1F* HChargeN;
109  bool isMasked;
110 };
111 
112 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
113  public:
114  explicit SiStripGainFromCalibTree(const edm::ParameterSet&);
115  ~SiStripGainFromCalibTree();
116 
117 
118  private:
119 
120 
121  virtual void algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
122  virtual void algoEndJob () override;
123  virtual void algoAnalyze (const edm::Event &, const edm::EventSetup &) override;
124 
125 
126  void algoAnalyzeTheTree();
127  void algoComputeMPVandGain();
128 
129  bool IsFarFromBorder(TrajectoryStateOnSurface* trajState, const uint32_t detid, const edm::EventSetup* iSetup);
130  void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
131  bool IsGoodLandauFit(double* FitResults);
132  void storeOnTree();
133  void MakeCalibrationMap();
134  bool produceTagFilter();
135 
136 
137 
138  SiStripApvGain* getNewObject() override;
139 
140  TFileService *tfs;
141  DQMStore* dbe;
142  bool harvestingMode;
143  double MinNrEntries;
144  double MaxMPVError;
145  double MaxChi2OverNDF;
146  double MinTrackMomentum;
147  double MaxTrackMomentum;
148  double MinTrackEta;
149  double MaxTrackEta;
150  unsigned int MaxNrStrips;
151  unsigned int MinTrackHits;
152  double MaxTrackChiOverNdf;
153  bool AllowSaturation;
154  bool FirstSetOfConstants;
155  bool Validation;
156  bool OldGainRemoving;
157  int CalibrationLevel;
158 
159  bool useCalibration;
160  string m_calibrationPath;
161 
162  double tagCondition_NClusters;
163  double tagCondition_GoodFrac;
164 
165  edm::InputTag theTracksLabel;
166  std::string AlgoMode;
167  std::string OutputGains;
168  vector<string> VInputFiles;
169 
170  MonitorElement* Charge_Vs_Index;
171  MonitorElement* Charge_Vs_Index_Absolute;
172  MonitorElement* Charge_Vs_PathlengthTIB;
173  MonitorElement* Charge_Vs_PathlengthTOB;
174  MonitorElement* Charge_Vs_PathlengthTIDP;
175  MonitorElement* Charge_Vs_PathlengthTIDM;
176  MonitorElement* Charge_Vs_PathlengthTECP1;
177  MonitorElement* Charge_Vs_PathlengthTECP2;
178  MonitorElement* Charge_Vs_PathlengthTECM1;
179  MonitorElement* Charge_Vs_PathlengthTECM2;
180 
181  unsigned int NEvent;
182  unsigned int NTrack;
183  unsigned int NCluster;
184  unsigned int SRun;
185  unsigned int ERun;
186  unsigned int GOOD;
187  unsigned int BAD;
188  unsigned int MASKED;
189 
190 
191 
192  private :
193  class isEqual{
194  public:
195  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
196  };
197 
198  std::vector<stAPVGain*> APVsCollOrdered;
199  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
200 };
201 
202 SiStripGainFromCalibTree::SiStripGainFromCalibTree(const edm::ParameterSet& iConfig) : ConditionDBWriter<SiStripApvGain>(iConfig)
203 {
204  OutputGains = iConfig.getParameter<std::string>("OutputGains");
205  theTracksLabel = iConfig.getUntrackedParameter<edm::InputTag>("Tracks");
206 
207  AlgoMode = iConfig.getUntrackedParameter<std::string>("AlgoMode", "CalibTree");
208  MinNrEntries = iConfig.getUntrackedParameter<double> ("minNrEntries" , 20);
209  MaxMPVError = iConfig.getUntrackedParameter<double> ("maxMPVError" , 500.0);
210  MaxChi2OverNDF = iConfig.getUntrackedParameter<double> ("maxChi2OverNDF" , 5.0);
211  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 3.0);
212  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
213  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
214  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
215  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 2);
216  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 8);
217  MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double> ("MaxTrackChiOverNdf" , 3);
218  AllowSaturation = iConfig.getUntrackedParameter<bool> ("AllowSaturation" , false);
219  FirstSetOfConstants = iConfig.getUntrackedParameter<bool> ("FirstSetOfConstants", true);
220  Validation = iConfig.getUntrackedParameter<bool> ("Validation" , false);
221  OldGainRemoving = iConfig.getUntrackedParameter<bool> ("OldGainRemoving" , false);
222 
223  CalibrationLevel = iConfig.getUntrackedParameter<int> ("CalibrationLevel" , 0);
224  VInputFiles = iConfig.getParameter<vector<string> > ("InputFiles");
225 
226 
227  useCalibration = iConfig.getUntrackedParameter<bool> ("UseCalibration" , false);
228  m_calibrationPath = iConfig.getUntrackedParameter<string> ("calibrationPath");
229  harvestingMode = iConfig.getUntrackedParameter<bool> ("harvestingMode" , false);
230 
231  tagCondition_NClusters = iConfig.getUntrackedParameter<double> ("NClustersForTagProd" , 2E8);
232  tagCondition_GoodFrac = iConfig.getUntrackedParameter<double> ("GoodFracForTagProd" , 0.95);
233 
234  dbe = edm::Service<DQMStore>().operator->();
235  dbe->setVerbose(10);
236 }
237 
238 void SiStripGainFromCalibTree::algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup)
239 {
240  cout << "algoBeginRun start" << endl;
241  if(!harvestingMode){
242  cout << " booking start" << endl;
243  dbe->setCurrentFolder("AlCaReco/SiStripGains/");
244  Charge_Vs_Index = dbe->book2D("Charge_Vs_Index" , "Charge_Vs_Index" , 72785, 0 , 72784,1000,0,2000);
245  Charge_Vs_Index_Absolute = dbe->book2D("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 72785, 0 , 72784, 500,0,2000);
246  Charge_Vs_PathlengthTIB = dbe->book2D("Charge_Vs_PathlengthTIB" , "Charge_Vs_PathlengthTIB" , 20 , 0.3 , 1.3 , 250,0,2000);
247  Charge_Vs_PathlengthTOB = dbe->book2D("Charge_Vs_PathlengthTOB" , "Charge_Vs_PathlengthTOB" , 20 , 0.3 , 1.3 , 250,0,2000);
248  Charge_Vs_PathlengthTIDP = dbe->book2D("Charge_Vs_PathlengthTIDP" , "Charge_Vs_PathlengthTIDP" , 20 , 0.3 , 1.3 , 250,0,2000);
249  Charge_Vs_PathlengthTIDM = dbe->book2D("Charge_Vs_PathlengthTIDM" , "Charge_Vs_PathlengthTIDM" , 20 , 0.3 , 1.3 , 250,0,2000);
250  Charge_Vs_PathlengthTECP1 = dbe->book2D("Charge_Vs_PathlengthTECP1", "Charge_Vs_PathlengthTECP1", 20 , 0.3 , 1.3 , 250,0,2000);
251  Charge_Vs_PathlengthTECP2 = dbe->book2D("Charge_Vs_PathlengthTECP2", "Charge_Vs_PathlengthTECP2", 20 , 0.3 , 1.3 , 250,0,2000);
252  Charge_Vs_PathlengthTECM1 = dbe->book2D("Charge_Vs_PathlengthTECM1", "Charge_Vs_PathlengthTECM1", 20 , 0.3 , 1.3 , 250,0,2000);
253  Charge_Vs_PathlengthTECM2 = dbe->book2D("Charge_Vs_PathlengthTECM2", "Charge_Vs_PathlengthTECM2", 20 , 0.3 , 1.3 , 250,0,2000);
254  }
255  edm::ESHandle<TrackerGeometry> tkGeom;
256  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
257  auto const & Det = tkGeom->dets();
258 
259  edm::ESHandle<SiStripGain> gainHandle;
260  iSetup.get<SiStripGainRcd>().get(gainHandle);
261  if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}
262 
263  for(unsigned int i=0;i<gainHandle->getNumberOfTags();i++){
264  printf("Reccord %i --> Rcd Name = %s Label Name = %s\n",i,gainHandle->getRcdName(i).c_str(), gainHandle->getLabelName(i).c_str());
265  }
266 
267  edm::ESHandle<SiStripQuality> SiStripQuality_;
268  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
269 
270  unsigned int Index=0;
271  for(unsigned int i=0;i<Det.size();i++){
272  DetId Detid = Det[i]->geographicalId();
273  int SubDet = Detid.subdetId();
274 
275  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
276  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
277 
278  auto DetUnit = dynamic_cast<const StripGeomDetUnit*> (Det[i]);
279  if(!DetUnit)continue;
280 
281  const StripTopology& Topo = DetUnit->specificTopology();
282  unsigned int NAPV = Topo.nstrips()/128;
283  for(unsigned int j=0;j<NAPV;j++){
284  stAPVGain* APV = new stAPVGain;
285  APV->Index = Index;
286  APV->Bin = -1;
287  APV->DetId = Detid.rawId();
288  APV->APVId = j;
289  APV->SubDet = SubDet;
290  APV->FitMPV = -1;
291  APV->FitMPVErr = -1;
292  APV->FitWidth = -1;
293  APV->FitWidthErr = -1;
294  APV->FitChi2 = -1;
295  APV->Gain = -1;
296  APV->PreviousGain = 1;
297  APV->x = DetUnit->position().basicVector().x();
298  APV->y = DetUnit->position().basicVector().y();
299  APV->z = DetUnit->position().basicVector().z();
300  APV->Eta = DetUnit->position().basicVector().eta();
301  APV->Phi = DetUnit->position().basicVector().phi();
302  APV->R = DetUnit->position().basicVector().transverse();
303  APV->Thickness = DetUnit->surface().bounds().thickness();
304  APV->NEntries = 0;
305  APV->isMasked = SiStripQuality_->IsApvBad(Detid.rawId(),j);
306 
307  if(!FirstSetOfConstants){
308  if(gainHandle->getNumberOfTags()!=2){printf("ERROR: NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n");fflush(stdout);exit(0);};
309  APV->PreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
310  //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f\n",APV->DetId,APV->APVId,APV->PreviousGain);
311  }
312 
313  APVsCollOrdered.push_back(APV);
314  APVsColl[(APV->DetId<<3) | APV->APVId] = APV;
315  Index++;
316  }
317  }
318  }
319 
320  MakeCalibrationMap();
321 
322  NEvent = 0;
323  NTrack = 0;
324  NCluster = 0;
325  SRun = 1<<31;
326  ERun = 0;
327  GOOD = 0;
328  BAD = 0;
329  MASKED = 0;
330 
331  cout << "algoBeginRun end" << endl;
332 }
333 
334 void
335 SiStripGainFromCalibTree::algoEndJob() {
336  if(AlgoMode == "PCL" && !harvestingMode)return;//nothing to do in that case
337 
338  if(AlgoMode=="CalibTree"){
339  // Loop on calibTrees to fill the 2D histograms
340  algoAnalyzeTheTree();
341  }else if(harvestingMode){
342  // Load the 2D histograms from the DQM objects
343  // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
344  // DQMStore so that they can be used in the fit
345  cout << " retrieving from DQMStore" << endl;
346  Charge_Vs_Index = dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index");
347  Charge_Vs_Index_Absolute = dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index_Absolute");
348  Charge_Vs_PathlengthTIB = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIB");
349  Charge_Vs_PathlengthTOB = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTOB");
350  Charge_Vs_PathlengthTIDP = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDP");
351  Charge_Vs_PathlengthTIDM = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDM");
352  Charge_Vs_PathlengthTECP1 = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP1");
353  Charge_Vs_PathlengthTECP2 = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP2");
354  Charge_Vs_PathlengthTECM1 = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM1");
355  Charge_Vs_PathlengthTECM2 = dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM2");
356  }
357 
358  // Now that we have the full statistics we can extract the information of the 2D histograms
359  algoComputeMPVandGain();
360 
361  //FIXME: for the moment the tree is disabled in PCL, among other reasons the fact that the TFileSevice is not available @ Tier0
362  if(AlgoMode != "PCL"){
363  storeOnTree();
364  //also save the 2D monitor elements to this file as TH2D tfs
365  tfs = edm::Service<TFileService>().operator->();
366  tfs->make<TH2F> (*Charge_Vs_Index->getTH2F());
367  tfs->make<TH2F> (*Charge_Vs_Index_Absolute->getTH2F());
368  tfs->make<TH2F> (*Charge_Vs_PathlengthTIB->getTH2F());
369  tfs->make<TH2F> (*Charge_Vs_PathlengthTOB->getTH2F());
370  tfs->make<TH2F> (*Charge_Vs_PathlengthTIDP->getTH2F());
371  tfs->make<TH2F> (*Charge_Vs_PathlengthTIDM->getTH2F());
372  tfs->make<TH2F> (*Charge_Vs_PathlengthTECP1->getTH2F());
373  tfs->make<TH2F> (*Charge_Vs_PathlengthTECP2->getTH2F());
374  tfs->make<TH2F> (*Charge_Vs_PathlengthTECM1->getTH2F());
375  tfs->make<TH2F> (*Charge_Vs_PathlengthTECM2->getTH2F());
376  }
377 }
378 
379 
380 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
381 {
382  FitResults[0] = -0.5; //MPV
383  FitResults[1] = 0; //MPV error
384  FitResults[2] = -0.5; //Width
385  FitResults[3] = 0; //Width error
386  FitResults[4] = -0.5; //Fit Chi2/NDF
387 
388  if( InputHisto->GetEntries() < MinNrEntries)return;
389 
390  // perform fit with standard landau
391  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
392  MyLandau->SetParameter(1,300);
393  InputHisto->Fit(MyLandau,"0QR WW");
394 
395  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
396  FitResults[0] = MyLandau->GetParameter(1); //MPV
397  FitResults[1] = MyLandau->GetParError(1); //MPV error
398  FitResults[2] = MyLandau->GetParameter(2); //Width
399  FitResults[3] = MyLandau->GetParError(2); //Width error
400  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
401 
402  delete MyLandau;
403 }
404 
405 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
406  if(FitResults[0] <= 0 )return false;
407 // if(FitResults[1] > MaxMPVError )return false;
408 // if(FitResults[4] > MaxChi2OverNDF)return false;
409  return true;
410 }
411 
412 
413 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
414 {
415  for(unsigned int i=0;i<VInputFiles.size();i++){
416  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
417  TChain* tree = new TChain("gainCalibrationTree/tree");
418  tree->Add(VInputFiles[i].c_str());
419 
420  TString EventPrefix("");
421  TString EventSuffix("");
422 
423  TString TrackPrefix("track");
424  TString TrackSuffix("");
425 
426  TString CalibPrefix("GainCalibration");
427  TString CalibSuffix("");
428 
429  unsigned int eventnumber = 0; tree->SetBranchAddress(EventPrefix + "event" + EventSuffix, &eventnumber , NULL);
430  unsigned int runnumber = 0; tree->SetBranchAddress(EventPrefix + "run" + EventSuffix, &runnumber , NULL);
431  std::vector<bool>* TrigTech = 0; tree->SetBranchAddress(EventPrefix + "TrigTech" + EventSuffix, &TrigTech , NULL);
432 
433  std::vector<double>* trackchi2ndof = 0; tree->SetBranchAddress(TrackPrefix + "chi2ndof" + TrackSuffix, &trackchi2ndof , NULL);
434  std::vector<float>* trackp = 0; tree->SetBranchAddress(TrackPrefix + "momentum" + TrackSuffix, &trackp , NULL);
435  std::vector<float>* trackpt = 0; tree->SetBranchAddress(TrackPrefix + "pt" + TrackSuffix, &trackpt , NULL);
436  std::vector<double>* tracketa = 0; tree->SetBranchAddress(TrackPrefix + "eta" + TrackSuffix, &tracketa , NULL);
437  std::vector<double>* trackphi = 0; tree->SetBranchAddress(TrackPrefix + "phi" + TrackSuffix, &trackphi , NULL);
438  std::vector<unsigned int>* trackhitsvalid = 0; tree->SetBranchAddress(TrackPrefix + "hitsvalid" + TrackSuffix, &trackhitsvalid, NULL);
439 
440  std::vector<int>* trackindex = 0; tree->SetBranchAddress(CalibPrefix + "trackindex" + CalibSuffix, &trackindex , NULL);
441  std::vector<unsigned int>* rawid = 0; tree->SetBranchAddress(CalibPrefix + "rawid" + CalibSuffix, &rawid , NULL);
442  std::vector<float>* localdirx = 0; tree->SetBranchAddress(CalibPrefix + "localdirx" + CalibSuffix, &localdirx , NULL);
443  std::vector<float>* localdiry = 0; tree->SetBranchAddress(CalibPrefix + "localdiry" + CalibSuffix, &localdiry , NULL);
444  std::vector<float>* localdirz = 0; tree->SetBranchAddress(CalibPrefix + "localdirz" + CalibSuffix, &localdirz , NULL);
445  std::vector<unsigned short>* firststrip = 0; tree->SetBranchAddress(CalibPrefix + "firststrip" + CalibSuffix, &firststrip , NULL);
446  std::vector<unsigned short>* nstrips = 0; tree->SetBranchAddress(CalibPrefix + "nstrips" + CalibSuffix, &nstrips , NULL);
447  std::vector<bool>* saturation = 0; tree->SetBranchAddress(CalibPrefix + "saturation" + CalibSuffix, &saturation , NULL);
448  std::vector<bool>* overlapping = 0; tree->SetBranchAddress(CalibPrefix + "overlapping" + CalibSuffix, &overlapping , NULL);
449  std::vector<bool>* farfromedge = 0; tree->SetBranchAddress(CalibPrefix + "farfromedge" + CalibSuffix, &farfromedge , NULL);
450  std::vector<unsigned int>* charge = 0; tree->SetBranchAddress(CalibPrefix + "charge" + CalibSuffix, &charge , NULL);
451  std::vector<float>* path = 0; tree->SetBranchAddress(CalibPrefix + "path" + CalibSuffix, &path , NULL);
452  std::vector<float>* chargeoverpath = 0; tree->SetBranchAddress(CalibPrefix + "chargeoverpath" + CalibSuffix, &chargeoverpath, NULL);
453  std::vector<unsigned char>* amplitude = 0; tree->SetBranchAddress(CalibPrefix + "amplitude" + CalibSuffix, &amplitude , NULL);
454  std::vector<double>* gainused = 0; tree->SetBranchAddress(CalibPrefix + "gainused" + CalibSuffix, &gainused , NULL);
455 
456 
457  printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));
458  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
459  printf("Looping on the Tree :");
460  int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
461  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
462  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
463  tree->GetEntry(ientry);
464 
465  if(runnumber<SRun)SRun=runnumber;
466  if(runnumber>ERun)ERun=runnumber;
467 
468  NEvent++;
469  NTrack+=(*trackp).size();
470 
471  unsigned int FirstAmplitude=0;
472  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
473  FirstAmplitude+=(*nstrips)[i];
474  int TI = (*trackindex)[i];
475  if((*tracketa )[TI] < MinTrackEta )continue;
476  if((*tracketa )[TI] > MaxTrackEta )continue;
477  if((*trackp )[TI] < MinTrackMomentum )continue;
478  if((*trackp )[TI] > MaxTrackMomentum )continue;
479  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
480  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
481  if((*farfromedge)[i] == false )continue;
482  if((*overlapping)[i] == true )continue;
483  if((*saturation )[i] && !AllowSaturation)continue;
484  if((*nstrips )[i] > MaxNrStrips )continue;
485  NCluster++;
486 
487  stAPVGain* APV = APVsColl[((*rawid)[i]<<3) | ((*firststrip)[i]/128)];
488  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
489 
490  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
491  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
492  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
493 
494  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
495  //printf("Amplitudes: ");
496  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
497  //printf("\n");
498 
499 
500  int Charge = 0;
501  if(useCalibration || !FirstSetOfConstants){
502  bool Saturation = false;
503  for(unsigned int s=0;s<(*nstrips)[i];s++){
504  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
505  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
506  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
507  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
508  if(StripCharge>1024){
509  StripCharge = 255;
510  Saturation = true;
511  }else if(StripCharge>254){
512  StripCharge = 254;
513  Saturation = true;
514  }
515  Charge += StripCharge;
516  }
517  if(Saturation && !AllowSaturation)continue;
518  }else{
519  Charge = (*charge)[i];
520  }
521 
522  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
523 
524  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
525  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
526  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
527  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
528  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
529 
530 
531  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill((*path)[i],Charge);
532  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill((*path)[i],Charge);
533  }else if(APV->SubDet==StripSubdetector::TID){
534  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
535  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
536  }
537  }else if(APV->SubDet==StripSubdetector::TEC){
538  if(APV->Eta<0){
539  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
540  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
541  }
542  }else if(APV->Eta>0){
543  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
544  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
545  }
546  }
547  }
548 
549  }// END OF ON-CLUSTER LOOP
550  }printf("\n");// END OF EVENT LOOP
551 
552  }
553 }
554 
555 
556 
557 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
558  unsigned int I=0;
559  TH1F* Proj = NULL;
560  double FitResults[5];
561  double MPVmean = 300;
562 
563  TH2F *chvsidx = Charge_Vs_Index->getTH2F();
564 
565 
566  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
567  printf("Fitting Charge Distribution :");
568  int TreeStep = APVsColl.size()/50;
569  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
570  if(I%TreeStep==0){printf(".");fflush(stdout);}
571  stAPVGain* APV = it->second;
572  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
573 
574  if(APV->isMasked){MASKED++; continue;}
575 
576  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
577  if(!Proj)continue;
578 
579  if(CalibrationLevel==0){
580  }else if(CalibrationLevel==1){
581  int SecondAPVId = APV->APVId;
582  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
583  stAPVGain* APV2 = APVsColl[(APV->DetId<<3) | SecondAPVId];
584  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
585  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
586  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
587  }else if(CalibrationLevel==2){
588  for(unsigned int i=0;i<6;i++){
589  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
590  tmpit = APVsColl.find((APV->DetId<<3) | i);
591  if(tmpit==APVsColl.end())continue;
592  stAPVGain* APV2 = tmpit->second;
593  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
594  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
595  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
596  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
597  }
598  }else{
599  CalibrationLevel = 0;
600  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
601  }
602 
603  getPeakOfLandau(Proj,FitResults);
604  APV->FitMPV = FitResults[0];
605  APV->FitMPVErr = FitResults[1];
606  APV->FitWidth = FitResults[2];
607  APV->FitWidthErr = FitResults[3];
608  APV->FitChi2 = FitResults[4];
609  APV->NEntries = Proj->GetEntries();
610 
611  if(IsGoodLandauFit(FitResults)){
612  APV->Gain = APV->FitMPV / MPVmean;
613  GOOD++;
614  }else{
615  APV->Gain = APV->PreviousGain;
616  BAD++;
617  }
618  if(APV->Gain<=0) APV->Gain = 1;
619 
620  //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);
621  delete Proj;
622  }printf("\n");
623 }
624 
625 
626 void SiStripGainFromCalibTree::storeOnTree()
627 {
628  tfs = edm::Service<TFileService>().operator->();
629 
630  unsigned int tree_Index;
631  unsigned int tree_Bin;
632  unsigned int tree_DetId;
633  unsigned char tree_APVId;
634  unsigned char tree_SubDet;
635  float tree_x;
636  float tree_y;
637  float tree_z;
638  float tree_Eta;
639  float tree_R;
640  float tree_Phi;
641  float tree_Thickness;
642  float tree_FitMPV;
643  float tree_FitMPVErr;
644  float tree_FitWidth;
645  float tree_FitWidthErr;
646  float tree_FitChi2NDF;
647  double tree_Gain;
648  double tree_PrevGain;
649  double tree_NEntries;
650  bool tree_isMasked;
651 
652  TTree* MyTree;
653  MyTree = tfs->make<TTree> ("APVGain","APVGain");
654  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
655  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
656  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
657  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
658  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
659  MyTree->Branch("x" ,&tree_x ,"x/F");
660  MyTree->Branch("y" ,&tree_y ,"y/F");
661  MyTree->Branch("z" ,&tree_z ,"z/F");
662  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
663  MyTree->Branch("R" ,&tree_R ,"R/F");
664  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
665  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
666  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
667  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
668  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
669  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
670  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
671  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
672  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
673  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
674  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
675 
676 
677  FILE* Gains = stdout;
678  if(AlgoMode!="PCL"){
679  fprintf(Gains,"NEvents = %i\n",NEvent);
680  fprintf(Gains,"NTracks = %i\n",NTrack);
681  fprintf(Gains,"NClusters = %i\n",NCluster);
682  fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
683  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
684  Gains=fopen(OutputGains.c_str(),"w");
685  }
686  fprintf(Gains,"NEvents = %i\n",NEvent);
687  fprintf(Gains,"NTracks = %i\n",NTrack);
688  fprintf(Gains,"NClusters = %i\n",NCluster);
689  fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
690  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
691 
692  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
693  stAPVGain* APV = APVsCollOrdered[a];
694  if(APV==NULL)continue;
695 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
696  fprintf(Gains,"%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
697 
698  tree_Index = APV->Index;
699  tree_Bin = Charge_Vs_Index->getTH2F()->GetXaxis()->FindBin(APV->Index);
700  tree_DetId = APV->DetId;
701  tree_APVId = APV->APVId;
702  tree_SubDet = APV->SubDet;
703  tree_x = APV->x;
704  tree_y = APV->y;
705  tree_z = APV->z;
706  tree_Eta = APV->Eta;
707  tree_R = APV->R;
708  tree_Phi = APV->Phi;
709  tree_Thickness = APV->Thickness;
710  tree_FitMPV = APV->FitMPV;
711  tree_FitMPVErr = APV->FitMPVErr;
712  tree_FitWidth = APV->FitWidth;
713  tree_FitWidthErr= APV->FitWidthErr;
714  tree_FitChi2NDF = APV->FitChi2;
715  tree_Gain = APV->Gain;
716  tree_PrevGain = APV->PreviousGain;
717  tree_NEntries = APV->NEntries;
718  tree_isMasked = APV->isMasked;
719 
720  MyTree->Fill();
721  }
722  if(AlgoMode!="PCL")fclose(Gains);
723 
724 
725 }
726 
727 bool SiStripGainFromCalibTree::produceTagFilter(){
728 
729  // The goal of this function is to check wether or not there is enough statistics to produce a meaningful tag for the DB or not
730  if(Charge_Vs_Index->getTH2F()->Integral() < tagCondition_NClusters) {
731  cout << "[SiStripGainFromCalibTree] produceTagFilter -> Return false: Statistics is too low : " << Charge_Vs_Index->getTH2F()->Integral() << endl;
732  return false;
733  }
734  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
735  cout << "[SiStripGainFromCalibTree] produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
736  return false;
737  }
738  return true;
739 }
740 
741 SiStripApvGain* SiStripGainFromCalibTree::getNewObject()
742 {
743  SiStripApvGain* obj = new SiStripApvGain();
744  if(!harvestingMode) return obj;
745 
746  if(!produceTagFilter()){
747  cout << "[SiStripGainFromCalibTree] getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
748  setDoStore(false);
749  return obj;
750  }
751 
752 
753  std::vector<float>* theSiStripVector = NULL;
754  unsigned int PreviousDetId = 0;
755  for(unsigned int a=0;a<APVsCollOrdered.size();a++)
756  {
757  stAPVGain* APV = APVsCollOrdered[a];
758  if(APV==NULL){ printf("Bug\n"); continue; }
759  if(APV->DetId != PreviousDetId){
760  if(theSiStripVector!=NULL){
761  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
762  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
763  }
764  theSiStripVector = new std::vector<float>;
765  PreviousDetId = APV->DetId;
766  }
767  theSiStripVector->push_back(APV->Gain);
768  }
769  if(theSiStripVector!=NULL){
770  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
771  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
772  }
773 
774  return obj;
775 }
776 
777 
778 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
779 {
780 }
781 
782 void SiStripGainFromCalibTree::MakeCalibrationMap(){
783  if(!useCalibration)return;
784 
785  TChain* t1 = new TChain("SiStripCalib/APVGain");
786  t1->Add(m_calibrationPath.c_str());
787 
788  unsigned int tree_DetId;
789  unsigned char tree_APVId;
790  double tree_Gain;
791 
792  t1->SetBranchAddress("DetId" ,&tree_DetId );
793  t1->SetBranchAddress("APVId" ,&tree_APVId );
794  t1->SetBranchAddress("Gain" ,&tree_Gain );
795 
796  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
797  t1->GetEntry(ientry);
798  stAPVGain* APV = APVsColl[(tree_DetId<<3) | (unsigned int)tree_APVId];
799  APV->CalibGain = tree_Gain;
800  }
801 }
802 
803 void
804 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
805 {
806  // in AlCaHarvesting mode we just need to run the logic in the endJob step
807  if(harvestingMode) return;
808 
809  if(AlgoMode=="CalibTree")return;
810 
811  if(NEvent==0){
812  SRun = iEvent.id().run();
813  }
814  ERun = iEvent.id().run();
815  NEvent++;
816 
817  //FROM SHALLOW GAIN
818  edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
819  edm::ESHandle<SiStripGain> gainHandle; iSetup.get<SiStripGainRcd>().get(gainHandle);
820  edm::Handle<edm::View<reco::Track> > tracks; iEvent.getByLabel(theTracksLabel, tracks);
821  edm::Handle<TrajTrackAssociationCollection> associations; iEvent.getByLabel(theTracksLabel, associations);
822 
823  for( TrajTrackAssociationCollection::const_iterator association = associations->begin(); association != associations->end(); association++) {
824  const Trajectory* traj = association->key.get();
825  const reco::Track* track = association->val.get();
826 
827  //clean on the tracks
828  if(fabs(track->eta()) < MinTrackEta )continue;
829  if(fabs(track->eta()) > MaxTrackEta )continue;
830  if(track->p() < MinTrackMomentum )continue;
831  if(track->p() > MaxTrackMomentum )continue;
832  if(track->numberOfValidHits() < MinTrackHits )continue;
833  if(track->chi2()/track->ndof() > MaxTrackChiOverNdf )continue;
834  NTrack++;
835 
836  vector<TrajectoryMeasurement> measurements = traj->measurements();
837  for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
838  TrajectoryStateOnSurface trajState = measurement_it->updatedState();
839  if( !trajState.isValid() ) continue;
840 
841  const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
842  const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
843  const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
844  const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
845 
846  const SiStripCluster* Cluster = NULL;
847  uint32_t DetId = 0;
848 
849  for(unsigned int h=0;h<2;h++){
850  if(!sistripmatchedhit && h==1){
851  continue;
852  }else if(sistripmatchedhit && h==0){
853  Cluster = &sistripmatchedhit->monoCluster();
854  DetId = sistripmatchedhit->monoId();
855  }else if(sistripmatchedhit && h==1){
856  Cluster = &sistripmatchedhit->stereoCluster();;
857  DetId = sistripmatchedhit->stereoId();
858  }else if(sistripsimplehit){
859  Cluster = (sistripsimplehit->cluster()).get();
860  DetId = sistripsimplehit->geographicalId().rawId();
861  }else if(sistripsimple1dhit){
862  Cluster = (sistripsimple1dhit->cluster()).get();
863  DetId = sistripsimple1dhit->geographicalId().rawId();
864  }else{
865  continue;
866  }
867 
868  LocalVector trackDirection = trajState.localDirection();
869  double cosine = trackDirection.z()/trackDirection.mag();
870  const auto& Ampls = Cluster->amplitudes();
871  int FirstStrip = Cluster->firstStrip();
872  int APVId = FirstStrip/128;
873  bool Saturation = false;
874  bool Overlapping = false;
875  unsigned int charge = 0;
876  double PrevGain = -1;
877  stAPVGain* APV = APVsColl[(DetId<<3) | (FirstStrip/128)];
878  double Path = (10.0*APV->Thickness)/fabs(cosine);
879 
880 
881  if(gainHandle.isValid()){
882  SiStripApvGain::Range detGainRange = gainHandle->getRange(DetId);
883  PrevGain = *(detGainRange.first + APVId);
884  }
885 
886  for(unsigned int a=0;a<Ampls.size();a++){
887  charge+=Ampls[a];
888  if(Ampls[a] >=254)Saturation =true;
889  }
890 
891  if(FirstStrip==0 )Overlapping=true;
892  if(FirstStrip==128 )Overlapping=true;
893  if(FirstStrip==256 )Overlapping=true;
894  if(FirstStrip==384 )Overlapping=true;
895  if(FirstStrip==512 )Overlapping=true;
896  if(FirstStrip==640 )Overlapping=true;
897 
898  if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlapping=true;
899  if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlapping=true;
900  if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlapping=true;
901  if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlapping=true;
902  if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlapping=true;
903 
904  if(FirstStrip+Ampls.size()==127 )Overlapping=true;
905  if(FirstStrip+Ampls.size()==255 )Overlapping=true;
906  if(FirstStrip+Ampls.size()==383 )Overlapping=true;
907  if(FirstStrip+Ampls.size()==511 )Overlapping=true;
908  if(FirstStrip+Ampls.size()==639 )Overlapping=true;
909  if(FirstStrip+Ampls.size()==767 )Overlapping=true;
910 
911  //cleaning on the cluster
912  if(IsFarFromBorder(&trajState,DetId, &iSetup) == false )continue;
913  if(Overlapping == true )continue;
914  if(Saturation && !AllowSaturation )continue;
915  if(Ampls.size() > MaxNrStrips )continue;
916  NCluster++;
917 
918 
919  int Charge = 0;
920  if(useCalibration || !FirstSetOfConstants){
921  bool Saturation = false;
922  for(unsigned int s=0;s<Ampls.size();s++){
923  int StripCharge = Ampls[s];
924  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
925  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
926  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
927  if(StripCharge>1024){
928  StripCharge = 255;
929  Saturation = true;
930  }else if(StripCharge>254){
931  StripCharge = 254;
932  Saturation = true;
933  }
934  Charge += StripCharge;
935  }
936  if(Saturation && !AllowSaturation)continue;
937  }else{
938  Charge = charge;
939  }
940 
941  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
942 
943  double ClusterChargeOverPath = ( (double) Charge )/Path ;
944  if(Validation) {ClusterChargeOverPath/=PrevGain;}
945  if(OldGainRemoving){ClusterChargeOverPath*=PrevGain;}
946  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
947  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
948 
949 
950  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill(Path,Charge);
951  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill(Path,Charge);
952  }else if(APV->SubDet==StripSubdetector::TID){
953  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill(Path,Charge);
954  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill(Path,Charge);
955  }
956  }else if(APV->SubDet==StripSubdetector::TEC){
957  if(APV->Eta<0){
958  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill(Path,Charge);
959  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill(Path,Charge);
960  }
961  }else if(APV->Eta>0){
962  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill(Path,Charge);
963  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill(Path,Charge);
964  }
965  }
966  }
967 
968  }//loop on clusters
969  }//loop on measurements
970  }//loop on tracks
971 
972 }
973 
974 
975 bool SiStripGainFromCalibTree::IsFarFromBorder(TrajectoryStateOnSurface* trajState, const uint32_t detid, const edm::EventSetup* iSetup)
976 {
977  edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
978 
979  LocalPoint HitLocalPos = trajState->localPosition();
980  LocalError HitLocalError = trajState->localError().positionError() ;
981 
982  const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
983  if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
984  std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
985  return false;
986  }
987 
988  const BoundPlane plane = it->surface();
989  const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
990  const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
991 
992  double DistFromBorder = 1.0;
993  double HalfLength = it->surface().bounds().length() /2.0;
994 
995  if(trapezoidalBounds)
996  {
997  std::array<const float, 4> const & parameters = (*trapezoidalBounds).parameters();
998  //std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
999  HalfLength = parameters[3];
1000  }else if(rectangularBounds){
1001  HalfLength = it->surface().bounds().length() /2.0;
1002  }else{return false;}
1003 
1004  if (fabs(HitLocalPos.y())+HitLocalError.yy() >= (HalfLength - DistFromBorder) ) return false;
1005 
1006  return true;
1007 }
1008 
1009 
1010 
1011 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
SiStripRecHit1D::cluster
ClusterRef cluster() const
Definition: SiStripRecHit1D.h:28
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:578
DeDxTools::IsFarFromBorder
bool IsFarFromBorder(const TrajectoryStateOnSurface &trajState, const GeomDetUnit *it)
Definition: DeDxTools.cc:311
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
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
Parameters.parameters
dictionary parameters
Definition: Parameters.py:2
SiStripMatchedRecHit2D::stereoId
unsigned int stereoId() const
Definition: SiStripMatchedRecHit2D.h:23
Bounds::length
virtual float length() const =0
SiStripGainFromCalibTree::getNewObject
SiStripApvGain * getNewObject() override
Definition: SiStripGainFromCalibTree.cc:741
SiStripGainFromCalibTree::algoBeginRun
virtual void algoBeginRun(const edm::Run &run, const edm::EventSetup &iSetup) override
Definition: SiStripGainFromCalibTree.cc:238
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
SiStripMatchedRecHit2D::monoCluster
SiStripCluster const & monoCluster() const
Definition: SiStripMatchedRecHit2D.h:39
TrajectoryStateOnSurface::localDirection
LocalVector localDirection() const
Definition: TrajectoryStateOnSurface.h:122
SiStripGainFromCalibTree::SiStripGainFromCalibTree
SiStripGainFromCalibTree(const edm::ParameterSet &)
Definition: SiStripGainFromCalibTree.cc:202
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
NULL
#define NULL
Definition: scimark2.h:8
Surface::bounds
const Bounds & bounds() const
Definition: Surface.h:128
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
SiStripCluster::firstStrip
uint16_t firstStrip() const
Definition: SiStripCluster.h:46
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:106
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
LocalTrajectoryError::positionError
LocalError positionError() const
Definition: LocalTrajectoryError.h:83
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:170
SiStripApvGain::put
bool put(const uint32_t &detID, Range input)
Definition: SiStripApvGain.cc:8
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
SiStripGainFromCalibTree::APVsCollOrdered
std::vector< stAPVGain * > APVsCollOrdered
Definition: SiStripGainFromCalibTree.cc:198
Trajectory::measurements
DataContainer const & measurements() const
Definition: Trajectory.h:203
SiStripGainFromCalibTree::APVsColl
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
Definition: SiStripGainFromCalibTree.cc:199
iEvent
int iEvent
Definition: GenABIO.cc:230
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:614
cmsHarvester.path
tuple path
else: Piece not in the list, fine.
Definition: cmsHarvester.py:1591
LocalError::yy
float yy() const
Definition: LocalError.h:26
reco::TrackBase::chi2
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:512
SiStripGainFromCalibTree::IsGoodLandauFit
bool IsGoodLandauFit(double *FitResults)
Definition: SiStripGainFromCalibTree.cc:405
reco::TrackBase::ndof
double ndof() const
number of degrees of freedom of the fit
Definition: TrackBase.h:518
edm::Path
Definition: Path.h:40
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
SiStripApvGain::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripApvGain.h:27
j
int j
Definition: DBlmapReader.cc:9
Exhume::I
const std::complex< double > I
Definition: I.h:8
SiStripRecHit2D::cluster
ClusterRef cluster() const
Definition: SiStripRecHit2D.h:29
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:783
TrajectoryStateOnSurface::localError
const LocalTrajectoryError & localError() const
Definition: TrajectoryStateOnSurface.h:125
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
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:413
SiStripGainFromCalibTree::algoEndJob
virtual void algoEndJob() override
Definition: SiStripGainFromCalibTree.cc:335
DetId
Definition: DetId.h:18
SiStripGainFromCalibTree::algoAnalyze
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
Definition: SiStripGainFromCalibTree.cc:804
SiStripGainFromCalibTree::IsFarFromBorder
bool IsFarFromBorder(TrajectoryStateOnSurface *trajState, const uint32_t detid, const edm::EventSetup *iSetup)
Definition: SiStripGainFromCalibTree.cc:975
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
SiStripGainFromCalibTree::getPeakOfLandau
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
Definition: SiStripGainFromCalibTree.cc:380
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
ConditionDBWriter< SiStripApvGain >::setDoStore
void setDoStore(const bool doStore)
When set to false the payload will not be written to the db.
Definition: ConditionDBWriter.h:372
SiStripMatchedRecHit2D::stereoCluster
SiStripCluster const & stereoCluster() const
Definition: SiStripMatchedRecHit2D.h:36
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:60
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
SiStripMatchedRecHit2D::monoId
unsigned int monoId() const
Definition: SiStripMatchedRecHit2D.h:24
TrackingRecHit::geographicalId
DetId geographicalId() const
Definition: TrackingRecHit.h:116
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
CSCOverlapsAlignmentAlgorithm_cff.Path
tuple Path
Definition: CSCOverlapsAlignmentAlgorithm_cff.py:25
MonitorElement::getTH2F
TH2F * getTH2F(void) const
Definition: MonitorElement.cc:1481
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:47
T
long double T
Definition: Basic3DVectorLD.h:57
SiStripCluster::amplitudes
const std::vector< uint8_t > & amplitudes() const
Definition: SiStripCluster.h:59
edm::Run
Definition: Run.h:41
AlcaSiStripGainsHarvester_cff.DQMStore
tuple DQMStore
Definition: AlcaSiStripGainsHarvester_cff.py:10
Generated for CMSSW Reference Manual by doxygen 1.8.5