CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuScleFit.cc
Go to the documentation of this file.
1 // \class MuScleFit
2 // Fitter of momentum scale and resolution from resonance decays to muon track pairs
3 //
4 // $Date: 2010/10/19 16:36:51 $
5 // $Revision: 1.96 $
6 // \author R. Bellan, C.Mariotti, S.Bolognesi - INFN Torino / T.Dorigo, M.De Mattia - INFN Padova
7 //
8 // Recent additions:
9 // - several parameters allow a more flexible use, tests, and control handles
10 // for the likelihood. In particular, a set of integers controls the order
11 // with which parameters are released; another controls which parameters are
12 // fixed. A function allows to smear momenta and angles of muons from the
13 // resonance before any correction, using a set of random numbers generated
14 // once and for all in the constructor (this way the smearing remains the same
15 // for any given muon no matter how many times one loops and what corrections
16 // one applies).
17 // For a correct use of these flags, please see the function minimizeLikelihood() in
18 // MuScleFitUtils.cc
19 // - the fit now allows to extract resolution functions simultaneously with the
20 // momentum scale. So far a simple parametrization
21 // of muon momentum resolution and angle resolution has been implemented, but
22 // extensions are straightforward.
23 // - It is however advisable to fit separately resolution and scale. The suggested
24 // course of action is:
25 // 1) fit the scale with a simple parametrization
26 // 2) check results, fit with more complicated forms
27 // 3) verify which is a sufficiently accurate description of the data
28 // 4) fix scale parameters and fit for the resolution
29 // 5) go back to fitting the scale with resolution parameters fixed to fitted values
30 // - Also note that resolution fits may fail to converge due to instability of the
31 // probability distribution to the limit of large widths. Improvements here are
32 // advisable.
33 // - The treatment of signal windows in the Y region
34 // has to be refined because of overlaps. More work is needed here, assigning a different
35 // weight to different hypothesis of the resonance producing a given mass, if there are
36 // multiple candidates.
37 // - Also, larger windows are to be allowed for fits to SA muons.
38 // - File Probs_1000.root contains the probability distribution of lorentzians convoluted
39 // with gaussian smearing functions, for the six resonances. A 1000x1000 grid
40 // in mass,sigma has been computed (using root macro Probs.C).
41 // A wider interval of masses for each resonance should be computed, to be used for standalone muons
42 //
43 //
44 // Notes on additions, TD 31/3/08
45 //
46 // - background model: at least a couple of different models, with two parameters,
47 // should be included in the fitting procedure such that the function massprob(),
48 // which produces the probability in the likelihood computation, incorporates the
49 // probability that the event is from background. That is, if the fitting function
50 // knows the shape of the mass spectrum ( say, a falling exponential plus a gaussian
51 // signal) it becomes possible to fit the scale together with the background shape
52 // and normalization parameters. Of course, one should do one thing at a time: first
53 // a scale fit, then a shape fit with scale parameters fixed, and then a combined
54 // fit. Convergence issues should be handled case by case.
55 // - The correct implementation of the above idea requires a reorganization of pass
56 // parameters (in the cfg) and fit parameters. The user has to be able to smear,
57 // bias, fix parameters, choose scale fitting functions, resolution fitting functions,
58 // and background functions. It should be possible to separate the fit functions from
59 // the biasing ones, which would allow a more thorough testing.
60 // - all the above can be obtained by making the .cfg instructions heavier. Since this
61 // is a routine operated by experts only, it is a sensible choice.
62 // - One should thus envision the following:
63 // 1) a set of parameters controlling the biasing function: parBias()
64 // 2) a set of parameters controlling the smearing function: parSmear()
65 // 3) a set of parameters to define resolution modeling and initial values: parResol()
66 // 3b) parResol() gets fix and order bits by parResolFix() and parResolOrder()
67 // 4) a set of parameters to define scale modeling and initial values: parScale()
68 // 4b) parScale() gets fix and order bits by parScaleFix() and parScaleOrder()
69 // 5) a set of parameters controlling the background shape and normalization: parNorm()
70 // 5b) parNorm() gets fix and order bits by parNormFix() and parNormOrder()
71 // The likelihood parameters then become a vector which is dynamically composed of
72 // sets 3), 4), and 5): parval() = parResol()+parScale()+parNorm()
73 // - In order to study better the likelihood behavior it would be advisable to introduce
74 // some histogram filling on the last iteration of the likelihood. It is not clear
75 // how best to achieve that: probably the simplest way is to make a histogram filling
76 // function run just after the likelihood computation, such that the best value of the
77 // fit parameters is used.
78 // - The muon pair which we call our resonance must be chosen in a way which does not
79 // bias our likelihood: we cannot just choose the pair closest to a resonance.
80 //
81 //
82 // Notes on additions, T.Dorigo 22/12/2008
83 // ---------------------------------------
84 //
85 // - File Probs_new_1000_CTEQ.root now contains a set of 24 additional two-dim histograms,
86 // defining the probability distribution of Z boson decays as a function of measured mass
87 // and expected sigma in 24 different bins of Z rapidity, extracted from CTEQ 6 PDF (at
88 // Leading Order) from the convolution in the factorization integral. See programs CTEQ.cpp
89 // and Fits.C.
90 // - The probability for Z boson events now thus depends on the measured rapidity of the dimuon
91 // system. All functions in file MuScleFitUtils.cc have been suitably changed.
92 //
93 // ----------------------------------------------------------------------------------
94 // Modifications by M. De Mattia 13/3/2009
95 // ---------------------------------------
96 // - The histograms map was moved to a base class (MuScleFitBase) from which this one inherits.
97 //
98 // Modifications by M. De Mattia 20/7/2009
99 // ---------------------------------------
100 // - Reworked background fit based on ranges. See comments in the code for more details.
101 // ---------------------------------------------------------------------------------------------
102 
103 #include "MuScleFit.h"
110 
111 #include "HepPDT/defs.h"
112 #include "HepPDT/TableBuilder.hh"
113 #include "HepPDT/ParticleDataTable.hh"
114 
115 #include "HepMC/GenParticle.h"
116 #include "HepMC/GenEvent.h"
117 // #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
119 
123 
124 #include "TFile.h"
125 #include "TTree.h"
126 #include "TMinuit.h"
127 #include <vector>
128 
131 
132 // To use callgrind for code profiling uncomment also the following define.
133 // #define USE_CALLGRIND
134 
135 #ifdef USE_CALLGRIND
136 #include "valgrind/callgrind.h"
137 #endif
138 
141 
142 // To read likelihood distributions from the database.
143  //#include "CondFormats/RecoMuonObjects/interface/MuScleFitLikelihoodPdf.h"
144  //#include "CondFormats/DataRecord/interface/MuScleFitLikelihoodPdfRcd.h"
145 
146 // Constructor
147 // -----------
149  MuScleFitBase( pset ),
150  totalEvents_(0)
151 {
153  if (debug_>0) std::cout << "[MuScleFit]: Constructor" << std::endl;
154 
155  if ((theMuonType_<-4 || theMuonType_>5) && theMuonType_<10) {
156  std::cout << "[MuScleFit]: Unknown muon type! Aborting." << std::endl;
157  abort();
158  }
159 
160  loopCounter = 0;
161 
162  // Boundaries for h-function computation (to be improved!)
163  // -------------------------------------------------------
164  minResMass_hwindow[0] = 76.;
165  maxResMass_hwindow[0] = 106.;
166  minResMass_hwindow[1] = 10.15;
167  maxResMass_hwindow[1] = 10.55;
168  minResMass_hwindow[2] = 9.8;
169  maxResMass_hwindow[2] = 10.2;
170  minResMass_hwindow[3] = 9.25;
171  maxResMass_hwindow[3] = 9.65;
172  minResMass_hwindow[4] = 3.58;
173  maxResMass_hwindow[4] = 3.78;
174  minResMass_hwindow[5] = 3.0;
175  maxResMass_hwindow[5] = 3.2;
176 
177  // Max number of loops (if > 2 then try to minimize likelihood more than once)
178  // ---------------------------------------------------------------------------
179  maxLoopNumber = pset.getUntrackedParameter<int>("maxLoopNumber", 2);
180  fastLoop = pset.getUntrackedParameter<bool>("FastLoop", true);
181 
182  // Selection of fits according to loop
183  MuScleFitUtils::doResolFit = pset.getParameter<std::vector<int> >("doResolFit");
184  MuScleFitUtils::doScaleFit = pset.getParameter<std::vector<int> >("doScaleFit");
185  MuScleFitUtils::doCrossSectionFit = pset.getParameter<std::vector<int> >("doCrossSectionFit");
186  MuScleFitUtils::doBackgroundFit = pset.getParameter<std::vector<int> >("doBackgroundFit");
187 
188  // Bias and smear types
189  // --------------------
190  int biasType = pset.getParameter<int>("BiasType");
191  MuScleFitUtils::BiasType = biasType;
192  // No error, the scale functions are used also for the bias
194  int smearType = pset.getParameter<int>("SmearType");
195  MuScleFitUtils::SmearType = smearType;
197 
198  // Fit types
199  // ---------
200  int resolFitType = pset.getParameter<int>("ResolFitType");
201  MuScleFitUtils::ResolFitType = resolFitType;
204  int scaleType = pset.getParameter<int>("ScaleFitType");
205  MuScleFitUtils::ScaleFitType = scaleType;
208 
209  // Initial parameters values
210  // -------------------------
211  MuScleFitUtils::parBias = pset.getParameter<std::vector<double> >("parBias");
212  MuScleFitUtils::parSmear = pset.getParameter<std::vector<double> >("parSmear");
213  MuScleFitUtils::parResol = pset.getParameter<std::vector<double> >("parResol");
214  MuScleFitUtils::parScale = pset.getParameter<std::vector<double> >("parScale");
215  MuScleFitUtils::parCrossSection = pset.getParameter<std::vector<double> >("parCrossSection");
216  MuScleFitUtils::parBgr = pset.getParameter<std::vector<double> >("parBgr");
217  MuScleFitUtils::parResolFix = pset.getParameter<std::vector<int> >("parResolFix");
218  MuScleFitUtils::parScaleFix = pset.getParameter<std::vector<int> >("parScaleFix");
219  MuScleFitUtils::parCrossSectionFix = pset.getParameter<std::vector<int> >("parCrossSectionFix");
220  MuScleFitUtils::parBgrFix = pset.getParameter<std::vector<int> >("parBgrFix");
221  MuScleFitUtils::parResolOrder = pset.getParameter<std::vector<int> >("parResolOrder");
222  MuScleFitUtils::parScaleOrder = pset.getParameter<std::vector<int> >("parScaleOrder");
223  MuScleFitUtils::parCrossSectionOrder = pset.getParameter<std::vector<int> >("parCrossSectionOrder");
224  MuScleFitUtils::parBgrOrder = pset.getParameter<std::vector<int> >("parBgrOrder");
225 
226  MuScleFitUtils::resfind = pset.getParameter<std::vector<int> >("resfind");
227  MuScleFitUtils::FitStrategy = pset.getParameter<int>("FitStrategy");
228 
229  // Option to skip unnecessary stuff
230  // --------------------------------
231  MuScleFitUtils::speedup = pset.getParameter<bool>("speedup");
232 
233  // Option to skip simTracks comparison
234  compareToSimTracks_ = pset.getParameter<bool>("compareToSimTracks");
235  simTracksCollection_ = pset.getUntrackedParameter<edm::InputTag>("SimTracksCollection", edm::InputTag("g4SimHits"));
236 
237  triggerResultsLabel_ = pset.getUntrackedParameter<std::string>("TriggerResultsLabel");
238  triggerResultsProcess_ = pset.getUntrackedParameter<std::string>("TriggerResultsProcess");
239  triggerPath_ = pset.getUntrackedParameter<std::string>("TriggerPath");
240  negateTrigger_ = pset.getUntrackedParameter<bool>("NegateTrigger", false);
241  saveAllToTree_ = pset.getUntrackedParameter<bool>("SaveAllToTree", false);
242 
243  PATmuons_ = pset.getUntrackedParameter<bool>("PATmuons", false);
244  genParticlesName_ = pset.getUntrackedParameter<std::string>("GenParticlesName", "genParticles");
245 
246  // Use the probability file or not. If not it will perform a simpler selection taking the muon pair with
247  // invariant mass closer to the pdf value and will crash if some fit is attempted.
248  MuScleFitUtils::useProbsFile_ = pset.getUntrackedParameter<bool>("UseProbsFile", true);
249 
250  // This must be set to true if using events generated with Sherpa
251  MuScleFitUtils::sherpa_ = pset.getUntrackedParameter<bool>("Sherpa", false);
252 
253  MuScleFitUtils::rapidityBinsForZ_ = pset.getUntrackedParameter<bool>("RapidityBinsForZ", true);
254 
255  // Set the cuts on muons to be used in the fit
256  MuScleFitUtils::maxMuonPt_ = pset.getUntrackedParameter<double>("MaxMuonPt", 100000000.);
257  MuScleFitUtils::minMuonPt_ = pset.getUntrackedParameter<double>("MinMuonPt", 0.);
258  MuScleFitUtils::minMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MinMuonEtaFirstRange", -6.);
259  MuScleFitUtils::maxMuonEtaFirstRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaFirstRange", 6.);
260  MuScleFitUtils::minMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MinMuonEtaSecondRange", -100.);
261  MuScleFitUtils::maxMuonEtaSecondRange_ = pset.getUntrackedParameter<double>("MaxMuonEtaSecondRange", 100.);
262 
263  MuScleFitUtils::debugMassResol_ = pset.getUntrackedParameter<bool>("DebugMassResol", false);
264  // MuScleFitUtils::massResolComponentsStruct MuScleFitUtils::massResolComponents;
265 
266  // Check for parameters consistency
267  // it will abort in case of errors.
268  checkParameters();
269 
270  // Generate array of gaussian-distributed numbers for smearing
271  // -----------------------------------------------------------
272  if (MuScleFitUtils::SmearType>0) {
273  std::cout << "[MuScleFit-Constructor]: Generating random values for smearing" << std::endl;
274  TF1 G("G", "[0]*exp(-0.5*pow(x,2))", -5., 5.);
275  double norm = 1/sqrt(2*TMath::Pi());
276  G.SetParameter (0,norm);
277  for (int i=0; i<10000; i++) {
278  for (int j=0; j<7; j++) {
279  MuScleFitUtils::x[j][i] = G.GetRandom();
280  }
281  }
282  }
284 
285  if(theMuonType_ > 0 && theMuonType_ < 4) {
288  }
289  else if(theMuonType_ == 4 || theMuonType_ >= 10 || theMuonType_==-1 || theMuonType_==-2 || theMuonType_==-3 || theMuonType_==-4) {
292  }
293  else{
294  std::cout<<"Wrong muon type "<<theMuonType_<<std::endl;
295  exit(1);
296  }
297 
298  // When using standalone muons switch to the single Z pdf
299  if( theMuonType_ == 2 ) {
300  MuScleFitUtils::rapidityBinsForZ_ = false;
301  }
302 
303  // Initialize ResMaxSigma And ResHalfWidth - 0 = global, 1 = SM, 2 = tracker
304  // -------------------------------------------------------------------------
323 
325  MuScleFitUtils::resfind,
326  MuScleFitUtils::speedup, genParticlesName_,
327  compareToSimTracks_, simTracksCollection_,
328  MuScleFitUtils::sherpa_, debug_));
329 
330  MuScleFitUtils::backgroundHandler = new BackgroundHandler( pset.getParameter<std::vector<int> >("BgrFitType"),
331  pset.getParameter<std::vector<double> >("LeftWindowBorder"),
332  pset.getParameter<std::vector<double> >("RightWindowBorder"),
335 
336  MuScleFitUtils::crossSectionHandler = new CrossSectionHandler( MuScleFitUtils::parCrossSection, MuScleFitUtils::resfind );
337 
338  // Build cross section scale factors
339  // MuScleFitUtils::resfind
340 
341  MuScleFitUtils::normalizeLikelihoodByEventNumber_ = pset.getUntrackedParameter<bool>("NormalizeLikelihoodByEventNumber", true);
342  if(debug_>0) std::cout << "End of MuScleFit constructor" << std::endl;
343 
344  inputRootTreeFileName_ = pset.getParameter<std::string>("InputRootTreeFileName");
345  outputRootTreeFileName_ = pset.getParameter<std::string>("OutputRootTreeFileName");
346  maxEventsFromRootTree_ = pset.getParameter<int>("MaxEventsFromRootTree");
347 
348  MuScleFitUtils::startWithSimplex_ = pset.getParameter<bool>("StartWithSimplex");
349  MuScleFitUtils::computeMinosErrors_ = pset.getParameter<bool>("ComputeMinosErrors");
350  MuScleFitUtils::minimumShapePlots_ = pset.getParameter<bool>("MinimumShapePlots");
351 
353 }
354 
355 // Destructor
356 // ----------
358  if (debug_>0) std::cout << "[MuScleFit]: Destructor" << std::endl;
359  std::cout << "Total number of analyzed events = " << totalEvents_ << std::endl;
360 
361  if( !(outputRootTreeFileName_.empty()) ) {
362  // Save the events to a root tree unless we are reading from the edm root file and the SavedPair size is different from the totalEvents_
363  if( !(inputRootTreeFileName_.empty() && (int(MuScleFitUtils::SavedPair.size()) != totalEvents_)) ) {
364  std::cout << "Saving muon pairs to root tree" << std::endl;
365  RootTreeHandler rootTreeHandler;
367  // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, 0, saveAllToTree_);
369  }
370  else {
371  // rootTreeHandler.writeTree(outputRootTreeFileName_, &(MuScleFitUtils::SavedPair), theMuonType_, &(MuScleFitUtils::genPair), saveAllToTree_ );
373  }
374  }
375  else {
376  std::cout << "ERROR: events in the vector = " << MuScleFitUtils::SavedPair.size() << " != totalEvents = " << totalEvents_ << std::endl;
377  }
378  }
379 }
380 
381 // Begin job
382 // ---------
384 // void MuScleFit::beginOfJob ()
385 // void MuScleFit::beginOfJob( const edm::EventSetup& eventSetup )
386 {
387  if (debug_>0) std::cout << "[MuScleFit]: beginOfJob" << std::endl;
388  //if(maxLoopNumber>1)
391  }
392 
393  if (debug_>0) std::cout << "[MuScleFit]: beginOfJob" << std::endl;
394 
395  // Create the root file
396  // --------------------
397  for (unsigned int i=0; i<(maxLoopNumber); i++) {
398  std::stringstream ss;
399  ss << i;
400  std::string rootFileName = ss.str() + "_" + theRootFileName_;
401  theFiles_.push_back (new TFile(rootFileName.c_str(), "RECREATE"));
402  }
403  if (debug_>0) std::cout << "[MuScleFit]: Root file created" << std::endl;
404 
405  std::cout << "creating plotter" << std::endl;
407  plotter->debug = debug_;
408 }
409 
410 // End of job method
411 // -----------------
413  if (debug_>0) std::cout << "[MuScleFit]: endOfJob" << std::endl;
414 }
415 
416 // New loop
417 // --------
418 void MuScleFit::startingNewLoop( unsigned int iLoop )
419 {
420  if (debug_>0) std::cout << "[MuScleFit]: Starting loop # " << iLoop << std::endl;
421 
422  // Number of muons used
423  // --------------------
425 
426  // Counters for problem std::cout-ing
427  // -----------------------------
429 
430  // Create the root file
431  // --------------------
432  fillHistoMap(theFiles_[iLoop], iLoop);
433 
434  loopCounter = iLoop;
436 
437  iev = 0;
439 
441 }
442 
443 // End of loop routine
444 // -------------------
445 edm::EDLooper::Status MuScleFit::endOfLoop( const edm::EventSetup& eventSetup, unsigned int iLoop )
446 {
447  unsigned int iFastLoop = 1;
448 
449  // Read the events from the root tree if requested
450  if( !(inputRootTreeFileName_.empty()) ) {
452  // When reading from local file all the loops are done here
454  iFastLoop = 0;
455  }
456  else {
457  endOfFastLoop(iLoop);
458  }
459 
460  // If a fastLoop is required we do all the remaining iterations here
461  if( fastLoop == true ) {
462  for( ; iFastLoop<maxLoopNumber; ++iFastLoop ) {
463 
464  std::cout << "Starting fast loop number " << iFastLoop << std::endl;
465 
466  // In the first loop is called by the framework
467  // if( iFastLoop > 0 ) {
468  startingNewLoop(iFastLoop);
469  // }
470 
471  // std::vector<std::pair<lorentzVector,lorentzVector> >::const_iterator it = MuScleFitUtils::SavedPair.begin();
472  // for( ; it != SavedPair.end(); ++it ) {
473  while( iev<totalEvents_ ) {
474  if( iev%1000 == 0 ) {
475  std::cout << "Fast looping on event number " << iev << std::endl;
476  }
477  // This reads muons from SavedPair using iev to keep track of the event
478  duringFastLoop();
479  }
480  std::cout << "End of fast loop number " << iFastLoop << ". Ran on " << iev << " events" << std::endl;
481  endOfFastLoop(iFastLoop);
482  }
483  }
484 
485  if (iFastLoop>=maxLoopNumber-1) {
486  return kStop;
487  } else {
488  return kContinue;
489  }
490 }
491 
492 void MuScleFit::endOfFastLoop( const unsigned int iLoop )
493 {
494  // std::cout<< "Inside endOfFastLoop, iLoop = " << iLoop << " and loopCounter = " << loopCounter << std::endl;
495 
496  if( loopCounter == 0 ) {
497  // plotter->writeHistoMap();
498  // The destructor will call the writeHistoMap after the cd to the output file
499  delete plotter;
500  }
501 
502  std::cout << "Ending loop # " << iLoop << std::endl;
503 
504  // Write the histos to file
505  // ------------------------
506  // theFiles_[iLoop]->cd();
507  writeHistoMap(iLoop);
508 
509  // Likelihood minimization to compute corrections
510  // ----------------------------------------------
511  // theFiles_[iLoop]->cd();
512  TDirectory * likelihoodDir = theFiles_[iLoop]->mkdir("likelihood");
513  likelihoodDir->cd();
515 
516  // ATTENTION, this was put BEFORE the minimizeLikelihood. Check for problems.
517  theFiles_[iLoop]->Close();
518  // ATTENTION: Check that this delete does not give any problem
519  delete theFiles_[iLoop];
520 
521  // Clear the histos
522  // ----------------
523  clearHistoMap();
524 }
525 
526 // Stuff to do during loop
527 // -----------------------
529 {
530  edm::Handle<edm::TriggerResults> triggerResults;
531  event.getByLabel(edm::InputTag(triggerResultsLabel_.c_str(), "", triggerResultsProcess_.c_str()), triggerResults);
532  //event.getByLabel(InputTag(triggerResultsLabel_),triggerResults);
533  bool isFired = false;
534 
535  if(triggerPath_ == "")
536  isFired = true;
537  else if(triggerPath_ == "All"){
538  isFired =triggerResults->accept();
539  if(debug_>0)
540  std::cout<<"Trigger "<<isFired<<std::endl;
541  }
542  else{
543  bool changed;
545  hltConfig.init(event.getRun(), eventSetup, triggerResultsProcess_, changed);
546  unsigned int triggerIndex( hltConfig.triggerIndex(triggerPath_) );
547  // triggerIndex must be less than the size of HLTR or you get a CMSException: _M_range_check
548  if (triggerIndex < triggerResults->size()) {
549  isFired = triggerResults->accept(triggerIndex);
550  if(debug_>0)
551  std::cout<<triggerPath_<<" "<<isFired<<std::endl;
552  }
553  }
554 
555  if( negateTrigger_ && isFired ) return kContinue;
556  else if( !(negateTrigger_) && !isFired ) return kContinue;
557 
558 #ifdef USE_CALLGRIND
559  CALLGRIND_START_INSTRUMENTATION;
560 #endif
561 
562  if (debug_>0) {
563  std::cout << "[MuScleFit-duringLoop]: loopCounter = " << loopCounter
564  << " Run: " << event.id().run() << " Event: " << event.id().event() << std::endl;
565  }
566 
567  // On the first iteration we read the bank, otherwise we fetch the information from the muon tree
568  // ------------------------------------ Important Note --------------------------------------- //
569  // The fillMuonCollection method applies any smearing or bias to the muons, so we NEVER use
570  // unbiased muons.
571  // ----------------------------------------------------------------------------------------------
572  if( loopCounter == 0 ) {
573 
574  if( !fastLoop || inputRootTreeFileName_.empty() ) {
575  if( debug_ > 0 ) std::cout << "Reading from edm event" << std::endl;
576  selectMuons(event);
577  duringFastLoop();
578  ++totalEvents_;
579  }
580  }
581 
582  return kContinue;
583 
584 #ifdef USE_CALLGRIND
585  CALLGRIND_STOP_INSTRUMENTATION;
586  CALLGRIND_DUMP_STATS;
587 #endif
588 }
589 
591 {
594 
595  std::vector<reco::LeafCandidate> muons;
596  muonSelector_->selectMuons(event, muons, genMuonPairs_, MuScleFitUtils::simPair, plotter);
597  plotter->fillRec(muons);
598 
599  // Find the two muons from the resonance, and set ResFound bool
600  // ------------------------------------------------------------
601  std::pair<reco::Particle::LorentzVector, reco::Particle::LorentzVector> recMuFromBestRes =
603 
605  if (debug_>0) {
606  std::cout <<std::setprecision(9)<< "Pt after findbestrecores: " << (recMuFromBestRes.first).Pt() << " "
607  << (recMuFromBestRes.second).Pt() << std::endl;
608  std::cout << "recMu1 = " << recMu1 << std::endl;
609  std::cout << "recMu2 = " << recMu2 << std::endl;
610  }
611  recMu1 = recMuFromBestRes.first;
612  recMu2 = recMuFromBestRes.second;
613  if (debug_>0) {
614  std::cout << "after recMu1 = " << recMu1 << std::endl;
615  std::cout << "after recMu2 = " << recMu2 << std::endl;
616  std::cout << "mu1.pt = " << recMu1.Pt() << std::endl;
617  std::cout << "mu2.pt = " << recMu2.Pt() << std::endl;
618  }
619  MuScleFitUtils::SavedPair.push_back( std::make_pair( recMu1, recMu2 ) );
620  } else {
621  MuScleFitUtils::SavedPair.push_back( std::make_pair( lorentzVector(0.,0.,0.,0.), lorentzVector(0.,0.,0.,0.) ) );
622  }
623  // Save the events also in the external tree so that it can be saved later
626  event.run(), event.id().event()));
627  // Fill the internal genPair tree from the external one
628  MuScleFitUtils::genPair.push_back(std::make_pair( genMuonPairs_.back().mu1, genMuonPairs_.back().mu2 ));
629 }
630 
631 void MuScleFit::selectMuons(const int maxEvents, const TString & treeFileName)
632 {
633  std::cout << "Reading muon pairs from Root Tree in " << treeFileName << std::endl;
634  RootTreeHandler rootTreeHandler;
637  }
638  else {
640  }
641  // Now loop on all the pairs and apply any smearing and bias if needed
642  std::vector<std::pair<lorentzVector,lorentzVector> >::iterator it = MuScleFitUtils::SavedPair.begin();
643  for( ; it != MuScleFitUtils::SavedPair.end(); ++it ) {
644 
645  // Apply any cut if requested
646  // Note that cuts here are only applied to already selected muons. They should not be used unless
647  // you are sure that the difference is negligible (e.g. the number of events with > 2 muons is negligible).
648  double pt1 = it->first.pt();
649  // std::cout << "pt1 = " << pt1 << std::endl;
650  double pt2 = it->second.pt();
651  // std::cout << "pt2 = " << pt2 << std::endl;
652  double eta1 = it->first.eta();
653  // std::cout << "eta1 = " << eta1 << std::endl;
654  double eta2 = it->second.eta();
655  // std::cout << "eta2 = " << eta2 << std::endl;
656  // If they don't pass the cuts set to null vectors
663  // std::cout << "removing muons not passing cuts" << std::endl;
664  it->first = reco::Particle::LorentzVector(0,0,0,0);
665  it->second = reco::Particle::LorentzVector(0,0,0,0);
666  }
667 
668  // First is always mu-, second mu+
669  if( (MuScleFitUtils::SmearType != 0) || (MuScleFitUtils::BiasType != 0) ) {
670  applySmearing(it->first);
671  applyBias(it->first, -1);
672  applySmearing(it->second);
673  applyBias(it->second, 1);
674  }
675  }
677  if( !(MuScleFitUtils::speedup) ) {
679  }
680 }
681 
683 {
684  // On loops>0 the two muons are directly obtained from the SavedMuon array
685  // -----------------------------------------------------------------------
686  MuScleFitUtils::ResFound = false;
689  if (recMu1.Pt()>0 && recMu2.Pt()>0) {
691  if (debug_>0) std::cout << "Ev = " << iev << ": found muons in tree with Pt = "
692  << recMu1.Pt() << " " << recMu2.Pt() << std::endl;
693  }
694 
695  if( debug_>0 ) std::cout << "About to start lik par correction and histo filling; ResFound is "
696  << MuScleFitUtils::ResFound << std::endl;
697  // If resonance found, do the hard work
698  // ------------------------------------
700 
701  // Find weight and reference mass for this muon pair
702  // -------------------------------------------------
703  // The last parameter = true means that we want to use always the background window to compute the weight,
704  // otherwise the probability will be filled only for the resonance region.
705  double weight = MuScleFitUtils::computeWeight( (recMu1+recMu2).mass(), iev, true );
706  if (debug_>0) {
707  std::cout << "Loop #" << loopCounter << "Event #" << iev << ": before correction Pt1 = "
708  << recMu1.Pt() << " Pt2 = " << recMu2.Pt() << std::endl;
709  }
710  // For successive iterations, correct the muons only if the previous iteration was a scale fit.
711  // --------------------------------------------------------------------------------------------
712  if ( loopCounter>0 ) {
716  }
717  }
718  if (debug_>0) {
719  std::cout << "Loop #" << loopCounter << "Event #" << iev << ": after correction Pt1 = "
720  << recMu1.Pt() << " Pt2 = " << recMu2.Pt() << std::endl;
721  }
722 
724 
725  //Fill histograms
726  //------------------
727  mapHisto_["hRecBestMu"]->Fill(recMu1);
728  mapHisto_["hRecBestMuVSEta"]->Fill(recMu1);
729  mapHisto_["hRecBestMu"]->Fill(recMu2);
730  mapHisto_["hRecBestMuVSEta"]->Fill(recMu2);
731  mapHisto_["hDeltaRecBestMu"]->Fill(recMu1, recMu2);
732  // Reconstructed resonance
733  mapHisto_["hRecBestRes"]->Fill(bestRecRes, weight);
734  mapHisto_["hRecBestResAllEvents"]->Fill(bestRecRes, 1.);
735  // Fill histogram of Res mass vs muon variables
736  mapHisto_["hRecBestResVSMu"]->Fill (recMu1, bestRecRes, -1);
737  mapHisto_["hRecBestResVSMu"]->Fill (recMu2, bestRecRes, +1);
738  // Fill histogram of Res mass vs Res variables
739  mapHisto_["hRecBestResVSRes"]->Fill (bestRecRes, bestRecRes, +1);
740 
741  std::vector<double> * parval;
742  std::vector<double> initpar;
743  // Store a pointer to the vector of parameters of the last iteration, or the initial
744  // parameters if this is the first iteration
745  if (loopCounter==0) {
746  initpar = MuScleFitUtils::parResol;
747  initpar.insert( initpar.end(), MuScleFitUtils::parScale.begin(), MuScleFitUtils::parScale.end() );
748  initpar.insert( initpar.end(), MuScleFitUtils::parCrossSection.begin(), MuScleFitUtils::parCrossSection.end() );
749  initpar.insert( initpar.end(), MuScleFitUtils::parBgr.begin(), MuScleFitUtils::parBgr.end() );
750  parval = &initpar;
751  } else {
752  parval = &(MuScleFitUtils::parvalue[loopCounter-1]);
753  }
754 
755  //Compute pt resolution w.r.t generated and simulated muons
756  //--------------------------------------------------------
757  if( !MuScleFitUtils::speedup ) {
758 
759  //first is always mu-, second is always mu+
762  }
765  }
766  if( compareToSimTracks_ ) {
767  //first is always mu-, second is always mu+
770  }
773  }
774  }
775  }
776 
777  // ATTENTION: this was done only when a matching was found. Moved it outside because, genInfo or not, we still want to see the resolution function
778  // Fill also the resolution histogramsm using the resolution functions:
779  // the parameters are those from the last iteration, as the muons up to this point have also the corrections from the same iteration.
780  // Need to use a different array (ForVec), containing functors able to operate on std::vector<double>
781  mapHisto_["hFunctionResolPt"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu1.Pt(), recMu1.Eta(), *parval ), -1 );
782  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu1.Pt(), recMu1.Eta(), *parval ), -1 );
783  mapHisto_["hFunctionResolPhi"]->Fill( recMu1, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu1.Pt(), recMu1.Eta(), *parval ), -1 );
784  mapHisto_["hFunctionResolPt"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPt(recMu2.Pt(), recMu2.Eta(), *parval ), +1 );
785  mapHisto_["hFunctionResolCotgTheta"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaCotgTh(recMu2.Pt(), recMu2.Eta(), *parval ), +1 );
786  mapHisto_["hFunctionResolPhi"]->Fill( recMu2, MuScleFitUtils::resolutionFunctionForVec->sigmaPhi(recMu2.Pt(), recMu2.Eta(), *parval ), +1 );
787 
788  // Compute likelihood histograms
789  // -----------------------------
790  if( debug_ > 0 ) std::cout << "mass = " << bestRecRes.mass() << std::endl;
791  if (weight!=0.) {
792  double massResol;
793  double prob;
794  double deltalike;
795  if (loopCounter==0) {
796  std::vector<double> initpar;
797  for (int i=0; i<(int)(MuScleFitUtils::parResol.size()); i++) {
798  initpar.push_back(MuScleFitUtils::parResol[i]);
799  }
800  for (int i=0; i<(int)(MuScleFitUtils::parScale.size()); i++) {
801  initpar.push_back(MuScleFitUtils::parScale[i]);
802  }
803 // for (int i=0; i<(int)(MuScleFitUtils::parCrossSection.size()); i++) {
804 // initpar.push_back(MuScleFitUtils::parCrossSection[i]);
805 // }
807 
808  for (int i=0; i<(int)(MuScleFitUtils::parBgr.size()); i++) {
809  initpar.push_back(MuScleFitUtils::parBgr[i]);
810  }
811  massResol = MuScleFitUtils::massResolution( recMu1, recMu2, initpar );
812  prob = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(), massResol, initpar, true );
813  } else {
816  prob = MuScleFitUtils::massProb( bestRecRes.mass(), bestRecRes.Eta(), bestRecRes.Rapidity(),
817  massResol, MuScleFitUtils::parvalue[loopCounter-1], true );
818  }
819  if( debug_ > 0 ) std::cout << "inside weight: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
820  if (prob>0) {
821  if( debug_ > 0 ) std::cout << "inside prob: mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
822 
823  deltalike = log(prob)*weight; // NB maximum likelihood --> deltalike is maximized
824  mapHisto_["hLikeVSMu"]->Fill(recMu1, deltalike);
825  mapHisto_["hLikeVSMu"]->Fill(recMu2, deltalike);
826  mapHisto_["hLikeVSMuMinus"]->Fill(recMu1, deltalike);
827  mapHisto_["hLikeVSMuPlus"]->Fill(recMu2, deltalike);
828 
829  double recoMass = (recMu1+recMu2).mass();
830  if( recoMass != 0 ) {
831  // IMPORTANT: massResol is not a relative resolution
832  mapHisto_["hResolMassVSMu"]->Fill(recMu1, massResol, -1);
833  mapHisto_["hResolMassVSMu"]->Fill(recMu2, massResol, +1);
834  mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu1, massResol/recoMass, -1);
835  mapHisto_["hFunctionResolMassVSMu"]->Fill(recMu2, massResol/recoMass, +1);
836  }
837 
839  mapHisto_["hdMdPt1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdpt1, -1);
840  mapHisto_["hdMdPt2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdpt2, +1);
841  mapHisto_["hdMdPhi1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdphi1, -1);
842  mapHisto_["hdMdPhi2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdphi2, +1);
843  mapHisto_["hdMdCotgTh1"]->Fill(recMu1, MuScleFitUtils::massResolComponents.dmdcotgth1, -1);
844  mapHisto_["hdMdCotgTh2"]->Fill(recMu2, MuScleFitUtils::massResolComponents.dmdcotgth2, +1);
845  }
846 
847  if( !MuScleFitUtils::speedup ) {
848  double genMass = (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second).mass();
849  // Fill the mass resolution (computed from MC), we use the covariance class to compute the variance
850  if( genMass != 0 ) {
852  mapHisto_["hGenResVSMu"]->Fill((MuScleFitUtils::genPair[iev].second), (MuScleFitUtils::genPair[iev].first + MuScleFitUtils::genPair[iev].second), +1);
853  double diffMass = (recoMass - genMass)/genMass;
854  // double diffMass = recoMass - genMass;
855  // Fill if for both muons
856  double pt1 = recMu1.pt();
857  double eta1 = recMu1.eta();
858  double pt2 = recMu2.pt();
859  double eta2 = recMu2.eta();
860  // This is to avoid nan
861  if( diffMass == diffMass ) {
862  // Mass relative difference vs Pt and Eta. To be used to extract the true mass resolution
863  mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt1, diffMass);
864  mapHisto_["hDeltaMassOverGenMassVsPt"]->Fill(pt2, diffMass);
865  mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta1, diffMass);
866  mapHisto_["hDeltaMassOverGenMassVsEta"]->Fill(eta2, diffMass);
867  // This is used for the covariance comparison
868  mapHisto_["hMassResolutionVsPtEta"]->Fill(pt1, eta1, diffMass, diffMass);
869  mapHisto_["hMassResolutionVsPtEta"]->Fill(pt2, eta2, diffMass, diffMass);
870  }
871  else {
872  std::cout << "Error, there is a nan: recoMass = " << recoMass << ", genMass = " << genMass << std::endl;
873  }
874  }
875  // Fill with mass resolution from resolution function
877  mapHisto_["hFunctionResolMass"]->Fill( recMu1, std::pow(massRes,2), -1 );
878  mapHisto_["hFunctionResolMass"]->Fill( recMu2, std::pow(massRes,2), +1 );
879  }
880 
881  mapHisto_["hMass_P"]->Fill(bestRecRes.mass(), prob);
882  if( debug_ > 0 ) std::cout << "mass = " << bestRecRes.mass() << ", prob = " << prob << std::endl;
883  mapHisto_["hMass_fine_P"]->Fill(bestRecRes.mass(), prob);
884 
885  mapHisto_["hMassProbVsRes"]->Fill(bestRecRes, bestRecRes, +1, prob);
886  mapHisto_["hMassProbVsMu"]->Fill(recMu1, bestRecRes, -1, prob);
887  mapHisto_["hMassProbVsMu"]->Fill(recMu2, bestRecRes, +1, prob);
888  mapHisto_["hMassProbVsRes_fine"]->Fill(bestRecRes, bestRecRes, +1, prob);
889  mapHisto_["hMassProbVsMu_fine"]->Fill(recMu1, bestRecRes, -1, prob);
890  mapHisto_["hMassProbVsMu_fine"]->Fill(recMu2, bestRecRes, +1, prob);
891  }
892  }
893  } // end if ResFound
894 
895  // Fill the pair
896  // -------------
897  if (loopCounter>0) {
898  if (debug_>0) std::cout << "[MuScleFit]: filling the pair" << std::endl;
899  MuScleFitUtils::SavedPair[iev] = std::make_pair( recMu1, recMu2 );
900  }
901 
902  iev++;
904 
905  // return kContinue;
906 }
907 
909  //first is always mu-, second is always mu+
910  double deltaR = sqrt(MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) * MuScleFitUtils::deltaPhi(recMu.Phi(),genMu.Phi()) +
911  ((recMu.Eta()-genMu.Eta()) * (recMu.Eta()-genMu.Eta())));
912  if(deltaR<0.01)
913  return true;
914  else if( debug_ > 0 ) {
915  std::cout<<"Reco muon "<<recMu<<" with eta "<<recMu.Eta()<<" and phi "<<recMu.Phi()<<std::endl
916  <<" DOES NOT MATCH with generated muon from resonance: "<<std::endl
917  <<genMu<<" with eta "<<genMu.Eta()<<" and phi "<<genMu.Phi()<<std::endl;
918  }
919  return false;
920 }
921 
923  const std::string & inputName, const int charge )
924 {
925  std::string name(inputName + "VSMu");
926  mapHisto_["hResolPt"+name]->Fill(recMu, (-genMu.Pt()+recMu.Pt())/genMu.Pt(), charge);
927  mapHisto_["hResolTheta"+name]->Fill(recMu, (-genMu.Theta()+recMu.Theta()), charge);
928  mapHisto_["hResolCotgTheta"+name]->Fill(recMu,(-cos(genMu.Theta())/sin(genMu.Theta())
929  +cos(recMu.Theta())/sin(recMu.Theta())), charge);
930  mapHisto_["hResolEta"+name]->Fill(recMu, (-genMu.Eta()+recMu.Eta()),charge);
931  mapHisto_["hResolPhi"+name]->Fill(recMu, MuScleFitUtils::deltaPhiNoFabs(recMu.Phi(), genMu.Phi()), charge);
932 
933  // Fill only if it was matched to a genMu and this muon is valid
934  if( (genMu.Pt() != 0) && (recMu.Pt() != 0) ) {
935  mapHisto_["hPtRecoVsPt"+inputName]->Fill(genMu.Pt(), recMu.Pt());
936  }
937 }
938 
940 {
941  if( MuScleFitUtils::SmearType>0 ) {
943  if (debug_>0) std::cout << "Muon #" << MuScleFitUtils::goodmuon
944  << ": after smearing Pt = " << mu.Pt() << std::endl;
945  }
946 }
947 
949 {
950  if( MuScleFitUtils::BiasType>0 ) {
951  mu = MuScleFitUtils::applyBias( mu, charge );
952  if (debug_>0) std::cout << "Muon #" << MuScleFitUtils::goodmuon
953  << ": after bias Pt = " << mu.Pt() << std::endl;
954  }
955 }
956 
957 // Simple method to check parameters consistency. It aborts the job if the parameters
958 // are not consistent.
960 
961  // Fits selection dimension check
963  std::cout << "[MuScleFit-Constructor]: wrong size of resolution fits selector = " << MuScleFitUtils::doResolFit.size() << std::endl;
964  std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
965  abort();
966  }
968  std::cout << "[MuScleFit-Constructor]: wrong size of scale fits selector = " << MuScleFitUtils::doScaleFit.size() << std::endl;
969  std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
970  abort();
971  }
973  std::cout << "[MuScleFit-Constructor]: wrong size of cross section fits selector = " << MuScleFitUtils::doCrossSectionFit.size() << std::endl;
974  std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
975  abort();
976  }
978  std::cout << "[MuScleFit-Constructor]: wrong size of background fits selector = " << MuScleFitUtils::doBackgroundFit.size() << std::endl;
979  std::cout << "it must have as many values as the number of loops, which is = " << maxLoopNumber << std::endl;
980  abort();
981  }
982 
983  // Bias parameters: dimension check
984  // --------------------------------
985  if ((MuScleFitUtils::BiasType==1 && MuScleFitUtils::parBias.size()!=2) || // linear in pt
986  (MuScleFitUtils::BiasType==2 && MuScleFitUtils::parBias.size()!=2) || // linear in |eta|
987  (MuScleFitUtils::BiasType==3 && MuScleFitUtils::parBias.size()!=4) || // sinusoidal in phi
988  (MuScleFitUtils::BiasType==4 && MuScleFitUtils::parBias.size()!=3) || // linear in pt and |eta|
989  (MuScleFitUtils::BiasType==5 && MuScleFitUtils::parBias.size()!=3) || // linear in pt and sinusoidal in phi
990  (MuScleFitUtils::BiasType==6 && MuScleFitUtils::parBias.size()!=3) || // linear in |eta| and sinusoidal in phi
991  (MuScleFitUtils::BiasType==7 && MuScleFitUtils::parBias.size()!=4) || // linear in pt and |eta| and sinusoidal in phi
992  (MuScleFitUtils::BiasType==8 && MuScleFitUtils::parBias.size()!=4) || // linear in pt and parabolic in |eta|
993  (MuScleFitUtils::BiasType==9 && MuScleFitUtils::parBias.size()!=2) || // exponential in pt
994  (MuScleFitUtils::BiasType==10 && MuScleFitUtils::parBias.size()!=3) || // parabolic in pt
995  (MuScleFitUtils::BiasType==11 && MuScleFitUtils::parBias.size()!=4) || // linear in pt and sin in phi with chg
996  (MuScleFitUtils::BiasType==12 && MuScleFitUtils::parBias.size()!=6) || // linear in pt and para in plus sin in phi with chg
997  (MuScleFitUtils::BiasType==13 && MuScleFitUtils::parBias.size()!=8) || // linear in pt and para in plus sin in phi with chg
998  MuScleFitUtils::BiasType<0 || MuScleFitUtils::BiasType>13) {
999  std::cout << "[MuScleFit-Constructor]: Wrong bias type or number of parameters: aborting!" << std::endl;
1000  abort();
1001  }
1002  // Smear parameters: dimension check
1003  // ---------------------------------
1011  MuScleFitUtils::SmearType<0 || MuScleFitUtils::SmearType>7) {
1012  std::cout << "[MuScleFit-Constructor]: Wrong smear type or number of parameters: aborting!" << std::endl;
1013  abort();
1014  }
1015  // Protect against bad size of parameters
1016  // --------------------------------------
1019  std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Resol: aborting!" << std::endl;
1020  abort();
1021  }
1024  std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Scale: aborting!" << std::endl;
1025  abort();
1026  }
1029  std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1030  abort();
1031  }
1034  std::cout << "[MuScleFit-Constructor]: Mismatch in number of parameters for Bgr: aborting!" << std::endl;
1035  abort();
1036  }
1037 
1038  // Protect against an incorrect number of resonances
1039  // -------------------------------------------------
1040  if (MuScleFitUtils::resfind.size()!=6) {
1041  std::cout << "[MuScleFit-Constructor]: resfind must have 6 elements (1 Z, 3 Y, 2 Psi): aborting!" << std::endl;
1042  abort();
1043  }
1044 }
1045 
1047 
1048  reco::TrackRef iTrack = aMuon->innerTrack();
1049  const reco::HitPattern& p = iTrack->hitPattern();
1050 
1051  reco::TrackRef gTrack = aMuon->globalTrack();
1052  const reco::HitPattern& q = gTrack->hitPattern();
1053 
1054  return (//isMuonInAccept(aMuon) &&// no acceptance cuts!
1055  iTrack->found() > 11 &&
1056  gTrack->chi2()/gTrack->ndof() < 20.0 &&
1057  q.numberOfValidMuonHits() > 0 &&
1058  iTrack->chi2()/iTrack->ndof() < 4.0 &&
1059  aMuon->muonID("TrackerMuonArbitrated") &&
1060  aMuon->muonID("TMLastStationAngTight") &&
1061  p.pixelLayersWithMeasurement() > 1 &&
1062  fabs(iTrack->dxy()) < 3.0 && //should be done w.r.t. PV!
1063  fabs(iTrack->dz()) < 15.0 );//should be done w.r.t. PV!
1064 }
1065 
1066 
1068 
1069  reco::TrackRef iTrack = aMuon->innerTrack();
1070  const reco::HitPattern& p = iTrack->hitPattern();
1071 
1072  return (//isMuonInAccept(aMuon) // no acceptance cuts!
1073  iTrack->found() > 11 &&
1074  iTrack->chi2()/iTrack->ndof() < 4.0 &&
1075  aMuon->muonID("TrackerMuonArbitrated") &&
1076  aMuon->muonID("TMLastStationAngTight") &&
1077  p.pixelLayersWithMeasurement() > 1 &&
1078  fabs(iTrack->dxy()) < 3.0 && //should be done w.r.t. PV!
1079  fabs(iTrack->dz()) < 15.0 );//should be done w.r.t. PV!
1080 
1081 }
static double deltaPhiNoFabs(const double &phi1, const double &phi2)
Without fabs at the end, used to have a symmetric distribution for the resolution fits and variance c...
static std::vector< std::pair< lorentzVector, lorentzVector > > simPair
const double Pi
static std::vector< int > doScaleFit
static std::vector< int > doResolFit
T getParameter(std::string const &) const
scaleFunctionBase< double * > * scaleFunctionService(const int identifier)
Service to build the scale functor corresponding to the passed identifier.
Definition: Functions.cc:3
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
static std::vector< double > parBias
void selectMuons(const edm::Event &event)
Definition: MuScleFit.cc:590
virtual edm::EDLooper::Status endOfLoop(const edm::EventSetup &eventSetup, unsigned int iLoop)
Definition: MuScleFit.cc:445
static std::vector< int > parCrossSectionOrder
std::vector< GenMuonPair > genMuonPairs_
Stores the genMuon pairs and the motherId prior to the creation of the internal tree.
Definition: MuScleFitBase.h:82
reco::TrackRef innerTrack() const
reference to Track reconstructed in the tracker only (reimplemented from reco::Muon) ...
Definition: Muon.h:74
static smearFunctionBase * smearFunction
static std::vector< int > parScaleOrder
edm::InputTag theMuonLabel_
Definition: MuScleFitBase.h:46
static std::vector< int > doCrossSectionFit
static void minimizeLikelihood()
static double maxMuonEtaSecondRange_
static std::vector< int > parBgrOrder
bool selGlobalMuon(const pat::Muon *aMuon)
Function for onia selections.
Definition: MuScleFit.cc:1046
static unsigned int loopCounter
unsigned int loopCounter
Definition: MuScleFit.h:140
bool muonID(const std::string &name) const
Definition: Muon.cc:280
int numberOfValidMuonHits() const
Definition: HitPattern.cc:346
int totalEvents_
Definition: MuScleFit.h:150
static bool startWithSimplex_
std::string triggerResultsProcess_
Definition: MuScleFit.h:165
static std::vector< int > doBackgroundFit
void fillRec(std::vector< reco::LeafCandidate > &muons)
static std::vector< double > parResol
static bool debugMassResol_
bool saveAllToTree_
Definition: MuScleFit.h:168
virtual void duringFastLoop()
Definition: MuScleFit.cc:682
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void checkParameters()
Definition: MuScleFit.cc:959
static BackgroundHandler * backgroundHandler
static double x[7][10000]
static int debug
static double ResMass[6]
static double massWindowHalfWidth[3][6]
static bool speedup
Run const & getRun() const
Definition: Event.cc:44
std::map< std::string, Histograms * > mapHisto_
The map of histograms.
Definition: MuScleFitBase.h:77
static bool ResFound
static scaleFunctionBase< std::vector< double > > * biasFunction
bool selTrackerMuon(const pat::Muon *aMuon)
Definition: MuScleFit.cc:1067
static std::vector< int > parBgrFix
static bool minimumShapePlots_
int pixelLayersWithMeasurement() const
Definition: HitPattern.cc:888
edm::InputTag simTracksCollection_
Definition: MuScleFit.h:153
static int MuonTypeForCheckMassWindow
static double minMuonEtaFirstRange_
static double massProb(const double &mass, const double &rapidity, const int ires, const double &massResol)
void fillComparisonHistograms(const reco::Particle::LorentzVector &genMu, const reco::Particle::LorentzVector &recoMu, const std::string &inputName, const int charge)
Fill the reco vs gen and reco vs sim comparison histograms.
Definition: MuScleFit.cc:922
static struct MuScleFitUtils::massResolComponentsStruct massResolComponents
double charge(const std::vector< uint8_t > &Ampls)
static int BiasType
static std::vector< int > parScaleFix
static bool computeMinosErrors_
reco::Particle::LorentzVector lorentzVector
Definition: GenMuonPair.h:8
static scaleFunctionBase< std::vector< double > > * scaleFunctionForVec
static std::vector< std::vector< double > > parvalue
unsigned int triggerIndex(const std::string &triggerName) const
slot position of trigger path in trigger table (0 to size-1)
U second(std::pair< T, U > const &p)
static double maxMuonPt_
static int ScaleFitType
static std::vector< std::pair< lorentzVector, lorentzVector > > genPair
std::string outputRootTreeFileName_
Definition: MuScleFit.h:160
std::string theGenInfoRootFileName_
Definition: MuScleFitBase.h:48
resolutionFunctionBase< std::vector< double > > * resolutionFunctionVecService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier when receiving a std::...
Definition: Functions.cc:116
virtual void startingNewLoop(unsigned int iLoop)
Definition: MuScleFit.cc:418
static double massResolution(const lorentzVector &mu1, const lorentzVector &mu2)
reco::Particle::LorentzVector recMu2
Definition: MuScleFit.h:148
bool PATmuons_
Definition: MuScleFit.h:154
virtual void endOfJob()
Definition: MuScleFit.cc:412
reco::Particle::LorentzVector recMu1
Definition: MuScleFit.h:148
static lorentzVector applyBias(const lorentzVector &muon, const int charge)
T sqrt(T t)
Definition: SSEVec.h:28
static std::vector< double > parBgr
unsigned int maxLoopNumber
Definition: MuScleFit.h:139
tuple norm
Definition: lumiNorm.py:78
static int SmearType
std::string genParticlesName_
Definition: MuScleFit.h:155
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
RunNumber_t run() const
Definition: Event.h:66
void clearHistoMap()
Clean the histograms map.
int j
Definition: DBlmapReader.cc:9
void readTree(const int maxEvents, const TString &fileName, MuonPairVector *savedPair, const int muonType, MuonPairVector *genPair=0)
bool compareToSimTracks_
Definition: MuScleFit.h:152
static double computeWeight(const double &mass, const int iev, const bool doUseBkgrWindow=false)
static std::vector< std::pair< lorentzVector, lorentzVector > > SavedPair
MuScleFit(const edm::ParameterSet &pset)
Definition: MuScleFit.cc:148
static lorentzVector applyScale(const lorentzVector &muon, const std::vector< double > &parval, const int charge)
static std::vector< int > parResolFix
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
bool first
Definition: L1TdeRCT.cc:79
static std::vector< double > parSmear
void beginOfJobInConstructor()
Definition: MuScleFit.cc:383
smearFunctionBase * smearFunctionService(const int identifier)
Service to build the smearing functor corresponding to the passed identifier.
Definition: Functions.cc:73
reco::TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector (reimplemented from reco::Muon) ...
Definition: Muon.h:82
static std::vector< int > parResolOrder
static bool sherpa_
std::string triggerPath_
Definition: MuScleFit.h:166
std::auto_ptr< MuScleFitMuonSelector > muonSelector_
Definition: MuScleFit.h:170
bool fastLoop
Definition: MuScleFit.h:142
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
void writeTree(const TString &fileName, const std::vector< MuonPair > *savedPair, const int muonType=0, const std::vector< GenMuonPair > *genPair=0, const bool saveAll=false)
scaleFunctionBase< std::vector< double > > * scaleFunctionVecService(const int identifier)
Service to build the scale functor corresponding to the passed identifier when receiving a std::vecto...
Definition: Functions.cc:38
void writeHistoMap(const unsigned int iLoop)
Save the histograms map to file.
static std::vector< double > parScale
static double oldNormalization_
void applySmearing(reco::Particle::LorentzVector &mu)
Apply the smearing if needed using the function in MuScleFitUtils.
Definition: MuScleFit.cc:939
Log< T >::type log(const T &t)
Definition: Log.h:22
virtual ~MuScleFit()
Definition: MuScleFit.cc:357
std::string theRootFileName_
Definition: MuScleFitBase.h:47
bool init(const edm::Run &iRun, const edm::EventSetup &iSetup, const std::string &processName, bool &changed)
std::string triggerResultsLabel_
Definition: MuScleFit.h:164
static bool rapidityBinsForZ_
MuScleFitPlotter * plotter
Definition: MuScleFit.h:144
void fillHistoMap(TFile *outputFile, unsigned int iLoop)
Create the histograms map.
Definition: MuScleFitBase.cc:7
std::vector< TFile * > theFiles_
The files were the histograms are saved.
Definition: MuScleFitBase.h:74
double minResMass_hwindow[6]
Definition: MuScleFit.h:134
static int ResolFitType
static int goodmuon
static double minMuonEtaSecondRange_
void fillGen(const std::vector< std::pair< reco::Particle::LorentzVector, reco::Particle::LorentzVector > > &genPairs)
static int iev_
static std::vector< double > parCrossSection
tuple muons
Definition: patZpeak.py:38
static lorentzVector applySmearing(const lorentzVector &muon)
static bool normalizeLikelihoodByEventNumber_
std::string inputRootTreeFileName_
Definition: MuScleFit.h:158
bool checkDeltaR(reco::Particle::LorentzVector &genMu, reco::Particle::LorentzVector &recMu)
Check if two lorentzVector are near in deltaR.
Definition: MuScleFit.cc:908
resolutionFunctionBase< double * > * resolutionFunctionService(const int identifier)
Service to build the resolution functor corresponding to the passed identifier.
Definition: Functions.cc:88
int maxEventsFromRootTree_
Definition: MuScleFit.h:162
double maxResMass_hwindow[6]
Definition: MuScleFit.h:135
static int MuonType
std::vector< MuonPair > muonPairs_
Used to store the muon pairs plus run and event number prior to the creation of the internal tree...
Definition: MuScleFitBase.h:80
static double deltaPhi(const double &phi1, const double &phi2)
tuple cout
Definition: gather_cfg.py:41
static double minMuonPt_
static scaleFunctionBase< double * > * scaleFunction
static bool useProbsFile_
static std::vector< int > resfind
void readProbabilityDistributionsFromFile()
Read probability distributions from a local root file.
static int counter_resprob
static resolutionFunctionBase< std::vector< double > > * resolutionFunctionForVec
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:26
static int FitStrategy
Analysis-level muon class.
Definition: Muon.h:51
static double maxMuonEtaFirstRange_
tuple size
Write out results.
static std::pair< lorentzVector, lorentzVector > findBestRecoRes(const std::vector< reco::LeafCandidate > &muons)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void applyBias(reco::Particle::LorentzVector &mu, const int charge)
Apply the bias if needed using the function in MuScleFitUtils.
Definition: MuScleFit.cc:948
bool negateTrigger_
Definition: MuScleFit.h:167
static CrossSectionHandler * crossSectionHandler
static std::vector< int > parCrossSectionFix
static resolutionFunctionBase< double * > * resolutionFunction
virtual edm::EDLooper::Status duringLoop(const edm::Event &event, const edm::EventSetup &eventSetup)
Definition: MuScleFit.cc:528
void addParameters(std::vector< double > &initpar)
Inputs the vars in a vector.
virtual void endOfFastLoop(const unsigned int iLoop)
Definition: MuScleFit.cc:492