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") storeOnTree();
363 }
364 
365 
366 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
367 {
368  FitResults[0] = -0.5; //MPV
369  FitResults[1] = 0; //MPV error
370  FitResults[2] = -0.5; //Width
371  FitResults[3] = 0; //Width error
372  FitResults[4] = -0.5; //Fit Chi2/NDF
373 
374  if( InputHisto->GetEntries() < MinNrEntries)return;
375 
376  // perform fit with standard landau
377  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
378  MyLandau->SetParameter(1,300);
379  InputHisto->Fit(MyLandau,"0QR WW");
380 
381  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
382  FitResults[0] = MyLandau->GetParameter(1); //MPV
383  FitResults[1] = MyLandau->GetParError(1); //MPV error
384  FitResults[2] = MyLandau->GetParameter(2); //Width
385  FitResults[3] = MyLandau->GetParError(2); //Width error
386  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
387 
388  delete MyLandau;
389 }
390 
391 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
392  if(FitResults[0] <= 0 )return false;
393 // if(FitResults[1] > MaxMPVError )return false;
394 // if(FitResults[4] > MaxChi2OverNDF)return false;
395  return true;
396 }
397 
398 
399 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
400 {
401  for(unsigned int i=0;i<VInputFiles.size();i++){
402  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
403  TChain* tree = new TChain("commonCalibrationTree/tree");
404  tree->Add(VInputFiles[i].c_str());
405 
406  TString EventPrefix("");
407  TString EventSuffix("");
408 
409  TString TrackPrefix("track");
410  TString TrackSuffix("");
411 
412  TString CalibPrefix("GainCalibration");
413  TString CalibSuffix("");
414 
415  unsigned int eventnumber = 0; tree->SetBranchAddress(EventPrefix + "event" + EventSuffix, &eventnumber , NULL);
416  unsigned int runnumber = 0; tree->SetBranchAddress(EventPrefix + "run" + EventSuffix, &runnumber , NULL);
417  std::vector<bool>* TrigTech = 0; tree->SetBranchAddress(EventPrefix + "TrigTech" + EventSuffix, &TrigTech , NULL);
418 
419  std::vector<double>* trackchi2ndof = 0; tree->SetBranchAddress(TrackPrefix + "chi2ndof" + TrackSuffix, &trackchi2ndof , NULL);
420  std::vector<float>* trackp = 0; tree->SetBranchAddress(TrackPrefix + "momentum" + TrackSuffix, &trackp , NULL);
421  std::vector<float>* trackpt = 0; tree->SetBranchAddress(TrackPrefix + "pt" + TrackSuffix, &trackpt , NULL);
422  std::vector<double>* tracketa = 0; tree->SetBranchAddress(TrackPrefix + "eta" + TrackSuffix, &tracketa , NULL);
423  std::vector<double>* trackphi = 0; tree->SetBranchAddress(TrackPrefix + "phi" + TrackSuffix, &trackphi , NULL);
424  std::vector<unsigned int>* trackhitsvalid = 0; tree->SetBranchAddress(TrackPrefix + "hitsvalid" + TrackSuffix, &trackhitsvalid, NULL);
425 
426  std::vector<int>* trackindex = 0; tree->SetBranchAddress(CalibPrefix + "trackindex" + CalibSuffix, &trackindex , NULL);
427  std::vector<unsigned int>* rawid = 0; tree->SetBranchAddress(CalibPrefix + "rawid" + CalibSuffix, &rawid , NULL);
428  std::vector<float>* localdirx = 0; tree->SetBranchAddress(CalibPrefix + "localdirx" + CalibSuffix, &localdirx , NULL);
429  std::vector<float>* localdiry = 0; tree->SetBranchAddress(CalibPrefix + "localdiry" + CalibSuffix, &localdiry , NULL);
430  std::vector<float>* localdirz = 0; tree->SetBranchAddress(CalibPrefix + "localdirz" + CalibSuffix, &localdirz , NULL);
431  std::vector<unsigned short>* firststrip = 0; tree->SetBranchAddress(CalibPrefix + "firststrip" + CalibSuffix, &firststrip , NULL);
432  std::vector<unsigned short>* nstrips = 0; tree->SetBranchAddress(CalibPrefix + "nstrips" + CalibSuffix, &nstrips , NULL);
433  std::vector<bool>* saturation = 0; tree->SetBranchAddress(CalibPrefix + "saturation" + CalibSuffix, &saturation , NULL);
434  std::vector<bool>* overlapping = 0; tree->SetBranchAddress(CalibPrefix + "overlapping" + CalibSuffix, &overlapping , NULL);
435  std::vector<bool>* farfromedge = 0; tree->SetBranchAddress(CalibPrefix + "farfromedge" + CalibSuffix, &farfromedge , NULL);
436  std::vector<unsigned int>* charge = 0; tree->SetBranchAddress(CalibPrefix + "charge" + CalibSuffix, &charge , NULL);
437  std::vector<float>* path = 0; tree->SetBranchAddress(CalibPrefix + "path" + CalibSuffix, &path , NULL);
438  std::vector<float>* chargeoverpath = 0; tree->SetBranchAddress(CalibPrefix + "chargeoverpath" + CalibSuffix, &chargeoverpath, NULL);
439  std::vector<unsigned char>* amplitude = 0; tree->SetBranchAddress(CalibPrefix + "amplitude" + CalibSuffix, &amplitude , NULL);
440  std::vector<double>* gainused = 0; tree->SetBranchAddress(CalibPrefix + "gainused" + CalibSuffix, &gainused , NULL);
441 
442 
443  printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));
444  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
445  printf("Looping on the Tree :");
446  int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
447  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
448  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
449  tree->GetEntry(ientry);
450 
451  if(runnumber<SRun)SRun=runnumber;
452  if(runnumber>ERun)ERun=runnumber;
453 
454  NEvent++;
455  NTrack+=(*trackp).size();
456 
457  unsigned int FirstAmplitude=0;
458  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
459  FirstAmplitude+=(*nstrips)[i];
460  int TI = (*trackindex)[i];
461  if((*tracketa )[TI] < MinTrackEta )continue;
462  if((*tracketa )[TI] > MaxTrackEta )continue;
463  if((*trackp )[TI] < MinTrackMomentum )continue;
464  if((*trackp )[TI] > MaxTrackMomentum )continue;
465  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
466  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
467  if((*farfromedge)[i] == false )continue;
468  if((*overlapping)[i] == true )continue;
469  if((*saturation )[i] && !AllowSaturation)continue;
470  if((*nstrips )[i] > MaxNrStrips )continue;
471  NCluster++;
472 
473  stAPVGain* APV = APVsColl[((*rawid)[i]<<3) | ((*firststrip)[i]/128)];
474  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
475 
476  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
477  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
478  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
479 
480  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
481  //printf("Amplitudes: ");
482  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
483  //printf("\n");
484 
485 
486  int Charge = 0;
487  if(useCalibration || !FirstSetOfConstants){
488  bool Saturation = false;
489  for(unsigned int s=0;s<(*nstrips)[i];s++){
490  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
491  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
492  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
493  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
494  if(StripCharge>1024){
495  StripCharge = 255;
496  Saturation = true;
497  }else if(StripCharge>254){
498  StripCharge = 254;
499  Saturation = true;
500  }
501  Charge += StripCharge;
502  }
503  if(Saturation && !AllowSaturation)continue;
504  }else{
505  Charge = (*charge)[i];
506  }
507 
508  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
509 
510  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
511  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
512  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
513  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
514  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
515 
516 
517  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill((*path)[i],Charge);
518  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill((*path)[i],Charge);
519  }else if(APV->SubDet==StripSubdetector::TID){
520  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
521  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
522  }
523  }else if(APV->SubDet==StripSubdetector::TEC){
524  if(APV->Eta<0){
525  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
526  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
527  }
528  }else if(APV->Eta>0){
529  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
530  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
531  }
532  }
533  }
534 
535  }// END OF ON-CLUSTER LOOP
536  }printf("\n");// END OF EVENT LOOP
537 
538  }
539 }
540 
541 
542 
543 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
544  unsigned int I=0;
545  TH1F* Proj = NULL;
546  double FitResults[5];
547  double MPVmean = 300;
548 
549  TH2F *chvsidx = Charge_Vs_Index->getTH2F();
550 
551 
552  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
553  printf("Fitting Charge Distribution :");
554  int TreeStep = APVsColl.size()/50;
555  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
556  if(I%TreeStep==0){printf(".");fflush(stdout);}
557  stAPVGain* APV = it->second;
558  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
559 
560  if(APV->isMasked){MASKED++; continue;}
561 
562  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
563  if(!Proj)continue;
564 
565  if(CalibrationLevel==0){
566  }else if(CalibrationLevel==1){
567  int SecondAPVId = APV->APVId;
568  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
569  stAPVGain* APV2 = APVsColl[(APV->DetId<<3) | SecondAPVId];
570  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
571  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
572  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
573  }else if(CalibrationLevel==2){
574  for(unsigned int i=0;i<6;i++){
575  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
576  tmpit = APVsColl.find((APV->DetId<<3) | i);
577  if(tmpit==APVsColl.end())continue;
578  stAPVGain* APV2 = tmpit->second;
579  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
580  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
581  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
582  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
583  }
584  }else{
585  CalibrationLevel = 0;
586  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
587  }
588 
589  getPeakOfLandau(Proj,FitResults);
590  APV->FitMPV = FitResults[0];
591  APV->FitMPVErr = FitResults[1];
592  APV->FitWidth = FitResults[2];
593  APV->FitWidthErr = FitResults[3];
594  APV->FitChi2 = FitResults[4];
595  APV->NEntries = Proj->GetEntries();
596 
597  if(IsGoodLandauFit(FitResults)){
598  APV->Gain = APV->FitMPV / MPVmean;
599  GOOD++;
600  }else{
601  APV->Gain = APV->PreviousGain;
602  BAD++;
603  }
604  if(APV->Gain<=0) APV->Gain = 1;
605 
606  //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);
607  delete Proj;
608  }printf("\n");
609 }
610 
611 
612 void SiStripGainFromCalibTree::storeOnTree()
613 {
614  tfs = edm::Service<TFileService>().operator->();
615 
616  unsigned int tree_Index;
617  unsigned int tree_Bin;
618  unsigned int tree_DetId;
619  unsigned char tree_APVId;
620  unsigned char tree_SubDet;
621  float tree_x;
622  float tree_y;
623  float tree_z;
624  float tree_Eta;
625  float tree_R;
626  float tree_Phi;
627  float tree_Thickness;
628  float tree_FitMPV;
629  float tree_FitMPVErr;
630  float tree_FitWidth;
631  float tree_FitWidthErr;
632  float tree_FitChi2NDF;
633  double tree_Gain;
634  double tree_PrevGain;
635  double tree_NEntries;
636  bool tree_isMasked;
637 
638  TTree* MyTree;
639  MyTree = tfs->make<TTree> ("APVGain","APVGain");
640  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
641  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
642  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
643  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
644  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
645  MyTree->Branch("x" ,&tree_x ,"x/F");
646  MyTree->Branch("y" ,&tree_y ,"y/F");
647  MyTree->Branch("z" ,&tree_z ,"z/F");
648  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
649  MyTree->Branch("R" ,&tree_R ,"R/F");
650  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
651  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
652  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
653  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
654  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
655  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
656  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
657  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
658  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
659  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
660  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
661 
662 
663  FILE* Gains = stdout;
664  if(AlgoMode!="PCL"){
665  fprintf(Gains,"NEvents = %i\n",NEvent);
666  fprintf(Gains,"NTracks = %i\n",NTrack);
667  fprintf(Gains,"NClusters = %i\n",NCluster);
668  fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
669  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
670  Gains=fopen(OutputGains.c_str(),"w");
671  }
672  fprintf(Gains,"NEvents = %i\n",NEvent);
673  fprintf(Gains,"NTracks = %i\n",NTrack);
674  fprintf(Gains,"NClusters = %i\n",NCluster);
675  fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
676  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
677 
678  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
679  stAPVGain* APV = APVsCollOrdered[a];
680  if(APV==NULL)continue;
681 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
682  fprintf(Gains,"%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
683 
684  tree_Index = APV->Index;
685  tree_Bin = Charge_Vs_Index->getTH2F()->GetXaxis()->FindBin(APV->Index);
686  tree_DetId = APV->DetId;
687  tree_APVId = APV->APVId;
688  tree_SubDet = APV->SubDet;
689  tree_x = APV->x;
690  tree_y = APV->y;
691  tree_z = APV->z;
692  tree_Eta = APV->Eta;
693  tree_R = APV->R;
694  tree_Phi = APV->Phi;
695  tree_Thickness = APV->Thickness;
696  tree_FitMPV = APV->FitMPV;
697  tree_FitMPVErr = APV->FitMPVErr;
698  tree_FitWidth = APV->FitWidth;
699  tree_FitWidthErr= APV->FitWidthErr;
700  tree_FitChi2NDF = APV->FitChi2;
701  tree_Gain = APV->Gain;
702  tree_PrevGain = APV->PreviousGain;
703  tree_NEntries = APV->NEntries;
704  tree_isMasked = APV->isMasked;
705 
706  MyTree->Fill();
707  }
708  if(AlgoMode!="PCL")fclose(Gains);
709 
710 
711 }
712 
713 bool SiStripGainFromCalibTree::produceTagFilter(){
714 
715  // 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
716  if(Charge_Vs_Index->getTH2F()->Integral() < tagCondition_NClusters) {
717  cout << "[SiStripGainFromCalibTree] produceTagFilter -> Return false: Statistics is too low : " << Charge_Vs_Index->getTH2F()->Integral() << endl;
718  return false;
719  }
720  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
721  cout << "[SiStripGainFromCalibTree] produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
722  return false;
723  }
724  return true;
725 }
726 
727 SiStripApvGain* SiStripGainFromCalibTree::getNewObject()
728 {
729  SiStripApvGain* obj = new SiStripApvGain();
730  if(!harvestingMode) return obj;
731 
732  if(!produceTagFilter()){
733  cout << "[SiStripGainFromCalibTree] getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
734  setDoStore(false);
735  return obj;
736  }
737 
738 
739  std::vector<float>* theSiStripVector = NULL;
740  unsigned int PreviousDetId = 0;
741  for(unsigned int a=0;a<APVsCollOrdered.size();a++)
742  {
743  stAPVGain* APV = APVsCollOrdered[a];
744  if(APV==NULL){ printf("Bug\n"); continue; }
745  if(APV->DetId != PreviousDetId){
746  if(theSiStripVector!=NULL){
747  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
748  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
749  }
750  theSiStripVector = new std::vector<float>;
751  PreviousDetId = APV->DetId;
752  }
753  theSiStripVector->push_back(APV->Gain);
754  }
755  if(theSiStripVector!=NULL){
756  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
757  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
758  }
759 
760  return obj;
761 }
762 
763 
764 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
765 {
766 }
767 
768 void SiStripGainFromCalibTree::MakeCalibrationMap(){
769  if(!useCalibration)return;
770 
771  TChain* t1 = new TChain("SiStripCalib/APVGain");
772  t1->Add(m_calibrationPath.c_str());
773 
774  unsigned int tree_DetId;
775  unsigned char tree_APVId;
776  double tree_Gain;
777 
778  t1->SetBranchAddress("DetId" ,&tree_DetId );
779  t1->SetBranchAddress("APVId" ,&tree_APVId );
780  t1->SetBranchAddress("Gain" ,&tree_Gain );
781 
782  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
783  t1->GetEntry(ientry);
784  stAPVGain* APV = APVsColl[(tree_DetId<<3) | (unsigned int)tree_APVId];
785  APV->CalibGain = tree_Gain;
786  }
787 }
788 
789 void
790 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
791 {
792  // in AlCaHarvesting mode we just need to run the logic in the endJob step
793  if(harvestingMode) return;
794 
795  if(AlgoMode=="CalibTree")return;
796 
797  if(NEvent==0){
798  SRun = iEvent.id().run();
799  }
800  ERun = iEvent.id().run();
801  NEvent++;
802 
803  //FROM SHALLOW GAIN
804  edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
805  edm::ESHandle<SiStripGain> gainHandle; iSetup.get<SiStripGainRcd>().get(gainHandle);
806  edm::Handle<edm::View<reco::Track> > tracks; iEvent.getByLabel(theTracksLabel, tracks);
807  edm::Handle<TrajTrackAssociationCollection> associations; iEvent.getByLabel(theTracksLabel, associations);
808 
809  for( TrajTrackAssociationCollection::const_iterator association = associations->begin(); association != associations->end(); association++) {
810  const Trajectory* traj = association->key.get();
811  const reco::Track* track = association->val.get();
812 
813  //clean on the tracks
814  if(fabs(track->eta()) < MinTrackEta )continue;
815  if(fabs(track->eta()) > MaxTrackEta )continue;
816  if(track->p() < MinTrackMomentum )continue;
817  if(track->p() > MaxTrackMomentum )continue;
818  if(track->numberOfValidHits() < MinTrackHits )continue;
819  if(track->chi2()/track->ndof() > MaxTrackChiOverNdf )continue;
820  NTrack++;
821 
822  vector<TrajectoryMeasurement> measurements = traj->measurements();
823  for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
824  TrajectoryStateOnSurface trajState = measurement_it->updatedState();
825  if( !trajState.isValid() ) continue;
826 
827  const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
828  const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
829  const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
830  const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
831 
832  const SiStripCluster* Cluster = NULL;
833  uint32_t DetId = 0;
834 
835  for(unsigned int h=0;h<2;h++){
836  if(!sistripmatchedhit && h==1){
837  continue;
838  }else if(sistripmatchedhit && h==0){
839  Cluster = &sistripmatchedhit->monoCluster();
840  DetId = sistripmatchedhit->monoId();
841  }else if(sistripmatchedhit && h==1){
842  Cluster = &sistripmatchedhit->stereoCluster();;
843  DetId = sistripmatchedhit->stereoId();
844  }else if(sistripsimplehit){
845  Cluster = (sistripsimplehit->cluster()).get();
846  DetId = sistripsimplehit->geographicalId().rawId();
847  }else if(sistripsimple1dhit){
848  Cluster = (sistripsimple1dhit->cluster()).get();
849  DetId = sistripsimple1dhit->geographicalId().rawId();
850  }else{
851  continue;
852  }
853 
854  LocalVector trackDirection = trajState.localDirection();
855  double cosine = trackDirection.z()/trackDirection.mag();
856  const vector<uint8_t>& Ampls = Cluster->amplitudes();
857  int FirstStrip = Cluster->firstStrip();
858  int APVId = FirstStrip/128;
859  bool Saturation = false;
860  bool Overlapping = false;
861  unsigned int charge = 0;
862  double PrevGain = -1;
863  stAPVGain* APV = APVsColl[(DetId<<3) | (FirstStrip/128)];
864  double Path = (10.0*APV->Thickness)/fabs(cosine);
865 
866 
867  if(gainHandle.isValid()){
868  SiStripApvGain::Range detGainRange = gainHandle->getRange(DetId);
869  PrevGain = *(detGainRange.first + APVId);
870  }
871 
872  for(unsigned int a=0;a<Ampls.size();a++){
873  charge+=Ampls[a];
874  if(Ampls[a] >=254)Saturation =true;
875  }
876 
877  if(FirstStrip==0 )Overlapping=true;
878  if(FirstStrip==128 )Overlapping=true;
879  if(FirstStrip==256 )Overlapping=true;
880  if(FirstStrip==384 )Overlapping=true;
881  if(FirstStrip==512 )Overlapping=true;
882  if(FirstStrip==640 )Overlapping=true;
883 
884  if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlapping=true;
885  if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlapping=true;
886  if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlapping=true;
887  if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlapping=true;
888  if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlapping=true;
889 
890  if(FirstStrip+Ampls.size()==127 )Overlapping=true;
891  if(FirstStrip+Ampls.size()==255 )Overlapping=true;
892  if(FirstStrip+Ampls.size()==383 )Overlapping=true;
893  if(FirstStrip+Ampls.size()==511 )Overlapping=true;
894  if(FirstStrip+Ampls.size()==639 )Overlapping=true;
895  if(FirstStrip+Ampls.size()==767 )Overlapping=true;
896 
897  //cleaning on the cluster
898  if(IsFarFromBorder(&trajState,DetId, &iSetup) == false )continue;
899  if(Overlapping == true )continue;
900  if(Saturation && !AllowSaturation )continue;
901  if(Ampls.size() > MaxNrStrips )continue;
902  NCluster++;
903 
904 
905  int Charge = 0;
906  if(useCalibration || !FirstSetOfConstants){
907  bool Saturation = false;
908  for(unsigned int s=0;s<Ampls.size();s++){
909  int StripCharge = Ampls[s];
910  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
911  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
912  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
913  if(StripCharge>1024){
914  StripCharge = 255;
915  Saturation = true;
916  }else if(StripCharge>254){
917  StripCharge = 254;
918  Saturation = true;
919  }
920  Charge += StripCharge;
921  }
922  if(Saturation && !AllowSaturation)continue;
923  }else{
924  Charge = charge;
925  }
926 
927  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
928 
929  double ClusterChargeOverPath = ( (double) Charge )/Path ;
930  if(Validation) {ClusterChargeOverPath/=PrevGain;}
931  if(OldGainRemoving){ClusterChargeOverPath*=PrevGain;}
932  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
933  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
934 
935 
936  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill(Path,Charge);
937  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill(Path,Charge);
938  }else if(APV->SubDet==StripSubdetector::TID){
939  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill(Path,Charge);
940  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill(Path,Charge);
941  }
942  }else if(APV->SubDet==StripSubdetector::TEC){
943  if(APV->Eta<0){
944  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill(Path,Charge);
945  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill(Path,Charge);
946  }
947  }else if(APV->Eta>0){
948  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill(Path,Charge);
949  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill(Path,Charge);
950  }
951  }
952  }
953 
954  }//loop on clusters
955  }//loop on measurements
956  }//loop on tracks
957 
958 }
959 
960 
961 bool SiStripGainFromCalibTree::IsFarFromBorder(TrajectoryStateOnSurface* trajState, const uint32_t detid, const edm::EventSetup* iSetup)
962 {
963  edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
964 
965  LocalPoint HitLocalPos = trajState->localPosition();
966  LocalError HitLocalError = trajState->localError().positionError() ;
967 
968  const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
969  if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
970  std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
971  return false;
972  }
973 
974  const BoundPlane plane = it->surface();
975  const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
976  const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
977 
978  double DistFromBorder = 1.0;
979  double HalfLength = it->surface().bounds().length() /2.0;
980 
981  if(trapezoidalBounds)
982  {
983  std::array<const float, 4> const & parameters = (*trapezoidalBounds).parameters();
984  //std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
985  HalfLength = parameters[3];
986  }else if(rectangularBounds){
987  HalfLength = it->surface().bounds().length() /2.0;
988  }else{return false;}
989 
990  if (fabs(HitLocalPos.y())+HitLocalError.yy() >= (HalfLength - DistFromBorder) ) return false;
991 
992  return true;
993 }
994 
995 
996 
997 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:591
DeDxTools::IsFarFromBorder
bool IsFarFromBorder(const TrajectoryStateOnSurface &trajState, const GeomDetUnit *it)
Definition: DeDxTools.cc:302
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
Bounds::length
virtual float length() const =0
SiStripGainFromCalibTree::getNewObject
SiStripApvGain * getNewObject() override
Definition: SiStripGainFromCalibTree.cc:727
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
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:35
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
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:164
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:627
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:525
SiStripGainFromCalibTree::IsGoodLandauFit
bool IsGoodLandauFit(double *FitResults)
Definition: SiStripGainFromCalibTree.cc:391
reco::TrackBase::ndof
double ndof() const
number of degrees of freedom of the fit
Definition: TrackBase.h:531
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
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:796
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:402
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:790
SiStripGainFromCalibTree::IsFarFromBorder
bool IsFarFromBorder(TrajectoryStateOnSurface *trajState, const uint32_t detid, const edm::EventSetup *iSetup)
Definition: SiStripGainFromCalibTree.cc:961
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:366
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
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
a
double a
Definition: hdecay.h:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
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:1551
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:48
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