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/SiStripRecHit2D.h"
37 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
38 #include "DataFormats/SiStripDetId/interface/SiStripSubStructure.h"
39 #include "DataFormats/DetId/interface/DetId.h"
40 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
41 #include "DataFormats/TrackReco/interface/DeDxHit.h"
42 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
43 
44 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
45 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
46 
47 #include "TrackingTools/PatternTools/interface/Trajectory.h"
48 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
49 
50 #include "DQMServices/Core/interface/DQMStore.h"
51 #include "DQMServices/Core/interface/MonitorElement.h"
52 #include "FWCore/ServiceRegistry/interface/Service.h"
53 #include "CommonTools/UtilAlgos/interface/TFileService.h"
54 
55 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
56 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
57 
58 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
59 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
60 
61 
62 #include "TFile.h"
63 #include "TObjString.h"
64 #include "TString.h"
65 #include "TH1F.h"
66 #include "TH2F.h"
67 #include "TProfile.h"
68 #include "TF1.h"
69 #include "TROOT.h"
70 #include "TTree.h"
71 #include "TChain.h"
72 
73 #include <ext/hash_map>
74 
75 
76 
77 using namespace edm;
78 using namespace reco;
79 using namespace std;
80 using namespace __gnu_cxx;
81 using __gnu_cxx::hash_map;
82 using __gnu_cxx::hash;
83 
84 struct stAPVGain{
85  unsigned int Index;
86  unsigned int Bin;
87  unsigned int DetId;
88  unsigned int APVId;
89  unsigned int SubDet;
90  float x;
91  float y;
92  float z;
93  float Eta;
94  float R;
95  float Phi;
96  float Thickness;
97  double FitMPV;
98  double FitMPVErr;
99  double FitWidth;
100  double FitWidthErr;
101  double FitChi2;
102  double Gain;
103  double CalibGain;
104  double PreviousGain;
105  double NEntries;
106  TH1F* HCharge;
107  TH1F* HChargeN;
108  bool isMasked;
109 };
110 
111 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
112  public:
113  explicit SiStripGainFromCalibTree(const edm::ParameterSet&);
114  ~SiStripGainFromCalibTree();
115 
116 
117  private:
118  virtual void algoBeginJob (const edm::EventSetup&) override;
119  virtual void algoEndJob () override;
120  virtual void algoAnalyze (const edm::Event &, const edm::EventSetup &) override;
121 
122  void algoAnalyzeTheTree();
123  void algoComputeMPVandGain();
124 
125  void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
126  bool IsGoodLandauFit(double* FitResults);
127  void storeOnTree();
128 
129  void MakeCalibrationMap();
130 
131 
132  SiStripApvGain* getNewObject() override;
133  edm::Service<TFileService> tfs;
134 
135  double MinNrEntries;
136  double MaxMPVError;
137  double MaxChi2OverNDF;
138  double MinTrackMomentum;
139  double MaxTrackMomentum;
140  double MinTrackEta;
141  double MaxTrackEta;
142  unsigned int MaxNrStrips;
143  unsigned int MinTrackHits;
144  double MaxTrackChiOverNdf;
145  bool AllowSaturation;
146  bool FirstSetOfConstants;
147  bool Validation;
148  bool OldGainRemoving;
149  int CalibrationLevel;
150 
151  bool useCalibration;
152  string m_calibrationPath;
153 
154  std::string OutputGains;
155  vector<string> VInputFiles;
156 
157  TH2F* Charge_Vs_Index;
158  TH2F* Charge_Vs_Index_Absolute;
159  TH2F* Charge_Vs_PathlengthTIB;
160  TH2F* Charge_Vs_PathlengthTOB;
161  TH2F* Charge_Vs_PathlengthTIDP;
162  TH2F* Charge_Vs_PathlengthTIDM;
163  TH2F* Charge_Vs_PathlengthTECP1;
164  TH2F* Charge_Vs_PathlengthTECP2;
165  TH2F* Charge_Vs_PathlengthTECM1;
166  TH2F* Charge_Vs_PathlengthTECM2;
167 
168  unsigned int NEvent;
169  unsigned int NTrack;
170  unsigned int NCluster;
171  unsigned int SRun;
172  unsigned int ERun;
173  unsigned int GOOD;
174  unsigned int BAD;
175 
176  private :
177  class isEqual{
178  public:
179  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
180  };
181 
182  std::vector<stAPVGain*> APVsCollOrdered;
183  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
184 };
185 
186 SiStripGainFromCalibTree::SiStripGainFromCalibTree(const edm::ParameterSet& iConfig) : ConditionDBWriter<SiStripApvGain>(iConfig)
187 {
188  OutputGains = iConfig.getParameter<std::string>("OutputGains");
189 
190  MinNrEntries = iConfig.getUntrackedParameter<double> ("minNrEntries" , 20);
191  MaxMPVError = iConfig.getUntrackedParameter<double> ("maxMPVError" , 500.0);
192  MaxChi2OverNDF = iConfig.getUntrackedParameter<double> ("maxChi2OverNDF" , 5.0);
193  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 3.0);
194  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
195  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
196  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
197  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 2);
198  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 8);
199  MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double> ("MaxTrackChiOverNdf" , 3);
200  AllowSaturation = iConfig.getUntrackedParameter<bool> ("AllowSaturation" , false);
201  FirstSetOfConstants = iConfig.getUntrackedParameter<bool> ("FirstSetOfConstants", true);
202  Validation = iConfig.getUntrackedParameter<bool> ("Validation" , false);
203  OldGainRemoving = iConfig.getUntrackedParameter<bool> ("OldGainRemoving" , false);
204 
205  CalibrationLevel = iConfig.getUntrackedParameter<int> ("CalibrationLevel" , 0);
206  VInputFiles = iConfig.getParameter<vector<string> > ("InputFiles");
207 
208 
209  useCalibration = iConfig.getUntrackedParameter<bool>("UseCalibration", false);
210  m_calibrationPath = iConfig.getUntrackedParameter<string>("calibrationPath");
211 
212 
213 }
214 
215 
216 
217 
218 void
219 SiStripGainFromCalibTree::algoBeginJob(const edm::EventSetup& iSetup)
220 {
221  Charge_Vs_Index = tfs->make<TH2F>("Charge_Vs_Index" , "Charge_Vs_Index" , 72785, 0 , 72784,1000,0,2000);
222  Charge_Vs_Index_Absolute = tfs->make<TH2F>("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 72785, 0 , 72784, 500,0,2000);
223  Charge_Vs_PathlengthTIB = tfs->make<TH2F>("Charge_Vs_PathlengthTIB" , "Charge_Vs_PathlengthTIB" , 20 , 0.3 , 1.3 , 250,0,2000);
224  Charge_Vs_PathlengthTOB = tfs->make<TH2F>("Charge_Vs_PathlengthTOB" , "Charge_Vs_PathlengthTOB" , 20 , 0.3 , 1.3 , 250,0,2000);
225  Charge_Vs_PathlengthTIDP = tfs->make<TH2F>("Charge_Vs_PathlengthTIDP" , "Charge_Vs_PathlengthTIDP" , 20 , 0.3 , 1.3 , 250,0,2000);
226  Charge_Vs_PathlengthTIDM = tfs->make<TH2F>("Charge_Vs_PathlengthTIDM" , "Charge_Vs_PathlengthTIDM" , 20 , 0.3 , 1.3 , 250,0,2000);
227  Charge_Vs_PathlengthTECP1 = tfs->make<TH2F>("Charge_Vs_PathlengthTECP1", "Charge_Vs_PathlengthTECP1", 20 , 0.3 , 1.3 , 250,0,2000);
228  Charge_Vs_PathlengthTECP2 = tfs->make<TH2F>("Charge_Vs_PathlengthTECP2", "Charge_Vs_PathlengthTECP2", 20 , 0.3 , 1.3 , 250,0,2000);
229  Charge_Vs_PathlengthTECM1 = tfs->make<TH2F>("Charge_Vs_PathlengthTECM1", "Charge_Vs_PathlengthTECM1", 20 , 0.3 , 1.3 , 250,0,2000);
230  Charge_Vs_PathlengthTECM2 = tfs->make<TH2F>("Charge_Vs_PathlengthTECM2", "Charge_Vs_PathlengthTECM2", 20 , 0.3 , 1.3 , 250,0,2000);
231 
232  edm::ESHandle<TrackerGeometry> tkGeom;
233  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
234  vector<GeomDet*> Det = tkGeom->dets();
235 
236  edm::ESHandle<SiStripGain> gainHandle;
237  iSetup.get<SiStripGainRcd>().get(gainHandle);
238  if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}
239 
240 // size_t numberOfTag = gainHandle->getNumberOfTags();
241  for(unsigned int i=0;i<gainHandle->getNumberOfTags();i++){
242  printf("Reccord %i --> Rcd Name = %s Label Name = %s\n",i,gainHandle->getRcdName(i).c_str(), gainHandle->getLabelName(i).c_str());
243  }
244 
245 
246  edm::ESHandle<SiStripQuality> SiStripQuality_;
247  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
248 
249 
250  unsigned int Index=0;
251  for(unsigned int i=0;i<Det.size();i++){
252  DetId Detid = Det[i]->geographicalId();
253  int SubDet = Detid.subdetId();
254 
255  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
256  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
257 
258  StripGeomDetUnit* DetUnit = dynamic_cast<StripGeomDetUnit*> (Det[i]);
259  if(!DetUnit)continue;
260 
261  const StripTopology& Topo = DetUnit->specificTopology();
262  unsigned int NAPV = Topo.nstrips()/128;
263  for(unsigned int j=0;j<NAPV;j++){
264  stAPVGain* APV = new stAPVGain;
265  APV->Index = Index;
266  APV->Bin = Charge_Vs_Index->GetXaxis()->FindBin(APV->Index);
267  APV->DetId = Detid.rawId();
268  APV->APVId = j;
269  APV->SubDet = SubDet;
270  APV->FitMPV = -1;
271  APV->FitMPVErr = -1;
272  APV->FitWidth = -1;
273  APV->FitWidthErr = -1;
274  APV->FitChi2 = -1;
275  APV->Gain = -1;
276  APV->PreviousGain = 1;
277  APV->x = DetUnit->position().basicVector().x();
278  APV->y = DetUnit->position().basicVector().y();
279  APV->z = DetUnit->position().basicVector().z();
280  APV->Eta = DetUnit->position().basicVector().eta();
281  APV->Phi = DetUnit->position().basicVector().phi();
282  APV->R = DetUnit->position().basicVector().transverse();
283  APV->Thickness = DetUnit->surface().bounds().thickness();
284  APV->NEntries = 0;
285  APV->isMasked = SiStripQuality_->IsApvBad(Detid.rawId(),j);
286 
287  if(!FirstSetOfConstants){
288  if(gainHandle->getNumberOfTags()!=2){printf("ERROR: NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n");fflush(stdout);exit(0);};
289  APV->PreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
290  printf("DETID = %7i APVID=%1i Previous Gain=%8.4f\n",APV->DetId,APV->APVId,APV->PreviousGain);
291  }
292 
293  APVsCollOrdered.push_back(APV);
294  APVsColl[(APV->DetId<<3) | APV->APVId] = APV;
295  Index++;
296  }
297  }
298  }
299 
300  MakeCalibrationMap();
301 
302 
303  NEvent = 0;
304  NTrack = 0;
305  NCluster = 0;
306  SRun = 1<<31;
307  ERun = 0;
308  GOOD = 0;
309  BAD = 0;
310 
311  algoAnalyzeTheTree();
312  algoComputeMPVandGain();
313 }
314 
315 
316 void
317 SiStripGainFromCalibTree::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
318 {
319 }
320 
321 void
322 SiStripGainFromCalibTree::algoEndJob() {
323 }
324 
325 
326 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
327 {
328  FitResults[0] = -0.5; //MPV
329  FitResults[1] = 0; //MPV error
330  FitResults[2] = -0.5; //Width
331  FitResults[3] = 0; //Width error
332  FitResults[4] = -0.5; //Fit Chi2/NDF
333 
334  if( InputHisto->GetEntries() < MinNrEntries)return;
335 
336  // perform fit with standard landau
337  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
338  MyLandau->SetParameter(1,300);
339  InputHisto->Fit("MyLandau","0QR WW");
340 
341  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
342  FitResults[0] = MyLandau->GetParameter(1); //MPV
343  FitResults[1] = MyLandau->GetParError(1); //MPV error
344  FitResults[2] = MyLandau->GetParameter(2); //Width
345  FitResults[3] = MyLandau->GetParError(2); //Width error
346  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
347 
348  delete MyLandau;
349 }
350 
351 bool SiStripGainFromCalibTree::IsGoodLandauFit(double* FitResults){
352  if(FitResults[0] < 2 )return false;
353  if(FitResults[1] > MaxMPVError )return false;
354  if(FitResults[4] > MaxChi2OverNDF)return false;
355  return true;
356 }
357 
358 
359 void SiStripGainFromCalibTree::algoAnalyzeTheTree()
360 {
361  for(unsigned int i=0;i<VInputFiles.size();i++){
362  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
363  TChain* tree = new TChain("gainCalibrationTree/tree");
364  tree->Add(VInputFiles[i].c_str());
365 
366  TString EventPrefix("");
367  TString EventSuffix("");
368 
369  TString TrackPrefix("track");
370  TString TrackSuffix("");
371 
372  TString CalibPrefix("GainCalibration");
373  TString CalibSuffix("");
374 
375  unsigned int eventnumber = 0; tree->SetBranchAddress(EventPrefix + "event" + EventSuffix, &eventnumber , NULL);
376  unsigned int runnumber = 0; tree->SetBranchAddress(EventPrefix + "run" + EventSuffix, &runnumber , NULL);
377  std::vector<bool>* TrigTech = 0; tree->SetBranchAddress(EventPrefix + "TrigTech" + EventSuffix, &TrigTech , NULL);
378 
379  std::vector<double>* trackchi2ndof = 0; tree->SetBranchAddress(TrackPrefix + "chi2ndof" + TrackSuffix, &trackchi2ndof , NULL);
380  std::vector<float>* trackp = 0; tree->SetBranchAddress(TrackPrefix + "momentum" + TrackSuffix, &trackp , NULL);
381  std::vector<float>* trackpt = 0; tree->SetBranchAddress(TrackPrefix + "pt" + TrackSuffix, &trackpt , NULL);
382  std::vector<double>* tracketa = 0; tree->SetBranchAddress(TrackPrefix + "eta" + TrackSuffix, &tracketa , NULL);
383  std::vector<double>* trackphi = 0; tree->SetBranchAddress(TrackPrefix + "phi" + TrackSuffix, &trackphi , NULL);
384  std::vector<unsigned int>* trackhitsvalid = 0; tree->SetBranchAddress(TrackPrefix + "hitsvalid" + TrackSuffix, &trackhitsvalid, NULL);
385 
386  std::vector<int>* trackindex = 0; tree->SetBranchAddress(CalibPrefix + "trackindex" + CalibSuffix, &trackindex , NULL);
387  std::vector<unsigned int>* rawid = 0; tree->SetBranchAddress(CalibPrefix + "rawid" + CalibSuffix, &rawid , NULL);
388  std::vector<float>* localdirx = 0; tree->SetBranchAddress(CalibPrefix + "localdirx" + CalibSuffix, &localdirx , NULL);
389  std::vector<float>* localdiry = 0; tree->SetBranchAddress(CalibPrefix + "localdiry" + CalibSuffix, &localdiry , NULL);
390  std::vector<float>* localdirz = 0; tree->SetBranchAddress(CalibPrefix + "localdirz" + CalibSuffix, &localdirz , NULL);
391  std::vector<unsigned short>* firststrip = 0; tree->SetBranchAddress(CalibPrefix + "firststrip" + CalibSuffix, &firststrip , NULL);
392  std::vector<unsigned short>* nstrips = 0; tree->SetBranchAddress(CalibPrefix + "nstrips" + CalibSuffix, &nstrips , NULL);
393  std::vector<bool>* saturation = 0; tree->SetBranchAddress(CalibPrefix + "saturation" + CalibSuffix, &saturation , NULL);
394  std::vector<bool>* overlapping = 0; tree->SetBranchAddress(CalibPrefix + "overlapping" + CalibSuffix, &overlapping , NULL);
395  std::vector<bool>* farfromedge = 0; tree->SetBranchAddress(CalibPrefix + "farfromedge" + CalibSuffix, &farfromedge , NULL);
396  std::vector<unsigned int>* charge = 0; tree->SetBranchAddress(CalibPrefix + "charge" + CalibSuffix, &charge , NULL);
397  std::vector<float>* path = 0; tree->SetBranchAddress(CalibPrefix + "path" + CalibSuffix, &path , NULL);
398  std::vector<float>* chargeoverpath = 0; tree->SetBranchAddress(CalibPrefix + "chargeoverpath" + CalibSuffix, &chargeoverpath, NULL);
399  std::vector<unsigned char>* amplitude = 0; tree->SetBranchAddress(CalibPrefix + "amplitude" + CalibSuffix, &amplitude , NULL);
400  std::vector<double>* gainused = 0; tree->SetBranchAddress(CalibPrefix + "gainused" + CalibSuffix, &gainused , NULL);
401 
402 
403  printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));
404  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
405  printf("Looping on the Tree :");
406  int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
407  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
408 // for (unsigned int ientry = 0; ientry < tree->GetEntries() && ientry<50000; ientry++) {
409  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
410  tree->GetEntry(ientry);
411 
412  if(runnumber<SRun)SRun=runnumber;
413  if(runnumber>ERun)ERun=runnumber;
414 
415  NEvent++;
416  NTrack+=(*trackp).size();
417 
418  unsigned int FirstAmplitude=0;
419  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
420  FirstAmplitude+=(*nstrips)[i];
421  int TI = (*trackindex)[i];
422  if((*tracketa )[TI] < MinTrackEta )continue;
423  if((*tracketa )[TI] > MaxTrackEta )continue;
424  if((*trackp )[TI] < MinTrackMomentum )continue;
425  if((*trackp )[TI] > MaxTrackMomentum )continue;
426  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
427  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
428  if((*farfromedge)[i] == false )continue;
429  if((*overlapping)[i] == true )continue;
430  if((*saturation )[i] && !AllowSaturation)continue;
431  if((*nstrips )[i] > MaxNrStrips )continue;
432  NCluster++;
433 
434  stAPVGain* APV = APVsColl[((*rawid)[i]<<3) | ((*firststrip)[i]/128)];
435 
436 // printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
437 
438 
439 
440  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
441  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
442  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
443 
444  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
445  //printf("Amplitudes: ");
446  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
447  //printf("\n");
448 
449 
450  int Charge = 0;
451  if(useCalibration || !FirstSetOfConstants){
452  bool Saturation = false;
453  for(unsigned int s=0;s<(*nstrips)[i];s++){
454  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
455  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
456  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
457  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
458  if(StripCharge>1024){
459  StripCharge = 255;
460  Saturation = true;
461  }else if(StripCharge>254){
462  StripCharge = 254;
463  Saturation = true;
464  }
465  Charge += StripCharge;
466  }
467  if(Saturation && !AllowSaturation)continue;
468  }else{
469  Charge = (*charge)[i];
470  }
471 
472  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
473 
474  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
475  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
476  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
477  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
478  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
479 
480 
481  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill((*path)[i],Charge);
482  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill((*path)[i],Charge);
483  }else if(APV->SubDet==StripSubdetector::TID){
484  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
485  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
486  }
487  }else if(APV->SubDet==StripSubdetector::TEC){
488  if(APV->Eta<0){
489  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
490  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
491  }
492  }else if(APV->Eta>0){
493  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
494  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
495  }
496  }
497  }
498 
499  }// END OF ON-CLUSTER LOOP
500  }printf("\n");// END OF EVENT LOOP
501 
502  }
503 }
504 
505 
506 
507 void SiStripGainFromCalibTree::algoComputeMPVandGain() {
508  unsigned int I=0;
509  TH1D* Proj = NULL;
510  double FitResults[5];
511  double MPVmean = 300;
512 
513  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
514  printf("Fitting Charge Distribution :");
515  int TreeStep = APVsColl.size()/50;
516  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
517  if(I%TreeStep==0){printf(".");fflush(stdout);}
518  //if(I>1000)break;
519  stAPVGain* APV = it->second;
520 
521  Proj = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV->Bin,APV->Bin,"e"));
522  if(!Proj)continue;
523 
524  if(CalibrationLevel==0){
525  }else if(CalibrationLevel==1){
526  int SecondAPVId = APV->APVId;
527  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
528  stAPVGain* APV2 = APVsColl[(APV->DetId<<3) | SecondAPVId];
529  TH1D* Proj2 = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
530  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
531  }else if(CalibrationLevel==2){
532  for(unsigned int i=0;i<6;i++){
533  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
534  tmpit = APVsColl.find((APV->DetId<<3) | i);
535  if(tmpit==APVsColl.end())continue;
536  stAPVGain* APV2 = tmpit->second;
537  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
538  TH1D* Proj2 = (TH1D*)(Charge_Vs_Index->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
539 // if(Proj2 && APV->DetId==369171124)printf("B) DetId %6i APVId %1i --> NEntries = %f\n",APV2->DetId, APV2->APVId, Proj2->GetEntries());
540  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
541  }
542  }else{
543  CalibrationLevel = 0;
544  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
545  }
546 
547  getPeakOfLandau(Proj,FitResults);
548  APV->FitMPV = FitResults[0];
549  APV->FitMPVErr = FitResults[1];
550  APV->FitWidth = FitResults[2];
551  APV->FitWidthErr = FitResults[3];
552  APV->FitChi2 = FitResults[4];
553  APV->NEntries = Proj->GetEntries();
554 
555  if(APV->FitMPV>0){
556  APV->Gain = APV->FitMPV / MPVmean;
557  GOOD++;
558  }else{
559  APV->Gain = 1;
560  BAD++;
561  }
562  if(APV->Gain<=0) APV->Gain = 1;
563 // if(!FirstSetOfConstants) APV->Gain *= APV->PreviousGain;
564 
565 
566  //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);
567  delete Proj;
568  }printf("\n");
569 
570  storeOnTree();
571 }
572 
573 
574 void SiStripGainFromCalibTree::storeOnTree()
575 {
576  unsigned int tree_Index;
577  unsigned int tree_Bin;
578  unsigned int tree_DetId;
579  unsigned char tree_APVId;
580  unsigned char tree_SubDet;
581  float tree_x;
582  float tree_y;
583  float tree_z;
584  float tree_Eta;
585  float tree_R;
586  float tree_Phi;
587  float tree_Thickness;
588  float tree_FitMPV;
589  float tree_FitMPVErr;
590  float tree_FitWidth;
591  float tree_FitWidthErr;
592  float tree_FitChi2NDF;
593  double tree_Gain;
594  double tree_PrevGain;
595  double tree_NEntries;
596  bool tree_isMasked;
597 
598  TTree* MyTree;
599  MyTree = tfs->make<TTree> ("APVGain","APVGain");
600  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
601  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
602  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
603  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
604  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
605  MyTree->Branch("x" ,&tree_x ,"x/F");
606  MyTree->Branch("y" ,&tree_y ,"y/F");
607  MyTree->Branch("z" ,&tree_z ,"z/F");
608  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
609  MyTree->Branch("R" ,&tree_R ,"R/F");
610  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
611  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
612  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
613  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
614  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
615  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
616  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
617  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
618  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
619  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
620  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
621 
622 
623  FILE* Gains = fopen(OutputGains.c_str(),"w");
624  fprintf(Gains,"NEvents = %i\n",NEvent);
625  fprintf(Gains,"NTracks = %i\n",NTrack);
626  fprintf(Gains,"NClusters = %i\n",NCluster);
627  fprintf(Gains,"Number of APVs = %lu\n",static_cast<unsigned long>(APVsColl.size()));
628  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD));
629 
630  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
631  stAPVGain* APV = APVsCollOrdered[a];
632  if(APV==NULL)continue;
633 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain);
634  fprintf(Gains,"%i | %i | PreviousGain = %7.5f NewGain = %7.5f\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain);
635 
636  tree_Index = APV->Index;
637  tree_Bin = APV->Bin;
638  tree_DetId = APV->DetId;
639  tree_APVId = APV->APVId;
640  tree_SubDet = APV->SubDet;
641  tree_x = APV->x;
642  tree_y = APV->y;
643  tree_z = APV->z;
644  tree_Eta = APV->Eta;
645  tree_R = APV->R;
646  tree_Phi = APV->Phi;
647  tree_Thickness = APV->Thickness;
648  tree_FitMPV = APV->FitMPV;
649  tree_FitMPVErr = APV->FitMPVErr;
650  tree_FitWidth = APV->FitWidth;
651  tree_FitWidthErr= APV->FitWidthErr;
652  tree_FitChi2NDF = APV->FitChi2;
653  tree_Gain = APV->Gain;
654  tree_PrevGain = APV->PreviousGain;
655  tree_NEntries = APV->NEntries;
656  tree_isMasked = APV->isMasked;
657 
658  MyTree->Fill();
659  }
660  fclose(Gains);
661 
662 
663 }
664 
665 
666 
667 SiStripApvGain* SiStripGainFromCalibTree::getNewObject()
668 {
669  SiStripApvGain* obj = new SiStripApvGain();
670  std::vector<float>* theSiStripVector = NULL;
671  unsigned int PreviousDetId = 0;
672  for(unsigned int a=0;a<APVsCollOrdered.size();a++)
673  {
674  stAPVGain* APV = APVsCollOrdered[a];
675  if(APV==NULL){ printf("Bug\n"); continue; }
676  if(APV->DetId != PreviousDetId){
677  if(theSiStripVector!=NULL){
678  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
679  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
680  }
681  theSiStripVector = new std::vector<float>;
682  PreviousDetId = APV->DetId;
683  }
684  theSiStripVector->push_back(APV->Gain);
685  }
686  if(theSiStripVector!=NULL){
687  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
688  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
689  }
690 
691  return obj;
692 }
693 
694 
695 SiStripGainFromCalibTree::~SiStripGainFromCalibTree()
696 {
697 }
698 
699 void SiStripGainFromCalibTree::MakeCalibrationMap(){
700  if(!useCalibration)return;
701 
702 
703  TChain* t1 = new TChain("SiStripCalib/APVGain");
704  t1->Add(m_calibrationPath.c_str());
705 
706  unsigned int tree_DetId;
707  unsigned char tree_APVId;
708  double tree_Gain;
709 
710  t1->SetBranchAddress("DetId" ,&tree_DetId );
711  t1->SetBranchAddress("APVId" ,&tree_APVId );
712  t1->SetBranchAddress("Gain" ,&tree_Gain );
713 
714  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
715  t1->GetEntry(ientry);
716  stAPVGain* APV = APVsColl[(tree_DetId<<3) | (unsigned int)tree_APVId];
717  APV->CalibGain = tree_Gain;
718  }
719 
720 }
721 
722 
723 
724 
725 DEFINE_FWK_MODULE(SiStripGainFromCalibTree);
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
Basic3DVector::y
T y() const
Cartesian y coordinate.
Definition: extBasic3DVector.h:102
SiStripGainFromCalibTree::tfs
edm::Service< TFileService > tfs
Definition: SiStripGainFromCalibTree.cc:133
Basic3DVector::x
T x() const
Cartesian x coordinate.
Definition: extBasic3DVector.h:99
SiStripGainFromCalibTree::getNewObject
SiStripApvGain * getNewObject() override
Definition: SiStripGainFromCalibTree.cc:667
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
SiStripGainFromCalibTree::SiStripGainFromCalibTree
SiStripGainFromCalibTree(const edm::ParameterSet &)
Definition: SiStripGainFromCalibTree.cc:186
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
Basic3DVector::phi
Geom::Phi< T > phi() const
Definition: extBasic3DVector.h:137
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
StripGeomDetUnit::specificTopology
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
Definition: StripGeomDetUnit.cc:19
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
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
Bounds::thickness
virtual float thickness() const =0
SiStripGainFromCalibTree::APVsCollOrdered
std::vector< stAPVGain * > APVsCollOrdered
Definition: SiStripGainFromCalibTree.cc:182
SiStripGainFromCalibTree::APVsColl
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
Definition: SiStripGainFromCalibTree.cc:183
SiStripGainFromCalibTree::algoBeginJob
virtual void algoBeginJob(const edm::EventSetup &) override
Definition: SiStripGainFromCalibTree.cc:219
iEvent
int iEvent
Definition: GenABIO.cc:243
GeomDet::position
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:41
Basic3DVector::z
T z() const
Cartesian z coordinate.
Definition: extBasic3DVector.h:105
SiStripGainFromCalibTree::IsGoodLandauFit
bool IsGoodLandauFit(double *FitResults)
Definition: SiStripGainFromCalibTree.cc:351
SiStripApvGain::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripApvGain.h:25
j
int j
Definition: DBlmapReader.cc:9
Exhume::I
const std::complex< double > I
Definition: I.h:8
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
SiStripGainFromCalibTree::algoEndJob
virtual void algoEndJob() override
Definition: SiStripGainFromCalibTree.cc:322
DetId
Definition: DetId.h:18
SiStripGainFromCalibTree::algoAnalyze
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
Definition: SiStripGainFromCalibTree.cc:317
SiStripGainFromCalibTree::getPeakOfLandau
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
Definition: SiStripGainFromCalibTree.cc:326
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
a
double a
Definition: hdecay.h:121
Basic3DVector::transverse
T transverse() const
Another name for perp()
Definition: extBasic3DVector.h:130
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:37
T
long double T
Definition: Basic3DVectorLD.h:59
PV3DBase::basicVector
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:56
Generated for CMSSW Reference Manual by doxygen 1.8.5