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 
15 
21 
28 
30 
33 
46 
49 
52 
57 
60 
63 
64 
65 #include "TFile.h"
66 #include "TObjString.h"
67 #include "TString.h"
68 #include "TH1F.h"
69 #include "TH2F.h"
70 #include "TProfile.h"
71 #include "TF1.h"
72 #include "TROOT.h"
73 #include "TTree.h"
74 #include "TChain.h"
75 
76 #include <ext/hash_map>
77 
78 
79 
80 using namespace edm;
81 using namespace reco;
82 using namespace std;
83 using namespace __gnu_cxx;
84 using __gnu_cxx::hash_map;
85 using __gnu_cxx::hash;
86 
87 struct stAPVGain{
88  unsigned int Index;
89  int Bin;
90  unsigned int DetId;
91  unsigned int APVId;
92  unsigned int SubDet;
93  float x;
94  float y;
95  float z;
96  float Eta;
97  float R;
98  float Phi;
99  float Thickness;
100  double FitMPV;
101  double FitMPVErr;
102  double FitWidth;
103  double FitWidthErr;
104  double FitChi2;
105  double FitNorm;
106  double Gain;
107  double CalibGain;
108  double PreviousGain;
110  double NEntries;
111  TH1F* HCharge;
112  TH1F* HChargeN;
113  bool isMasked;
114 };
115 
116 class SiStripGainFromCalibTree : public ConditionDBWriter<SiStripApvGain> {
117  public:
120 
121 
122  private:
123 
124 
125  virtual void algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) override;
126  virtual void algoEndRun (const edm::Run& run, const edm::EventSetup& iSetup) override;
127  virtual void algoBeginJob (const edm::EventSetup& iSetup) override;
128  virtual void algoEndJob () override;
129  virtual void algoAnalyze (const edm::Event &, const edm::EventSetup &) override;
130 
131  void merge(TH2* A, TH2* B); //needed to add histograms with different number of bins
132  void algoAnalyzeTheTree();
133  void algoComputeMPVandGain();
134 
135  bool IsFarFromBorder(TrajectoryStateOnSurface* trajState, const uint32_t detid, const edm::EventSetup* iSetup);
136  void getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange=50, double HighRange=5400);
137  bool IsGoodLandauFit(double* FitResults);
138  void storeOnTree(TFileService* tfs);
139  void MakeCalibrationMap();
140  bool produceTagFilter();
141 
142 
143 
144  SiStripApvGain* getNewObject() override;
145 
149  double MinNrEntries;
150  double MaxMPVError;
154  double MinTrackEta;
155  double MaxTrackEta;
156  unsigned int MaxNrStrips;
157  unsigned int MinTrackHits;
164 
168 
171 
175  vector<string> VInputFiles;
176 
187 
188  unsigned int NEvent;
189  unsigned int NTrack;
190  unsigned int NClusterStrip;
191  unsigned int NClusterPixel;
194  unsigned int SRun;
195  unsigned int ERun;
196  unsigned int GOOD;
197  unsigned int BAD;
198  unsigned int MASKED;
199 
202 
203  private :
204  class isEqual{
205  public:
206  template <class T> bool operator () (const T& PseudoDetId1, const T& PseudoDetId2) { return PseudoDetId1==PseudoDetId2; }
207  };
208 
209  std::vector<stAPVGain*> APVsCollOrdered;
210  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual > APVsColl;
211 };
212 
214  if(A->GetNbinsX() == B->GetNbinsX()){
215  A->Add(B);
216  }else{
217  for(int x=0;x<=B->GetNbinsX()+1; x++){
218  for(int y=0;y<=B->GetNbinsY()+1; y++){
219  A->SetBinContent(x,y,A->GetBinContent(x,y)+B->GetBinContent(x,y));
220  }}
221  }
222 }
223 
224 
225 
227 {
228  Charge_Vs_Index = NULL;//make sure this is initialized to NULL
229  OutputGains = iConfig.getParameter<std::string>("OutputGains");
231 
232  AlgoMode = iConfig.getUntrackedParameter<std::string>("AlgoMode", "CalibTree");
233  MinNrEntries = iConfig.getUntrackedParameter<double> ("minNrEntries" , 20);
234  MaxMPVError = iConfig.getUntrackedParameter<double> ("maxMPVError" , 500.0);
235  MaxChi2OverNDF = iConfig.getUntrackedParameter<double> ("maxChi2OverNDF" , 5.0);
236  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 3.0);
237  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
238  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
239  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
240  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 2);
241  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 8);
242  MaxTrackChiOverNdf = iConfig.getUntrackedParameter<double> ("MaxTrackChiOverNdf" , 3);
243  AllowSaturation = iConfig.getUntrackedParameter<bool> ("AllowSaturation" , false);
244  FirstSetOfConstants = iConfig.getUntrackedParameter<bool> ("FirstSetOfConstants", true);
245  Validation = iConfig.getUntrackedParameter<bool> ("Validation" , false);
246  OldGainRemoving = iConfig.getUntrackedParameter<bool> ("OldGainRemoving" , false);
247 
248  CalibrationLevel = iConfig.getUntrackedParameter<int> ("CalibrationLevel" , 0);
249  VInputFiles = iConfig.getUntrackedParameter<vector<string> > ("InputFiles");
250 
251 
252  useCalibration = iConfig.getUntrackedParameter<bool> ("UseCalibration" , false);
253  m_calibrationPath = iConfig.getUntrackedParameter<string> ("calibrationPath");
254  harvestingMode = iConfig.getUntrackedParameter<bool> ("harvestingMode" , false);
255 
256  tagCondition_NClusters = iConfig.getUntrackedParameter<double> ("NClustersForTagProd" , 2E8);
257  tagCondition_GoodFrac = iConfig.getUntrackedParameter<double> ("GoodFracForTagProd" , 0.95);
258 
259  saveSummary = iConfig.getUntrackedParameter<bool> ("saveSummary" , false);
260 
262  dbe->setVerbose(10);
263 
264  tracksToken = consumes<edm::View<reco::Track>>(theTracksLabel);
265  associationsToken = consumes<TrajTrackAssociationCollection>(theTracksLabel);
266 }
267 
268 
270 {
271  dbe->setCurrentFolder("AlCaReco/SiStripGains/");
272  if(AlgoMode != "PCL" or harvestingMode)dbe->setCurrentFolder("AlCaReco/SiStripGainsHarvesting/");
273 
274  Charge_Vs_Index = dbe->book2D("Charge_Vs_Index" , "Charge_Vs_Index" , 88625, 0 , 88624,2000,0,4000);
275  Charge_Vs_Index_Absolute = dbe->book2D("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 88625, 0 , 88624,1000,0,4000);
276 // Charge_Vs_Index = dbe->book2D("Charge_Vs_Index" , "Charge_Vs_Index" , 72785, 0 , 72784,1000,0,2000);
277 // Charge_Vs_Index_Absolute = dbe->book2D("Charge_Vs_Index_Absolute" , "Charge_Vs_Index_Absolute" , 72785, 0 , 72784, 500,0,2000);
278  Charge_Vs_PathlengthTIB = dbe->book2D("Charge_Vs_PathlengthTIB" , "Charge_Vs_PathlengthTIB" , 20 , 0.3 , 1.3 , 250,0,2000);
279  Charge_Vs_PathlengthTOB = dbe->book2D("Charge_Vs_PathlengthTOB" , "Charge_Vs_PathlengthTOB" , 20 , 0.3 , 1.3 , 250,0,2000);
280  Charge_Vs_PathlengthTIDP = dbe->book2D("Charge_Vs_PathlengthTIDP" , "Charge_Vs_PathlengthTIDP" , 20 , 0.3 , 1.3 , 250,0,2000);
281  Charge_Vs_PathlengthTIDM = dbe->book2D("Charge_Vs_PathlengthTIDM" , "Charge_Vs_PathlengthTIDM" , 20 , 0.3 , 1.3 , 250,0,2000);
282  Charge_Vs_PathlengthTECP1 = dbe->book2D("Charge_Vs_PathlengthTECP1", "Charge_Vs_PathlengthTECP1", 20 , 0.3 , 1.3 , 250,0,2000);
283  Charge_Vs_PathlengthTECP2 = dbe->book2D("Charge_Vs_PathlengthTECP2", "Charge_Vs_PathlengthTECP2", 20 , 0.3 , 1.3 , 250,0,2000);
284  Charge_Vs_PathlengthTECM1 = dbe->book2D("Charge_Vs_PathlengthTECM1", "Charge_Vs_PathlengthTECM1", 20 , 0.3 , 1.3 , 250,0,2000);
285  Charge_Vs_PathlengthTECM2 = dbe->book2D("Charge_Vs_PathlengthTECM2", "Charge_Vs_PathlengthTECM2", 20 , 0.3 , 1.3 , 250,0,2000);
286 
288  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
289  auto const & Det = tkGeom->dets();
290 
291  NPixelDets = 0;
292  NStripAPVs = 0;
293  unsigned int Index=0;
294  for(unsigned int i=0;i<Det.size();i++){
295  DetId Detid = Det[i]->geographicalId();
296  int SubDet = Detid.subdetId();
297 
298  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
299  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
300 
301  auto DetUnit = dynamic_cast<const StripGeomDetUnit*> (Det[i]);
302  if(!DetUnit)continue;
303 
304  const StripTopology& Topo = DetUnit->specificTopology();
305  unsigned int NAPV = Topo.nstrips()/128;
306 
307  for(unsigned int j=0;j<NAPV;j++){
308  stAPVGain* APV = new stAPVGain;
309  APV->Index = Index;
310  APV->Bin = -1;
311  APV->DetId = Detid.rawId();
312  APV->APVId = j;
313  APV->SubDet = SubDet;
314  APV->FitMPV = -1;
315  APV->FitMPVErr = -1;
316  APV->FitWidth = -1;
317  APV->FitWidthErr = -1;
318  APV->FitChi2 = -1;
319  APV->FitNorm = -1;
320  APV->Gain = -1;
321  APV->PreviousGain = 1;
322  APV->PreviousGainTick = 1;
323  APV->x = DetUnit->position().basicVector().x();
324  APV->y = DetUnit->position().basicVector().y();
325  APV->z = DetUnit->position().basicVector().z();
326  APV->Eta = DetUnit->position().basicVector().eta();
327  APV->Phi = DetUnit->position().basicVector().phi();
328  APV->R = DetUnit->position().basicVector().transverse();
329  APV->Thickness = DetUnit->surface().bounds().thickness();
330  APV->NEntries = 0;
331  APV->isMasked = false;
332 
333  APVsCollOrdered.push_back(APV);
334  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
335  Index++;
336  NStripAPVs++;
337  }
338  }
339  }
340 
341  for(unsigned int i=0;i<Det.size();i++){ //Make two loop such that the Pixel information is added at the end --> make transition simpler
342  DetId Detid = Det[i]->geographicalId();
343  int SubDet = Detid.subdetId();
345  auto DetUnit = dynamic_cast<const PixelGeomDetUnit*> (Det[i]);
346  if(!DetUnit) continue;
347 
348  const PixelTopology& Topo = DetUnit->specificTopology();
349  unsigned int NROCRow = Topo.nrows()/(80.);
350  unsigned int NROCCol = Topo.ncolumns()/(52.);
351 
352  for(unsigned int j=0;j<NROCRow;j++){
353  for(unsigned int i=0;i<NROCCol;i++){
354 
355  stAPVGain* APV = new stAPVGain;
356  APV->Index = Index;
357  APV->Bin = -1;
358  APV->DetId = Detid.rawId();
359  APV->APVId = (j<<3 | i);
360  APV->SubDet = SubDet;
361  APV->FitMPV = -1;
362  APV->FitMPVErr = -1;
363  APV->FitWidth = -1;
364  APV->FitWidthErr = -1;
365  APV->FitChi2 = -1;
366  APV->Gain = -1;
367  APV->PreviousGain = 1;
368  APV->x = DetUnit->position().basicVector().x();
369  APV->y = DetUnit->position().basicVector().y();
370  APV->z = DetUnit->position().basicVector().z();
371  APV->Eta = DetUnit->position().basicVector().eta();
372  APV->Phi = DetUnit->position().basicVector().phi();
373  APV->R = DetUnit->position().basicVector().transverse();
374  APV->Thickness = DetUnit->surface().bounds().thickness();
375  APV->isMasked = false; //SiPixelQuality_->IsModuleBad(Detid.rawId());
376  APV->NEntries = 0;
377 
378  APVsCollOrdered.push_back(APV);
379  APVsColl[(APV->DetId<<4) | APV->APVId] = APV;
380  Index++;
381  NPixelDets++;
382  }}
383  }
384  }
385 
386 
388 
389  NEvent = 0;
390  NTrack = 0;
391  NClusterStrip = 0;
392  NClusterPixel = 0;
393  SRun = 1<<31;
394  ERun = 0;
395  GOOD = 0;
396  BAD = 0;
397  MASKED = 0;
398 }
399 
400 
401 
403 {
404  edm::ESHandle<SiStripGain> gainHandle;
405  iSetup.get<SiStripGainRcd>().get(gainHandle);
406  if(!gainHandle.isValid()){edm::LogError("SiStripGainFromCalibTree")<< "gainHandle is not valid\n"; exit(0);}
407 
408  edm::ESHandle<SiStripQuality> SiStripQuality_;
409  iSetup.get<SiStripQualityRcd>().get(SiStripQuality_);
410 
411  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
413 
414  APV->isMasked = SiStripQuality_->IsApvBad(APV->DetId,APV->APVId);
415 // if(!FirstSetOfConstants){
416  if(gainHandle->getNumberOfTags()!=2){edm::LogError("SiStripGainFromCalibTree")<< "NUMBER OF GAIN TAG IS EXPECTED TO BE 2\n";fflush(stdout);exit(0);};
417  float newPreviousGain = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 1),1);
418  if(APV->PreviousGain!=1 and newPreviousGain!=APV->PreviousGain)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: ParticleGain in the global tag changed\n";
419  APV->PreviousGain = newPreviousGain;
420 
421  float newPreviousGainTick = gainHandle->getApvGain(APV->APVId,gainHandle->getRange(APV->DetId, 0),1);
422  if(APV->PreviousGainTick!=1 and newPreviousGainTick!=APV->PreviousGainTick)edm::LogWarning("SiStripGainFromCalibTree")<< "WARNING: TickMarkGain in the global tag changed\n";
423  APV->PreviousGainTick = newPreviousGainTick;
424 
425 
426  //printf("DETID = %7i APVID=%1i Previous Gain=%8.4f (G1) x %8.4f (G2)\n",APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain);
427 
428 
429 
430 // }
431  }
432 }
433 
435  if(AlgoMode == "PCL" && !harvestingMode)return;//nothing to do in that case
436 
437  if(AlgoMode == "PCL" and harvestingMode){
438  // Load the 2D histograms from the DQM objects
439  // When running in AlCaHarvesting mode the histos are already booked and should be just retrieved from
440  // DQMStore so that they can be used in the fit
441 
442 // edm::LogInfo("SiStripGainFromCalibTree")<< "Harvesting " << (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F()->GetEntries() << " more clusters\n";
443 
444  merge(Charge_Vs_Index ->getTH2F(), (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F() );
445 // Charge_Vs_Index ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index"))->getTH2F());
446  merge(Charge_Vs_Index_Absolute ->getTH2F(), (dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index_Absolute"))->getTH2F() );
447 // Charge_Vs_Index_Absolute ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_Index_Absolute"))->getTH2F());
448  Charge_Vs_PathlengthTIB ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIB"))->getTH2F());
449  Charge_Vs_PathlengthTOB ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTOB"))->getTH2F());
450  Charge_Vs_PathlengthTIDP ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDP"))->getTH2F());
451  Charge_Vs_PathlengthTIDM ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTIDM"))->getTH2F());
452  Charge_Vs_PathlengthTECP1 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP1"))->getTH2F());
453  Charge_Vs_PathlengthTECP2 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECP2"))->getTH2F());
454  Charge_Vs_PathlengthTECM1 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM1"))->getTH2F());
455  Charge_Vs_PathlengthTECM2 ->getTH2F()->Add((dbe->get("AlCaReco/SiStripGains/Charge_Vs_PathlengthTECM2"))->getTH2F());
456  }
457 }
458 void
460  if(AlgoMode == "PCL" && !harvestingMode)return;//nothing to do in that case
461 
462  if(AlgoMode == "CalibTree"){
463  // Loop on calibTrees to fill the 2D histograms
465  }else if(harvestingMode){
466  NClusterStrip = Charge_Vs_Index->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 );
467  NClusterPixel = Charge_Vs_Index->getTH2F()->Integral(NStripAPVs+2, NStripAPVs+NPixelDets+2, 0, 99999 );
468  }
469 
470  // Now that we have the full statistics we can extract the information of the 2D histograms
472 
473  if(AlgoMode != "PCL" or saveSummary){
474  //also save the 2D monitor elements to this file as TH2D tfs
476  tfs->make<TH2F> (*Charge_Vs_Index->getTH2F());
486 
487 
488  storeOnTree(tfs);
489  }
490 }
491 
492 
493 void SiStripGainFromCalibTree::getPeakOfLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
494 {
495  FitResults[0] = -0.5; //MPV
496  FitResults[1] = 0; //MPV error
497  FitResults[2] = -0.5; //Width
498  FitResults[3] = 0; //Width error
499  FitResults[4] = -0.5; //Fit Chi2/NDF
500  FitResults[5] = 0; //Normalization
501 
502  if( InputHisto->GetEntries() < MinNrEntries)return;
503 
504  // perform fit with standard landau
505  TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
506  MyLandau->SetParameter(1,300);
507  InputHisto->Fit(MyLandau,"0QR WW");
508 
509  // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
510  FitResults[0] = MyLandau->GetParameter(1); //MPV
511  FitResults[1] = MyLandau->GetParError(1); //MPV error
512  FitResults[2] = MyLandau->GetParameter(2); //Width
513  FitResults[3] = MyLandau->GetParError(2); //Width error
514  FitResults[4] = MyLandau->GetChisquare() / MyLandau->GetNDF(); //Fit Chi2/NDF
515  FitResults[5] = MyLandau->GetParameter(0);
516 
517  delete MyLandau;
518 }
519 
521  if(FitResults[0] <= 0 )return false;
522 // if(FitResults[1] > MaxMPVError )return false;
523 // if(FitResults[4] > MaxChi2OverNDF)return false;
524  return true;
525 }
526 
527 
529 {
530  for(unsigned int i=0;i<VInputFiles.size();i++){
531  printf("Openning file %3i/%3i --> %s\n",i+1, (int)VInputFiles.size(), (char*)(VInputFiles[i].c_str())); fflush(stdout);
532  TChain* tree = new TChain("gainCalibrationTree/tree");
533  tree->Add(VInputFiles[i].c_str());
534 
535  TString EventPrefix("");
536  TString EventSuffix("");
537 
538  TString TrackPrefix("track");
539  TString TrackSuffix("");
540 
541  TString CalibPrefix("GainCalibration");
542  TString CalibSuffix("");
543 
544  unsigned int eventnumber = 0; tree->SetBranchAddress(EventPrefix + "event" + EventSuffix, &eventnumber , NULL);
545  unsigned int runnumber = 0; tree->SetBranchAddress(EventPrefix + "run" + EventSuffix, &runnumber , NULL);
546  std::vector<bool>* TrigTech = 0; tree->SetBranchAddress(EventPrefix + "TrigTech" + EventSuffix, &TrigTech , NULL);
547 
548  std::vector<double>* trackchi2ndof = 0; tree->SetBranchAddress(TrackPrefix + "chi2ndof" + TrackSuffix, &trackchi2ndof , NULL);
549  std::vector<float>* trackp = 0; tree->SetBranchAddress(TrackPrefix + "momentum" + TrackSuffix, &trackp , NULL);
550  std::vector<float>* trackpt = 0; tree->SetBranchAddress(TrackPrefix + "pt" + TrackSuffix, &trackpt , NULL);
551  std::vector<double>* tracketa = 0; tree->SetBranchAddress(TrackPrefix + "eta" + TrackSuffix, &tracketa , NULL);
552  std::vector<double>* trackphi = 0; tree->SetBranchAddress(TrackPrefix + "phi" + TrackSuffix, &trackphi , NULL);
553  std::vector<unsigned int>* trackhitsvalid = 0; tree->SetBranchAddress(TrackPrefix + "hitsvalid" + TrackSuffix, &trackhitsvalid, NULL);
554 
555  std::vector<int>* trackindex = 0; tree->SetBranchAddress(CalibPrefix + "trackindex" + CalibSuffix, &trackindex , NULL);
556  std::vector<unsigned int>* rawid = 0; tree->SetBranchAddress(CalibPrefix + "rawid" + CalibSuffix, &rawid , NULL);
557  std::vector<float>* localdirx = 0; tree->SetBranchAddress(CalibPrefix + "localdirx" + CalibSuffix, &localdirx , NULL);
558  std::vector<float>* localdiry = 0; tree->SetBranchAddress(CalibPrefix + "localdiry" + CalibSuffix, &localdiry , NULL);
559  std::vector<float>* localdirz = 0; tree->SetBranchAddress(CalibPrefix + "localdirz" + CalibSuffix, &localdirz , NULL);
560  std::vector<unsigned short>* firststrip = 0; tree->SetBranchAddress(CalibPrefix + "firststrip" + CalibSuffix, &firststrip , NULL);
561  std::vector<unsigned short>* nstrips = 0; tree->SetBranchAddress(CalibPrefix + "nstrips" + CalibSuffix, &nstrips , NULL);
562  std::vector<bool>* saturation = 0; tree->SetBranchAddress(CalibPrefix + "saturation" + CalibSuffix, &saturation , NULL);
563  std::vector<bool>* overlapping = 0; tree->SetBranchAddress(CalibPrefix + "overlapping" + CalibSuffix, &overlapping , NULL);
564  std::vector<bool>* farfromedge = 0; tree->SetBranchAddress(CalibPrefix + "farfromedge" + CalibSuffix, &farfromedge , NULL);
565  std::vector<unsigned int>* charge = 0; tree->SetBranchAddress(CalibPrefix + "charge" + CalibSuffix, &charge , NULL);
566  std::vector<float>* path = 0; tree->SetBranchAddress(CalibPrefix + "path" + CalibSuffix, &path , NULL);
567  std::vector<float>* chargeoverpath = 0; tree->SetBranchAddress(CalibPrefix + "chargeoverpath" + CalibSuffix, &chargeoverpath, NULL);
568  std::vector<unsigned char>* amplitude = 0; tree->SetBranchAddress(CalibPrefix + "amplitude" + CalibSuffix, &amplitude , NULL);
569  std::vector<double>* gainused = 0; tree->SetBranchAddress(CalibPrefix + "gainused" + CalibSuffix, &gainused , NULL);
570 
571 
572  printf("Number of Events = %i + %i = %i\n",NEvent,(unsigned int)tree->GetEntries(),(unsigned int)(NEvent+tree->GetEntries()));
573  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
574  printf("Looping on the Tree :");
575  int TreeStep = tree->GetEntries()/50;if(TreeStep<=1)TreeStep=1;
576  for (unsigned int ientry = 0; ientry < tree->GetEntries(); ientry++) {
577  if(ientry%TreeStep==0){printf(".");fflush(stdout);}
578  tree->GetEntry(ientry);
579 
580  if(runnumber<SRun)SRun=runnumber;
581  if(runnumber>ERun)ERun=runnumber;
582 
583  NEvent++;
584  NTrack+=(*trackp).size();
585 
586  unsigned int FirstAmplitude=0;
587  for(unsigned int i=0;i<(*chargeoverpath).size();i++){
588  FirstAmplitude+=(*nstrips)[i];
589  int TI = (*trackindex)[i];
590  //printf("%i - %i - %i %i %i\n", (int)(*rawid)[i], (int)(*firststrip)[i]/128, (int)(*farfromedge)[i], (int)(*overlapping)[i], (int)(*saturation )[i] );
591  if((*tracketa )[TI] < MinTrackEta )continue;
592  if((*tracketa )[TI] > MaxTrackEta )continue;
593  if((*trackp )[TI] < MinTrackMomentum )continue;
594  if((*trackp )[TI] > MaxTrackMomentum )continue;
595  if((*trackhitsvalid)[TI] < MinTrackHits )continue;
596  if((*trackchi2ndof )[TI] > MaxTrackChiOverNdf )continue;
597 
598  stAPVGain* APV = APVsColl[((*rawid)[i]<<4) | ((*firststrip)[i]/128)]; //works for both strip and pixel thanks to firstStrip encoding for pixel in the calibTree
599 
600  if(APV->SubDet>2 && (*farfromedge)[i] == false )continue;
601  if(APV->SubDet>2 && (*overlapping)[i] == true )continue;
602  if(APV->SubDet>2 && (*saturation )[i] && !AllowSaturation)continue;
603  if(APV->SubDet>2 && (*nstrips )[i] > MaxNrStrips )continue;
604 
605 
606  //printf("detId=%7i run=%7i event=%9i charge=%5i cs=%3i\n",(*rawid)[i],runnumber,eventnumber,(*charge)[i],(*nstrips)[i]);
607 
608  //double trans = atan2((*localdiry)[i],(*localdirx)[i])*(180/3.14159265);
609  //double alpha = acos ((*localdirx)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
610  //double beta = acos ((*localdiry)[i] / sqrt( pow((*localdirx)[i],2) + pow((*localdirz)[i],2) ) ) * (180/3.14159265);
611 
612  //printf("NStrip = %i : Charge = %i --> Path = %f --> ChargeOverPath=%f\n",(*nstrips)[i],(*charge)[i],(*path)[i],(*chargeoverpath)[i]);
613  //printf("Amplitudes: ");
614  //for(unsigned int a=0;a<(*nstrips)[i];a++){printf("%i ",(*amplitude)[FirstAmplitude+a]);}
615  //printf("\n");
616 
617  if(APV->SubDet>2){NClusterStrip++;}else{NClusterPixel++;}
618 
619  int Charge = 0;
620  if(APV->SubDet>2 && (useCalibration || !FirstSetOfConstants)){
621  bool Saturation = false;
622  for(unsigned int s=0;s<(*nstrips)[i];s++){
623  int StripCharge = (*amplitude)[FirstAmplitude-(*nstrips)[i]+s];
624  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
625  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
626  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
627  if(StripCharge>1024){
628  StripCharge = 255;
629  Saturation = true;
630  }else if(StripCharge>254){
631  StripCharge = 254;
632  Saturation = true;
633  }
634  Charge += StripCharge;
635  }
636  if(Saturation && !AllowSaturation)continue;
637  }else if(APV->SubDet>2){
638  Charge = (*charge)[i];
639  }else{
640  Charge = (*charge)[i]/265.0; //expected scale factor between pixel and strip charge
641  }
642 
643  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],Charge,APV->CalibGain);
644 
645  double ClusterChargeOverPath = ( (double) Charge )/(*path)[i] ;
646  if(APV->SubDet>2){
647  if(Validation) {ClusterChargeOverPath/=(*gainused)[i];}
648  if(OldGainRemoving){ClusterChargeOverPath*=(*gainused)[i];}
649  }
650  Charge_Vs_Index_Absolute->Fill(APV->Index,Charge);
651  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
652 
653  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill((*path)[i],Charge);
654  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill((*path)[i],Charge);
655  }else if(APV->SubDet==StripSubdetector::TID){
656  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill((*path)[i],Charge);
657  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill((*path)[i],Charge);
658  }
659  }else if(APV->SubDet==StripSubdetector::TEC){
660  if(APV->Eta<0){
661  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill((*path)[i],Charge);
662  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill((*path)[i],Charge);
663  }
664  }else if(APV->Eta>0){
665  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill((*path)[i],Charge);
666  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill((*path)[i],Charge);
667  }
668  }
669  }
670 
671  }// END OF ON-CLUSTER LOOP
672  }printf("\n");// END OF EVENT LOOP
673 
674  }
675 }
676 
677 
678 
680  unsigned int I=0;
681  TH1F* Proj = NULL;
682  double FitResults[6];
683  double MPVmean = 300;
684 
685  TH2F *chvsidx = Charge_Vs_Index->getTH2F();
686 
687 
688  printf("Progressing Bar :0%% 20%% 40%% 60%% 80%% 100%%\n");
689  printf("Fitting Charge Distribution :");
690  int TreeStep = APVsColl.size()/50;
691  for(__gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator it = APVsColl.begin();it!=APVsColl.end();it++,I++){
692  if(I%TreeStep==0){printf(".");fflush(stdout);}
693  stAPVGain* APV = it->second;
694  if(APV->Bin<0) APV->Bin = chvsidx->GetXaxis()->FindBin(APV->Index);
695 
696  if(APV->isMasked){APV->Gain=APV->PreviousGain; MASKED++; continue;}
697 
698  Proj = (TH1F*)(chvsidx->ProjectionY("",chvsidx->GetXaxis()->FindBin(APV->Index),chvsidx->GetXaxis()->FindBin(APV->Index),"e"));
699  if(!Proj)continue;
700 
701  if(CalibrationLevel==0){
702  }else if(CalibrationLevel==1){
703  int SecondAPVId = APV->APVId;
704  if(SecondAPVId%2==0){ SecondAPVId = SecondAPVId+1; }else{ SecondAPVId = SecondAPVId-1; }
705  stAPVGain* APV2 = APVsColl[(APV->DetId<<4) | SecondAPVId];
706  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
707  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
708  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
709  }else if(CalibrationLevel==2){
710  for(unsigned int i=0;i<16;i++){ //loop up to 6APV for Strip and up to 16 for Pixels
711  __gnu_cxx::hash_map<unsigned int, stAPVGain*, __gnu_cxx::hash<unsigned int>, isEqual >::iterator tmpit;
712  tmpit = APVsColl.find((APV->DetId<<4) | i);
713  if(tmpit==APVsColl.end())continue;
714  stAPVGain* APV2 = tmpit->second;
715  if(APV2->DetId != APV->DetId || APV2->APVId == APV->APVId)continue;
716  if(APV2->Bin<0) APV2->Bin = chvsidx->GetXaxis()->FindBin(APV2->Index);
717  TH1F* Proj2 = (TH1F*)(chvsidx->ProjectionY("",APV2->Bin,APV2->Bin,"e"));
718  if(Proj2){Proj->Add(Proj2,1);delete Proj2;}
719  }
720  }else{
721  CalibrationLevel = 0;
722  printf("Unknown Calibration Level, will assume %i\n",CalibrationLevel);
723  }
724 
725  getPeakOfLandau(Proj,FitResults);
726  APV->FitMPV = FitResults[0];
727  APV->FitMPVErr = FitResults[1];
728  APV->FitWidth = FitResults[2];
729  APV->FitWidthErr = FitResults[3];
730  APV->FitChi2 = FitResults[4];
731  APV->FitNorm = FitResults[5];
732  APV->NEntries = Proj->GetEntries();
733 
734  if(IsGoodLandauFit(FitResults)){
735  APV->Gain = APV->FitMPV / MPVmean;
736  if(APV->SubDet>2)GOOD++;
737  }else{
738  APV->Gain = APV->PreviousGain;
739  if(APV->SubDet>2)BAD++;
740  }
741  if(APV->Gain<=0) APV->Gain = 1;
742 
743  //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);
744  delete Proj;
745  }printf("\n");
746 }
747 
748 
750 {
751  unsigned int tree_Index;
752  unsigned int tree_Bin;
753  unsigned int tree_DetId;
754  unsigned char tree_APVId;
755  unsigned char tree_SubDet;
756  float tree_x;
757  float tree_y;
758  float tree_z;
759  float tree_Eta;
760  float tree_R;
761  float tree_Phi;
762  float tree_Thickness;
763  float tree_FitMPV;
764  float tree_FitMPVErr;
765  float tree_FitWidth;
766  float tree_FitWidthErr;
767  float tree_FitChi2NDF;
768  float tree_FitNorm;
769  double tree_Gain;
770  double tree_PrevGain;
771  double tree_PrevGainTick;
772  double tree_NEntries;
773  bool tree_isMasked;
774 
775  TTree* MyTree;
776  MyTree = tfs->make<TTree> ("APVGain","APVGain");
777  MyTree->Branch("Index" ,&tree_Index ,"Index/i");
778  MyTree->Branch("Bin" ,&tree_Bin ,"Bin/i");
779  MyTree->Branch("DetId" ,&tree_DetId ,"DetId/i");
780  MyTree->Branch("APVId" ,&tree_APVId ,"APVId/b");
781  MyTree->Branch("SubDet" ,&tree_SubDet ,"SubDet/b");
782  MyTree->Branch("x" ,&tree_x ,"x/F");
783  MyTree->Branch("y" ,&tree_y ,"y/F");
784  MyTree->Branch("z" ,&tree_z ,"z/F");
785  MyTree->Branch("Eta" ,&tree_Eta ,"Eta/F");
786  MyTree->Branch("R" ,&tree_R ,"R/F");
787  MyTree->Branch("Phi" ,&tree_Phi ,"Phi/F");
788  MyTree->Branch("Thickness" ,&tree_Thickness ,"Thickness/F");
789  MyTree->Branch("FitMPV" ,&tree_FitMPV ,"FitMPV/F");
790  MyTree->Branch("FitMPVErr" ,&tree_FitMPVErr ,"FitMPVErr/F");
791  MyTree->Branch("FitWidth" ,&tree_FitWidth ,"FitWidth/F");
792  MyTree->Branch("FitWidthErr" ,&tree_FitWidthErr,"FitWidthErr/F");
793  MyTree->Branch("FitChi2NDF" ,&tree_FitChi2NDF ,"FitChi2NDF/F");
794  MyTree->Branch("FitNorm" ,&tree_FitNorm ,"FitNorm/F");
795  MyTree->Branch("Gain" ,&tree_Gain ,"Gain/D");
796  MyTree->Branch("PrevGain" ,&tree_PrevGain ,"PrevGain/D");
797  MyTree->Branch("PrevGainTick" ,&tree_PrevGainTick,"PrevGainTick/D");
798  MyTree->Branch("NEntries" ,&tree_NEntries ,"NEntries/D");
799  MyTree->Branch("isMasked" ,&tree_isMasked ,"isMasked/O");
800 
801 
802  FILE* Gains = stdout;
803  fprintf(Gains,"NEvents = %i\n",NEvent);
804  fprintf(Gains,"NTracks = %i\n",NTrack);
805  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
806  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
807  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
808  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
809  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
810 
811  Gains=fopen(OutputGains.c_str(),"w");
812  fprintf(Gains,"NEvents = %i\n",NEvent);
813  fprintf(Gains,"NTracks = %i\n",NTrack);
814  fprintf(Gains,"NClustersPixel = %i\n",NClusterPixel);
815  fprintf(Gains,"NClustersStrip = %i\n",NClusterStrip);
816  fprintf(Gains,"Number of Strip APVs = %lu\n",static_cast<unsigned long>(NStripAPVs));
817  fprintf(Gains,"Number of Pixel Dets = %lu\n",static_cast<unsigned long>(NPixelDets));
818  fprintf(Gains,"GoodFits = %i BadFits = %i ratio = %f%% (MASKED=%i)\n",GOOD,BAD,(100.0*GOOD)/(GOOD+BAD), MASKED);
819 
820  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
822  if(APV==NULL)continue;
823 // printf( "%i | %i | PreviousGain = %7.5f NewGain = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGain,APV->Gain, APV->NEntries);
824  fprintf(Gains,"%i | %i | PreviousGain = %7.5f(tick) x %7.5f(particle) NewGain (particle) = %7.5f (#clusters=%8.0f)\n", APV->DetId,APV->APVId,APV->PreviousGainTick, APV->PreviousGain,APV->Gain, APV->NEntries);
825 
826  tree_Index = APV->Index;
827  tree_Bin = Charge_Vs_Index->getTH2F()->GetXaxis()->FindBin(APV->Index);
828  tree_DetId = APV->DetId;
829  tree_APVId = APV->APVId;
830  tree_SubDet = APV->SubDet;
831  tree_x = APV->x;
832  tree_y = APV->y;
833  tree_z = APV->z;
834  tree_Eta = APV->Eta;
835  tree_R = APV->R;
836  tree_Phi = APV->Phi;
837  tree_Thickness = APV->Thickness;
838  tree_FitMPV = APV->FitMPV;
839  tree_FitMPVErr = APV->FitMPVErr;
840  tree_FitWidth = APV->FitWidth;
841  tree_FitWidthErr= APV->FitWidthErr;
842  tree_FitChi2NDF = APV->FitChi2;
843  tree_FitNorm = APV->FitNorm;
844  tree_Gain = APV->Gain;
845  tree_PrevGain = APV->PreviousGain;
846  tree_PrevGainTick = APV->PreviousGainTick;
847  tree_NEntries = APV->NEntries;
848  tree_isMasked = APV->isMasked;
849 
850 
851  if(tree_DetId==402673324){
852  printf("%i | %i : %f --> %f (%f)\n", tree_DetId, tree_APVId, tree_PrevGain, tree_Gain, tree_NEntries);
853  }
854 
855 
856  MyTree->Fill();
857  }
858  if(Gains)fclose(Gains);
859 
860 
861 }
862 
864 
865  // 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
866  if(Charge_Vs_Index->getTH2F()->Integral(0,NStripAPVs+1, 0, 99999 ) < tagCondition_NClusters) {
867  edm::LogWarning("SiStripGainFromCalibTree")<< "produceTagFilter -> Return false: Statistics is too low : " << Charge_Vs_Index->getTH2F()->Integral() << endl;
868  return false;
869  }
870  if((1.0 * GOOD) / (GOOD+BAD) < tagCondition_GoodFrac) {
871  edm::LogWarning("SiStripGainFromCalibTree")<< "produceTagFilter -> Return false: ratio of GOOD/TOTAL is too low: " << (1.0 * GOOD) / (GOOD+BAD) << endl;
872  return false;
873  }
874  return true;
875 }
876 
878 {
880  if(!harvestingMode) return obj;
881 
882  if(!produceTagFilter()){
883  edm::LogWarning("SiStripGainFromCalibTree")<< "getNewObject -> will not produce a paylaod because produceTagFilter returned false " << endl;
884  setDoStore(false);
885  return obj;
886  }
887 
888 
889  std::vector<float>* theSiStripVector = NULL;
890  unsigned int PreviousDetId = 0;
891  for(unsigned int a=0;a<APVsCollOrdered.size();a++){
893  if(APV==NULL){ printf("Bug\n"); continue; }
894  if(APV->SubDet<=2)continue;
895  if(APV->DetId != PreviousDetId){
896  if(theSiStripVector!=NULL){
897  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
898  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
899  }
900  theSiStripVector = new std::vector<float>;
901  PreviousDetId = APV->DetId;
902  }
903  theSiStripVector->push_back(APV->Gain);
904  }
905  if(theSiStripVector!=NULL){
906  SiStripApvGain::Range range(theSiStripVector->begin(),theSiStripVector->end());
907  if ( !obj->put(PreviousDetId,range) ) printf("Bug to put detId = %i\n",PreviousDetId);
908  }
909 
910  return obj;
911 }
912 
913 
915 {
916 }
917 
919  if(!useCalibration)return;
920 
921  TChain* t1 = new TChain("SiStripCalib/APVGain");
922  t1->Add(m_calibrationPath.c_str());
923 
924  unsigned int tree_DetId;
925  unsigned char tree_APVId;
926  double tree_Gain;
927 
928  t1->SetBranchAddress("DetId" ,&tree_DetId );
929  t1->SetBranchAddress("APVId" ,&tree_APVId );
930  t1->SetBranchAddress("Gain" ,&tree_Gain );
931 
932  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
933  t1->GetEntry(ientry);
934  stAPVGain* APV = APVsColl[(tree_DetId<<4) | (unsigned int)tree_APVId];
935  APV->CalibGain = tree_Gain;
936  }
937 }
938 
939 void
941 {
942  // in AlCaHarvesting mode we just need to run the logic in the endJob step
943  if(harvestingMode) return;
944 
945  if(AlgoMode=="CalibTree")return;
946 
947  if(NEvent==0){
948  SRun = iEvent.id().run();
949  }
950  ERun = iEvent.id().run();
951  NEvent++;
952 
953  //FROM SHALLOW GAIN
954  edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
955  edm::ESHandle<SiStripGain> gainHandle; iSetup.get<SiStripGainRcd>().get(gainHandle);
958 
959  for( TrajTrackAssociationCollection::const_iterator association = associations->begin(); association != associations->end(); association++) {
960  const Trajectory* traj = association->key.get();
961  const reco::Track* track = association->val.get();
962 
963  //clean on the tracks
964  if(fabs(track->eta()) < MinTrackEta )continue;
965  if(fabs(track->eta()) > MaxTrackEta )continue;
966  if(track->p() < MinTrackMomentum )continue;
967  if(track->p() > MaxTrackMomentum )continue;
968  if(track->numberOfValidHits() < MinTrackHits )continue;
969  if(track->chi2()/track->ndof() > MaxTrackChiOverNdf )continue;
970  NTrack++;
971 
972  vector<TrajectoryMeasurement> measurements = traj->measurements();
973  for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
974  TrajectoryStateOnSurface trajState = measurement_it->updatedState();
975  if( !trajState.isValid() ) continue;
976 
977  const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
978  const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
979  const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
980  const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
981  const SiPixelRecHit* sipixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
982 
983  const SiPixelCluster* PixelCluster = NULL;
984  const SiStripCluster* StripCluster = NULL;
985  uint32_t DetId = 0;
986 
987  for(unsigned int h=0;h<2;h++){
988  if(!sistripmatchedhit && h==1){
989  continue;
990  }else if(sistripmatchedhit && h==0){
991  StripCluster = &sistripmatchedhit->monoCluster();
992  DetId = sistripmatchedhit->monoId();
993  }else if(sistripmatchedhit && h==1){
994  StripCluster = &sistripmatchedhit->stereoCluster();;
995  DetId = sistripmatchedhit->stereoId();
996  }else if(sistripsimplehit){
997  StripCluster = (sistripsimplehit->cluster()).get();
998  DetId = sistripsimplehit->geographicalId().rawId();
999  }else if(sistripsimple1dhit){
1000  StripCluster = (sistripsimple1dhit->cluster()).get();
1001  DetId = sistripsimple1dhit->geographicalId().rawId();
1002  }else if(sipixelhit && h==0){
1003  PixelCluster = (sipixelhit->cluster()).get();
1004  DetId = sipixelhit->geographicalId().rawId();
1005  }else{
1006  continue;
1007  }
1008 
1009  LocalVector trackDirection = trajState.localDirection();
1010  double cosine = trackDirection.z()/trackDirection.mag();
1011  int APVId = 0;
1012  bool Saturation = false;
1013  bool Overlapping = false;
1014  unsigned int charge = 0;
1015  double PrevGain = 1;
1016  double Path = 0;
1017 
1018  stAPVGain* APV = NULL;
1019 
1020 
1021 
1022  if(StripCluster){
1023  const auto& Ampls = StripCluster->amplitudes();
1024  int FirstStrip = StripCluster->firstStrip();
1025  APVId = FirstStrip/128;
1026  APV = APVsColl[(DetId<<4) | (FirstStrip/128)];
1027  Path = (10.0*APV->Thickness)/fabs(cosine);
1028 
1029 
1030  if(gainHandle.isValid()){
1031  SiStripApvGain::Range detGainRange = gainHandle->getRange(DetId);
1032  PrevGain = *(detGainRange.first + APVId);
1033  }
1034 
1035  for(unsigned int a=0;a<Ampls.size();a++){
1036  charge+=Ampls[a];
1037  if(Ampls[a] >=254)Saturation =true;
1038  }
1039 
1040  if(FirstStrip==0 )Overlapping=true;
1041  if(FirstStrip==128 )Overlapping=true;
1042  if(FirstStrip==256 )Overlapping=true;
1043  if(FirstStrip==384 )Overlapping=true;
1044  if(FirstStrip==512 )Overlapping=true;
1045  if(FirstStrip==640 )Overlapping=true;
1046 
1047  if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlapping=true;
1048  if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlapping=true;
1049  if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlapping=true;
1050  if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlapping=true;
1051  if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlapping=true;
1052 
1053  if(FirstStrip+Ampls.size()==127 )Overlapping=true;
1054  if(FirstStrip+Ampls.size()==255 )Overlapping=true;
1055  if(FirstStrip+Ampls.size()==383 )Overlapping=true;
1056  if(FirstStrip+Ampls.size()==511 )Overlapping=true;
1057  if(FirstStrip+Ampls.size()==639 )Overlapping=true;
1058  if(FirstStrip+Ampls.size()==767 )Overlapping=true;
1059 
1060  //cleaning on the cluster
1061  if(IsFarFromBorder(&trajState,DetId, &iSetup) == false )continue;
1062  if(Overlapping == true )continue;
1063  if(Saturation && !AllowSaturation )continue;
1064  if(Ampls.size() > MaxNrStrips )continue;
1065  NClusterStrip++;
1066 
1068  bool Saturation = false;
1069  charge = 0;
1070  for(unsigned int s=0;s<Ampls.size();s++){
1071  int StripCharge = Ampls[s];
1072  if(useCalibration && !FirstSetOfConstants){ StripCharge=(int)(StripCharge*(APV->PreviousGain/APV->CalibGain));
1073  }else if(useCalibration){ StripCharge=(int)(StripCharge/APV->CalibGain);
1074  }else if(!FirstSetOfConstants){ StripCharge=(int)(StripCharge*APV->PreviousGain);}
1075  if(StripCharge>1024){
1076  StripCharge = 255;
1077  Saturation = true;
1078  }else if(StripCharge>254){
1079  StripCharge = 254;
1080  Saturation = true;
1081  }
1082  charge += StripCharge;
1083  }
1084  if(Saturation && !AllowSaturation)continue;
1085  }
1086  }else if(PixelCluster){
1087  const auto& Ampls = PixelCluster->pixelADC();
1088  int FirstRow = PixelCluster->minPixelRow();
1089  int FirstCol = PixelCluster->minPixelCol();
1090  APVId = ((FirstRow/80)<<3 | (FirstCol/52));
1091  APV = APVsColl[(DetId<<4) | APVId];
1092  Path = (10.0*APV->Thickness)/fabs(cosine);
1093 
1094  for(unsigned int a=0;a<Ampls.size();a++){
1095  charge+=Ampls[a];
1096  if(Ampls[a] >=254)Saturation =true;
1097  }
1098  charge/=265.0; //expected scale factor between pixel and strip charge
1099  }
1100 
1101  //printf("ChargeDifference = %i Vs %i with Gain = %f\n",(*charge)[i],charge,APV->CalibGain);
1102 
1103  double ClusterChargeOverPath = ( (double) charge )/Path ;
1104  if(Validation) {ClusterChargeOverPath/=PrevGain;}
1105  if(OldGainRemoving){ClusterChargeOverPath*=PrevGain;}
1106  Charge_Vs_Index_Absolute->Fill(APV->Index,charge);
1107  Charge_Vs_Index ->Fill(APV->Index,ClusterChargeOverPath);
1108 
1109 
1110  if(APV->SubDet==StripSubdetector::TIB){ Charge_Vs_PathlengthTIB ->Fill(Path,charge);
1111  }else if(APV->SubDet==StripSubdetector::TOB){ Charge_Vs_PathlengthTOB ->Fill(Path,charge);
1112  }else if(APV->SubDet==StripSubdetector::TID){
1113  if(APV->Eta<0){ Charge_Vs_PathlengthTIDM ->Fill(Path,charge);
1114  }else if(APV->Eta>0){ Charge_Vs_PathlengthTIDP ->Fill(Path,charge);
1115  }
1116  }else if(APV->SubDet==StripSubdetector::TEC){
1117  if(APV->Eta<0){
1118  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECM1->Fill(Path,charge);
1119  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECM2->Fill(Path,charge);
1120  }
1121  }else if(APV->Eta>0){
1122  if(APV->Thickness<0.04){ Charge_Vs_PathlengthTECP1->Fill(Path,charge);
1123  }else if(APV->Thickness>0.04){ Charge_Vs_PathlengthTECP2->Fill(Path,charge);
1124  }
1125  }
1126  }
1127 
1128  }//loop on clusters
1129  }//loop on measurements
1130  }//loop on tracks
1131 
1132 }
1133 
1134 
1136 {
1137  edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
1138 
1139  LocalPoint HitLocalPos = trajState->localPosition();
1140  LocalError HitLocalError = trajState->localError().positionError() ;
1141 
1142  const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
1143  if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
1144  edm::LogWarning("SiStripGainFromCalibTree")<< "this detID doesn't seem to belong to the Tracker" << std::endl;
1145  return false;
1146  }
1147 
1148  const BoundPlane plane = it->surface();
1149  const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
1150  const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
1151 
1152  double DistFromBorder = 1.0;
1153  double HalfLength = it->surface().bounds().length() /2.0;
1154 
1155  if(trapezoidalBounds)
1156  {
1157  std::array<const float, 4> const & parameters = (*trapezoidalBounds).parameters();
1158  //std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
1159  HalfLength = parameters[3];
1160  }else if(rectangularBounds){
1161  HalfLength = it->surface().bounds().length() /2.0;
1162  }else{return false;}
1163 
1164  if (fabs(HitLocalPos.y())+HitLocalError.yy() >= (HalfLength - DistFromBorder) ) return false;
1165 
1166  return true;
1167 }
1168 
1169 
1170 
RunNumber_t run() const
Definition: EventID.h:39
ClusterRef cluster() const
double p() const
momentum vector magnitude
Definition: TrackBase.h:602
bool IsFarFromBorder(const TrajectoryStateOnSurface &trajState, const GeomDetUnit *it)
Definition: DeDxTools.cc:311
T getParameter(std::string const &) const
virtual int nstrips() const =0
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
dictionary parameters
Definition: Parameters.py:2
int minPixelCol() const
unsigned int stereoId() const
virtual float length() const =0
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventIDconst &, edm::Timestampconst & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
SiStripApvGain * getNewObject() override
virtual int ncolumns() const =0
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
virtual void algoBeginRun(const edm::Run &run, const edm::EventSetup &iSetup) override
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
SiStripCluster const & monoCluster() const
LocalVector localDirection() const
virtual void algoEndRun(const edm::Run &run, const edm::EventSetup &iSetup) override
SiStripGainFromCalibTree(const edm::ParameterSet &)
T y() const
Definition: PV3DBase.h:63
virtual int nrows() const =0
#define NULL
Definition: scimark2.h:8
const Bounds & bounds() const
Definition: Surface.h:128
void storeOnTree(TFileService *tfs)
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
uint16_t firstStrip() const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
LocalError positionError() const
void Fill(long long x)
bool put(const uint32_t &detID, Range input)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
T x() const
Cartesian x coordinate.
std::vector< stAPVGain * > APVsCollOrdered
DataContainer const & measurements() const
Definition: Trajectory.h:203
__gnu_cxx::hash_map< unsigned int, stAPVGain *, __gnu_cxx::hash< unsigned int >, isEqual > APVsColl
int iEvent
Definition: GenABIO.cc:230
T mag() const
Definition: PV3DBase.h:67
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:638
int minPixelRow() const
float yy() const
Definition: LocalError.h:26
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:536
bool IsGoodLandauFit(double *FitResults)
double ndof() const
number of degrees of freedom of the fit
Definition: TrackBase.h:542
Definition: Path.h:40
T z() const
Definition: PV3DBase.h:64
edm::EDGetTokenT< edm::View< reco::Track > > tracksToken
virtual void algoBeginJob(const edm::EventSetup &iSetup) override
std::pair< ContainerIterator, ContainerIterator > Range
int j
Definition: DBlmapReader.cc:9
edm::EDGetTokenT< TrajTrackAssociationCollection > associationsToken
const std::complex< double > I
Definition: I.h:8
ClusterRef cluster() const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:807
const LocalTrajectoryError & localError() const
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
const std::vector< uint16_t > & pixelADC() const
virtual void algoEndJob() override
Definition: DetId.h:18
virtual void algoAnalyze(const edm::Event &, const edm::EventSetup &) override
bool IsFarFromBorder(TrajectoryStateOnSurface *trajState, const uint32_t detid, const edm::EventSetup *iSetup)
tuple tracks
Definition: testEve_cfg.py:39
bool merge(LuminosityBlockRange &lh, LuminosityBlockRange &rh)
void getPeakOfLandau(TH1 *InputHisto, double *FitResults, double LowRange=50, double HighRange=5400)
ClusterRef cluster() const
Definition: SiPixelRecHit.h:49
const T & get() const
Definition: EventSetup.h:56
void setDoStore(const bool doStore)
When set to false the payload will not be written to the db.
SiStripCluster const & stereoCluster() const
edm::EventID id() const
Definition: EventBase.h:60
Pixel cluster – collection of neighboring pixels above threshold.
double a
Definition: hdecay.h:121
TH2F * getTH2F(std::string name, std::string process, std::string rootfolder, DQMStore *dbe_, bool verb, bool clone)
unsigned int monoId() const
DetId geographicalId() const
volatile std::atomic< bool > shutdown_flag false
TH2F * getTH2F(void) const
bool isValid() const
Definition: ESHandle.h:47
long double T
const std::vector< uint8_t > & amplitudes() const
Definition: Run.h:43
Our base class.
Definition: SiPixelRecHit.h:23