src/smear/Detector.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::Detector.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-19
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/smear/Detector.h"
0011
0012 #include <algorithm>
0013 #include <functional>
0014 #include <list>
0015 #include <memory>
0016 #include <vector>
0017
0018 #include "eicsmear/erhic/EventDis.h"
0019 #include "eicsmear/smear/EventSmear.h"
0020 #include "eicsmear/erhic/Kinematics.h"
0021 #include "eicsmear/smear/ParticleMCS.h"
0022 #include "eicsmear/smear/Smearer.h"
0023 #include "eicsmear/erhic/VirtualParticle.h"
0024
0025 using std::cout;
0026 using std::cerr;
0027 using std::endl;
0028
0029 namespace Smear {
0030
0031 Detector::Detector()
0032 : useNM(false)
0033 , useJB(false)
0034 , useDA(false) {
0035 }
0036
0037 Detector::Detector(const Detector& other)
0038 : TObject(other) {
0039   useNM = other.useNM;
0040   useJB = other.useJB;
0041   useDA = other.useDA;
0042   Devices = other.CopyDevices();
0043   LegacyMode = other.GetLegacyMode();
0044 }
0045
0046 Detector& Detector::operator=(const Detector& that) {
0047   if (this != &that) {
0048     useNM = that.useNM;
0049     useJB = that.useJB;
0050     useDA = that.useDA;
0051     Devices = that.CopyDevices();
0052     LegacyMode = that.GetLegacyMode();
0053   }  // if
0054   return *this;
0055 }
0056
0057 Detector::~Detector() {
0058   DeleteAllDevices();
0059 }
0060
0061 void Detector::DeleteAllDevices() {
0062   for (unsigned i(0); i < GetNDevices(); i++) {
0063     delete Devices.at(i);
0064     Devices.at(i) = NULL;
0065   }  // for
0066   Devices.clear();
0067 }
0068
0069 void Detector::AddDevice(Smearer& dev) {
0070   Devices.push_back(dev.Clone());
0071 }
0072
0073 void Detector::SetEventKinematicsCalculator(TString s) {
0074   s.ToLower();
0075   useNM = s.Contains("nm") || s.Contains("null");
0076   useJB = s.Contains("jb") || s.Contains("jacquet");
0077   useDA = s.Contains("da") || s.Contains("double");
0078 }
0079
0080 Smearer* Detector::GetDevice(int n) {
0081   Smearer* smearer(NULL);
0082   if (unsigned(n) < Devices.size()) {
0083     smearer = Devices.at(n);
0084   }  // if
0085   return smearer;
0086 }
0087
0088 void Detector::FillEventKinematics(Event* eventS) {
0089   if (!(useNM || useJB || useDA)) {
0090     return;
0091   }  // if
0092      // Need a bit of jiggery-pokery here, as the incident beam info isn't
0093      // associated with the smeared event.
0094      // So, get the beam info from the MC event, but replace the scattered
0095      // electron with the smeared version.
0096      // Then we can use the standard JB/DA algorithms on the smeared event.
0097   const ParticleMCS* scattered = eventS->ScatteredLepton();
0098   typedef std::unique_ptr<erhic::DisKinematics> KinPtr;
0099   if (useNM && scattered) {
0100     KinPtr kin(erhic::LeptonKinematicsComputer(*eventS).Calculate());
0101     if (kin.get()) {
0102       eventS->SetLeptonKinematics(*kin);
0103     }  // if
0104   } else {
0105     eventS->SetLeptonKinematics( erhic::DisKinematics(-1., -1., -1., -1., -1.));
0106   }  // if
0107   if (useJB) {
0108     KinPtr kin(erhic::JacquetBlondelComputer(*eventS).Calculate());
0109     if (kin.get()) {
0110       eventS->SetJacquetBlondelKinematics(*kin);
0111     }  // if
0112   }  // if
0113   if (useDA && scattered) {
0114     KinPtr kin(erhic::DoubleAngleComputer(*eventS).Calculate());
0115     if (kin.get()) {
0116       eventS->SetDoubleAngleKinematics(*kin);
0117     }  // if
0118   }  // if
0119 }
0120
0121 std::list<Smearer*> Detector::Accept(const erhic::VirtualParticle& p) const {
0122   std::list<Smearer*> devices;
0123   // Only accept final-state particles, so skip the check against each
0124   // devices for non-final-state particles.
0125   if (p.GetStatus() == 1) {
0126     std::vector<Smearer*>::const_iterator iter;
0127     for (iter = Devices.begin(); iter != Devices.end(); ++iter) {
0128       // Store each device that accepts the particle.
0129       if ((*iter)->Accept.Is(p)) {
0130         devices.push_back(*iter);
0131       }  // if
0132     }  // for
0133   }  // if
0134   return devices;
0135 }
0136
0137 ParticleMCS* Detector::Smear(const erhic::VirtualParticle& prt) const {
0138   // Does the particle fall in the acceptance of any device?
0139   // If so, we smear it, if not, we skip it (store a NULL pointer).
0140   std::list<Smearer*> devices = Accept(prt);
0141   ParticleMCS* prtOut(NULL);
0142   if (!devices.empty()) {
0143     // It passes through at least one device, so smear it.
0144     // Devices in which it doesn't pass won't smear it.
0145     prtOut = new ParticleMCS();
0146     prtOut->SetSmeared();
0147     std::list<Smearer*>::iterator iter;
0148     for (iter = devices.begin(); iter != devices.end(); ++iter) {
0149       (*iter)->Smear(prt, *prtOut);
0150     }  // for
0151
0152     if (LegacyMode){
0153       // Compute derived momentum components.
0154       prtOut->SetPx( prtOut->GetP() * sin(prtOut->GetTheta()) * cos(prtOut->GetPhi()));
0155       prtOut->SetPy( prtOut->GetP() * sin(prtOut->GetTheta()) * sin(prtOut->GetPhi()));
0156       prtOut->SetPt( sqrt(pow(prtOut->GetPx(), 2.) + pow(prtOut->GetPy(), 2.)));
0157       prtOut->SetPz( prtOut->GetP() * cos(prtOut->GetTheta()));
0158
0159     } else { // not in LegacyMode
0160       // Sanity check and computation of derived momentum components.
0161       // ---------------------------------------------------------------------
0162       // - this cannot happen...
0163       if ( !prtOut->IsSmeared() ) throw std::runtime_error ("particle seems to be not smeared?!");
0164
0165       // Is momentum smeared at all?
0166       int MomComponentsChanged = prtOut->IsPSmeared() + prtOut->IsPtSmeared() + prtOut->IsPxSmeared() + prtOut->IsPySmeared() + prtOut->IsPzSmeared();
0167       int AngComponentsChanged = prtOut->IsPhiSmeared() + prtOut->IsThetaSmeared();
0168
0169       if ( MomComponentsChanged==0 ){
0170     // Momentum is untouched (pure calo measurement)
0171     // all we have to do is ensure phi and theta are explicitly smeared
0172     if ( !prtOut->IsPhiSmeared() ) {
0173       cerr << "Phi always needs to be smeared (at least with sigma=0)" << endl;
0174       cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0175       throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0176     }
0177     if ( !prtOut->IsThetaSmeared() ){
0178       cerr << "Theta always needs to be smeared (at least with sigma=0)" << endl;
0179       cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0180       throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0181     }
0182       } else if ( MomComponentsChanged + AngComponentsChanged != 3){
0183     // - Momentum is exactly defined by three independent quantities, including theta and phi
0184     cerr << "Expected 0 (excluding phi, theta) or exactly 3 (excluding phi, theta) smeared momentum quantities." << endl;
0185     cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0186     cerr << "MomComponentsChanged = " << MomComponentsChanged << endl;
0187     cerr << prt.GetEta() << endl;
0188     cerr << prtOut->IsPSmeared() << endl;
0189     cerr << prtOut->IsPtSmeared() << endl;
0190     cerr << prtOut->IsPxSmeared() << endl;
0191     cerr << prtOut->IsPySmeared() << endl;
0192     cerr << prtOut->IsPzSmeared() << endl;
0193     cerr << "AngComponentsChanged = " << AngComponentsChanged << endl;
0194     cerr << " pid : " << prt.Id() << endl;
0195     throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0196       } else {
0197     // We now have exactly three out of P, px, py, pz, pt, phi, theta. Compute the rest.
0198     // That's 35 total cases, but luckily many of them are not independent, or nonsensical in a detector.
0199     // NOTE: We need p^2 >= pT^2, pz^2.
0200     // Smearing (P and Pz) or (P and pT) is obscure enough to where we just make the ad-hoc decision
0201     // to adjust P in such a case.
0202
0203     // - px, py, pt, phi are intimately related. Let's condense.
0204     if ( prtOut->IsPxSmeared() ^ prtOut->IsPySmeared() ) {// "^" = XOR
0205       cerr << "Smearing only one out of px, py is not supported. Please contact the authors." << endl;
0206       cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0207       throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0208     } // illegal px, py --> removes 2 * ( 5 choose 2 ) = 20 combinations
0209
0210     if ( prtOut->IsPxSmeared() && prtOut->IsPySmeared() ) {
0211       if ( prtOut->IsPhiSmeared() ) {
0212         cerr << "Smearing px, py, and phi is inconsistent" << endl;
0213         cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0214         throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0215       }// 1, running 21
0216       if ( prtOut->IsPtSmeared() ) {
0217         cerr << "Smearing px, py, and pt is inconsistent" << endl;
0218         cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0219         throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0220       }  // 1, running 22
0221
0222       // Can set phi and pt now
0223       prtOut->SetPhi( std::atan2( prtOut->GetPy(),prtOut->GetPx() ) );
0224       prtOut->SetPt( std::sqrt( std::pow( prtOut->GetPx(),2) + std::pow( prtOut->GetPy(),2) ) );
0225
0226       // legal options remaining: pz, P, theta
0227       if ( prtOut->IsPzSmeared() ){ // pz is smeared
0228         prtOut->SetTheta( std::atan2(prtOut->GetPt() ,prtOut->GetPz() ) );
0229         prtOut->SetP( std::sqrt( std::pow( prtOut->GetPt(),2) + std::pow( prtOut->GetPz(),2) ) );
0230       } // 1, running 23
0231
0232       if ( prtOut->IsPSmeared() ){ // p is smeared
0233         if ( prtOut->GetP() < prtOut->GetPt() ) prtOut->SetP(prtOut->GetPt(), false);
0234         prtOut->SetTheta( std::atan2(prtOut->GetPt() ,prtOut->GetPz() ) );
0235         prtOut->SetPz( std::sqrt( std::pow(prtOut->GetP(), 2.) - std::pow(prtOut->GetPt(), 2.)) );
0236       } // 1, running 24
0237
0238       if ( prtOut->IsThetaSmeared() ){ // theta is smeared
0239         // Note: Here (as in other cases), we rely on the device to deliver physically sound
0240         // numbers (see and adapt HandleBogusValues). But for now we explicitly fault on division by zero
0241         assert ( fabs(std::tan(prtOut->GetTheta())) > 1e-8 );
0242         prtOut->SetPz( prtOut->GetPt() / std::tan(prtOut->GetTheta()) );
0243         prtOut->SetP( std::sqrt( std::pow( prtOut->GetPt(),2) + std::pow( prtOut->GetPz(),2) ) );
0244       } // 1, running 25
0245
0246     } // know both px and py --> knocked off 5 cases, 25 running
0247
0248     // We're left with px and py unsmeared, and choose 3 out of P, pt, pz, phi, theta = 10 combinations.
0249     // But we can reject the case of unsmeared phi since without px, py we can't reconstruct azimuth
0250     // No need to protect with another if() because we touched Phi in the previous clause
0251     if ( !prtOut->IsPhiSmeared() ) {
0252       cerr << "Momentum components are smeared, but neither phi nor px and py are." << endl;
0253       cerr << "For legacy smear scripts, use det.SetLegacyMode ( true );" << endl;
0254       throw std::runtime_error ("Failed consistency check in Detector::Smear()");
0255     } // 4 cases (4 choose 3), 29 running
0256
0257     // Leaves (choose 2 out of P, pz, pt, theta) = 6 combinations, math checks out.
0258     // Using a complementary if clause instead of another else for readability
0259     if ( !prtOut->IsPxSmeared() && !prtOut->IsPySmeared() ) {
0260       // NOTE: Currently, this is the only loop that should matter
0261       // since we don't allow px and py smearing in the formulas,
0262       // but there's no reason not to allow it in the future
0263       if(prtOut->IsPSmeared() && prtOut->IsThetaSmeared() ) /*(1100)*/ { // P, theta
0264
0265         prtOut->SetPt ( prtOut->GetP() * std::sin(prtOut->GetTheta()));
0266         prtOut->SetPz ( prtOut->GetP() * std::cos(prtOut->GetTheta()));
0267
0268       } else if( prtOut->IsPSmeared() && prtOut->IsPtSmeared() ) /*(1010)*/ { // P, pt
0269         if ( prtOut->GetP() < prtOut->GetPt() ) prtOut->SetP(prtOut->GetPt(), false);
0270         prtOut->SetPz( std::sqrt(std::pow(prtOut->GetP(), 2) - std::pow(prtOut->GetPt(), 2)));
0271         prtOut->SetTheta ( std::atan2(prtOut->GetPt(),prtOut->GetPz()));
0272
0273       } else if(prtOut->IsPzSmeared() && prtOut->IsPtSmeared() ) /*(0011)*/ { // pt, pz
0274
0275         prtOut->SetP( std::sqrt(std::pow(prtOut->GetPt(), 2) + std::pow(prtOut->GetPz(), 2)));
0276         prtOut->SetTheta ( std::atan2(prtOut->GetPt(),prtOut->GetPz()));
0277
0278       } else if(prtOut->IsThetaSmeared() && prtOut->IsPtSmeared()) /*(0110)*/ { // pt, theta
0279         // Note: Here (as in other cases), we rely on the device to deliver physically sound
0280         // numbers (see and adapt HandleBogusValues). But for now we explicitly fault on division by zero
0281         assert ( fabs(std::tan(prtOut->GetTheta())) > 1e-8 );
0282         prtOut->SetPz( prtOut->GetPt() / std::tan(prtOut->GetTheta()) );
0283         prtOut->SetP( std::sqrt(std::pow(prtOut->GetPt(), 2) + std::pow(prtOut->GetPz(), 2)));
0284
0285       } else if(prtOut->IsPSmeared() && prtOut->IsPzSmeared()) /*(1001)*/ { // P, pz
0286         if ( prtOut->GetP() < std::abs(prtOut->GetPz()) ) prtOut->SetP( std::abs(prtOut->GetPz()), false);
0287         prtOut->SetPt( std::sqrt(std::pow(prtOut->GetP(), 2) - std::pow(prtOut->GetPz(), 2)));
0288         prtOut->SetTheta( std::atan2(prtOut->GetPt() ,prtOut->GetPz() ) );
0289
0290       } else if(prtOut->IsThetaSmeared() && prtOut->IsPzSmeared())  /*(0101)*/ { // theta, pz
0291
0292         prtOut->SetPt( prtOut->GetPz() * std::tan(prtOut->GetTheta()));
0293         prtOut->SetP( std::sqrt( std::pow( prtOut->GetPt(),2) + std::pow( prtOut->GetPz(),2) ) );
0294       }
0295
0296       prtOut->SetPx (prtOut->GetP() * std::sin(prtOut->GetTheta()) * std::cos(prtOut->GetPhi()));
0297       prtOut->SetPy (prtOut->GetP() * std::sin(prtOut->GetTheta()) * std::sin(prtOut->GetPhi()));
0298
0299     } // Know px and py are unsmeared, phi IS smeared
0300
0301       } // case treatment for momentum components changed
0302
0303     } // LegacyMode
0304   } // if devices not empty
0305
0306   // Done.
0307   return prtOut;
0308
0309 }
0310
0311 std::vector<Smearer*> Detector::CopyDevices() const {
0312   std::vector<Smearer*> copies;
0313   for ( std::vector<Smearer*>::const_iterator it = Devices.begin();
0314     it!=Devices.end();
0315     ++it ){
0316     copies.push_back ( (*it)->Clone(""));
0317   }
0318   return copies;
0319 }
0320
0321 void Detector::Print(Option_t* o) const {
0322   for (unsigned i(0); i < GetNDevices(); ++i) {
0323     Devices.at(i)->Print(o);
0324   }  // for
0325 }
0326
0327   void Detector::SetLegacyMode( const bool mode ){
0328     LegacyMode = mode;
0329     if ( LegacyMode ){
0330       std::cout << "Warning: Turning on legacy mode, i.e. deactivating consistency checks and momentum regularization in Smear(). Use only for legacy smear scripts from earlier versions (<~1.0.4)" << endl;
0331     }
0332   }
0333
0334   bool Detector::GetLegacyMode() const {
0335     return LegacyMode;
0336   }
0337
0338 }  // namespace Smear