Core/src/FWMagField.cc

0001 #include <iostream>
0002
0003 #include "TH1F.h"
0004 #include "Fireworks/Core/interface/FWMagField.h"
0005 #include "Fireworks/Core/interface/fwLog.h"
0006 #include "DataFormats/FWLite/interface/Handle.h"
0007 #include "DataFormats/Common/interface/ConditionsInEdm.h"
0008 #include "DataFormats/Scalers/interface/DcsStatus.h"
0009 #include "DataFormats/FWLite/interface/Event.h"
0010
0011 FWMagField::FWMagField()
0012     : TEveMagField(),
0013
0014       m_source(kNone),
0015       m_userField(-1),
0016       m_eventField(-1),
0017
0018       m_reverse(true),
0019       m_simpleModel(false),
0020
0021       m_guessValHist(nullptr),
0022       m_numberOfFieldIsOnEstimates(0),
0023       m_numberOfFieldEstimates(0),
0024       m_updateFieldEstimate(true),
0025       m_guessedField(0) {
0026   m_guessValHist = new TH1F("FieldEstimations", "Field estimations from tracks and muons", 200, -4.5, 4.5);
0027   m_guessValHist->SetDirectory(nullptr);
0028 }
0029
0030 FWMagField::~FWMagField() { delete m_guessValHist; }
0031
0032 //______________________________________________________________________________
0033
0034 TEveVector FWMagField::GetField(Float_t x, Float_t y, Float_t z) const {
0035   // Virtual method of TEveMagField class.
0036
0037   Float_t R = sqrt(x * x + y * y);
0038   Float_t field = m_reverse ? -GetFieldMag() : GetFieldMag();
0039
0040   //barrel
0041   if (TMath::Abs(z) < 724) {
0042     //inside solenoid
0043     if (R < 300)
0044       return TEveVector(0, 0, field);
0045     // outside solinoid
0046     if (m_simpleModel || (R > 461.0 && R < 490.5) || (R > 534.5 && R < 597.5) || (R > 637.0 && R < 700.0))
0047       return TEveVector(0, 0, -field / 3.8 * 1.2);
0048
0049   } else {
0050     // endcaps
0051     if (m_simpleModel) {
0052       if (R < 50)
0053         return TEveVector(0, 0, field);
0054       if (z > 0)
0055         return TEveVector(x / R * field / 3.8 * 2.0, y / R * field / 3.8 * 2.0, 0);
0056       else
0057         return TEveVector(-x / R * field / 3.8 * 2.0, -y / R * field / 3.8 * 2.0, 0);
0058     }
0059     // proper model
0060     if (((TMath::Abs(z) > 724) && (TMath::Abs(z) < 786)) || ((TMath::Abs(z) > 850) && (TMath::Abs(z) < 910)) ||
0061         ((TMath::Abs(z) > 975) && (TMath::Abs(z) < 1003))) {
0062       if (z > 0)
0063         return TEveVector(x / R * field / 3.8 * 2.0, y / R * field / 3.8 * 2.0, 0);
0064       else
0065         return TEveVector(-x / R * field / 3.8 * 2.0, -y / R * field / 3.8 * 2.0, 0);
0066     }
0067   }
0068   return TEveVector(0, 0, 0);
0069 }
0070
0071 //______________________________________________________________________________
0072
0073 Float_t FWMagField::GetFieldMag() const {
0074   float res;
0075   switch (m_source) {
0076     case kEvent: {
0077       res = m_eventField;
0078       break;
0079     }
0080     case kUser: {
0081       res = m_userField;
0082       break;
0083     }
0084     default: {
0085       if (m_updateFieldEstimate) {
0086         if (m_guessValHist->GetEntries() > 2 && m_guessValHist->GetRMS() < 0.5) {
0087           m_guessedField = m_guessValHist->GetMean();
0088
0089           // std::cout << "FWMagField::GetFieldMag(), get average "
0090           //  << m_guessValHist->GetMean() << " guessed value: RMS= "<< m_guessValHist->GetRMS()
0091           //  <<" samples "<< m_guessValHist->GetEntries() << std::endl;
0092
0093         } else if (m_numberOfFieldIsOnEstimates > m_numberOfFieldEstimates / 2 || m_numberOfFieldEstimates == 0) {
0094           m_guessedField = 3.8;
0095           // fwLog(fwlog::kDebug) << "FWMagField::GetFieldMag() get default field, number estimates "
0096           //  << m_numberOfFieldEstimates << " number fields is on  m_numberOfFieldIsOnEstimates" <<std::endl;
0097         } else {
0098           m_guessedField = 0;
0099           //  fwLog(fwlog::kDebug) << "Update field estimate, guess field is OFF." <<std::endl;
0100         }
0101         m_updateFieldEstimate = false;
0102       }
0103
0104       res = m_guessedField;
0105     }
0106   }
0107
0108   return res;
0109 }
0110
0111 Float_t FWMagField::GetMaxFieldMag() const {
0112   // Runge-Kutta stepper does not like this to be zero.
0113   // Will be fixed in root.
0114   // The return value should definitley not be negative -- so Abs
0115   // should stay!
0116
0117   return TMath::Max(TMath::Abs(GetFieldMag()), 0.01f);
0118 }
0119
0120 //______________________________________________________________________________
0121
0122 void FWMagField::guessFieldIsOn(bool isOn) const {
0123   if (isOn)
0124     ++m_numberOfFieldIsOnEstimates;
0125   ++m_numberOfFieldEstimates;
0126   m_updateFieldEstimate = true;
0127 }
0128
0129 void FWMagField::guessField(float val) const {
0130   // fwLog(filedDebug) <<  "FWMagField::guessField "<< val << std::endl;
0131   m_guessValHist->Fill(val);
0132   m_updateFieldEstimate = true;
0133 }
0134
0135 void FWMagField::resetFieldEstimate() const {
0136   m_guessValHist->Reset();
0137   m_guessValHist->SetAxisRange(-4, 4);
0138   m_numberOfFieldIsOnEstimates = 0;
0139   m_numberOfFieldEstimates = 0;
0140   m_updateFieldEstimate = true;
0141 }
0142
0143 //______________________________________________________________________________
0144 void FWMagField::checkFieldInfo(const edm::EventBase* event) {
0145   const static float currentToField = 3.8 / 18160;
0146   bool available = false;
0147   try {
0148     edm::InputTag conditionsTag("conditionsInEdm");
0149     edm::Handle<edm::ConditionsInRunBlock> runCond;
0150     // FIXME: ugly hack to avoid exposing an fwlite::Event. Need to ask
0151     //        Chris / Ianna how to get mag field from an EventBase.
0152     const fwlite::Event* fwEvent = dynamic_cast<const fwlite::Event*>(event);
0153     if (!fwEvent)
0154       return;
0155
0156     m_source = kNone;
0157     fwEvent->getRun().getByLabel(conditionsTag, runCond);
0158
0159     if (runCond.isValid()) {
0160       available = true;
0161       m_eventField = currentToField * runCond->BAvgCurrent;
0162       m_source = kEvent;
0163       fwLog(fwlog::kDebug) << "Magnetic field info found in ConditionsInEdm branch : " << m_eventField << std::endl;
0164     } else {
0165       edm::InputTag dcsTag("scalersRawToDigi");
0166       edm::Handle<std::vector<DcsStatus> > dcsStatus;
0167       event->getByLabel(dcsTag, dcsStatus);
0168
0169       if (dcsStatus.isValid() && !dcsStatus->empty()) {
0170         float sum = 0;
0171         for (std::vector<DcsStatus>::const_iterator i = dcsStatus->begin(); i < dcsStatus->end(); ++i)
0172           sum += (*i).magnetCurrent();
0173
0174         available = true;
0175         m_eventField = currentToField * sum / dcsStatus->size();
0176         m_source = kEvent;
0177         fwLog(fwlog::kDebug) << "Magnetic field info found in DcsStatus branch: " << m_eventField << std::endl;
0178       }
0179     }
0180   } catch (cms::Exception&) {
0181     fwLog(fwlog::kDebug) << "Cought exception in FWMagField::checkFieldInfo\n";
0182   }
0183
0184   if (!available) {
0185     fwLog(fwlog::kDebug) << "No magnetic field info available in Event\n";
0186   }
0187 }
0188
0189 //______________________________________________________________________________
0190 void FWMagField::setFFFieldMag(float mag) {
0191   // AMT this is a workaround for seting FF in FFLooper
0192   // Correct imeplementation is having a base class of  FWMagField amd do implementation for FF and FWLite version
0193
0194   m_source = kEvent;
0195   m_eventField = mag;
0196 }