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