loon: FillMadId.cxx Source File

00001 // $Id: FillMadId.cxx,v 1.3 2008/04/29 21:07:16 rustem Exp $
00002 
00003 // C++
00004 #include <fstream>
00005 
00006 // MINOS
00007 #include "CandNtupleSR/NtpSREvent.h"
00008 #include "CandNtupleSR/NtpSRTrack.h"
00009 #include "MessageService/MsgService.h"
00010 #include "MinosObjectMap/MomNavigator.h"
00011 #include "Registry/Registry.h"
00012 #include "StandardNtuple/NtpStRecord.h"
00013 
00014 // MINOS
00015 #include "Mad/MadAbID.h"
00016 #include "Mad/MadDpID.h"
00017 #include "Registry/Registry.h"
00018 
00019 // Local
00020 #include "PhysicsNtuple/Default.h"
00021 #include "PhysicsNtuple/Factory.h"
00022 #include "FillMadId.h"
00023 
00024 REGISTER_ANP_OBJECT(AlgStore,FillMadId)
00025 
00026 CVSID("$Id: FillMadId.cxx,v 1.3 2008/04/29 21:07:16 rustem Exp $");
00027 
00028 using namespace std;
00029 
00030 //------------------------------------------------------------------------------------------
00031 Anp::FillMadId::FillMadId()
00032    :fKeyBaseA(7090),
00033     fKeyBaseD(7080),
00034     fInit(false),
00035     fFillAll(false),
00036     fModuleA(new MadAbID()),
00037     fModuleD(new MadDpID()),
00038     fBeam(BeamType::kUnknown)
00039 { 
00040 }
00041 
00042 //------------------------------------------------------------------------------------------
00043 Anp::FillMadId::~FillMadId()
00044 {
00045    delete fModuleA;
00046    delete fModuleD;
00047 }
00048 
00049 //-----------------------------------------------------------------------------
00050 bool Anp::FillMadId::Run(Record& record, TObject *ptr)
00051 {
00052    if(!Reset(record))
00053    {
00054       return true;
00055    }
00056 
00057    NtpStRecord *ntprec = dynamic_cast<NtpStRecord *>(ptr);
00058    if(!ntprec)
00059    {
00060       const MomNavigator *mom = dynamic_cast<const MomNavigator *> (ptr);
00061       if(mom)
00062       {
00063          ntprec = dynamic_cast<NtpStRecord *>(mom -> GetFragment("NtpStRecord")); 
00064       }
00065       else
00066       {
00067          MSG("FillAlg", Msg::kError) << "Failed to find MomNavigator pointer" << endl;
00068          return false;
00069       }      
00070    }
00071 
00072    if(!ntprec)
00073    {
00074       MSG("FillAlg", Msg::kError) << "Failed to get NtpStRecord pointer" << endl;
00075       return false;
00076    }
00077 
00078    TClonesArray *event_array = ntprec -> evt;
00079    TClonesArray *track_array = ntprec -> trk;
00080    if(!event_array || !track_array)
00081    {           
00082       MSG("FillAlg", Msg::kWarning) << "TClonesArray of events or tracks is invalid" << endl;
00083       return false;
00084    }   
00085 
00086    map<int, NtpSRTrack *> tmap;
00087 
00088    const int ntrack = track_array -> GetEntries();
00089    for(int i = 0; i < ntrack; ++i)
00090    {
00091       NtpSRTrack *ntptrk  = dynamic_cast<NtpSRTrack *>(track_array -> At(i));
00092       if(!ntptrk)
00093       {
00094          MSG("FillAlg", Msg::kWarning) << "Failed to find NtpSRTrack at index " << i << endl;
00095          continue;
00096       }
00097       
00098       if(!tmap.insert(map<int, NtpSRTrack *>::value_type(ntptrk->index, ntptrk)).second)
00099       {
00100          MSG("FillAlg", Msg::kError) << "Failed to insert track: " << ntptrk->index << endl;
00101       }
00102    }
00103 
00104    const int nevent = event_array -> GetEntries();
00105 
00106    for(int ievent = 0; ievent < nevent; ++ievent)
00107    {      
00108       NtpSREvent* ntpevt = dynamic_cast<NtpSREvent *>(event_array -> At(ievent));   
00109       if(!ntpevt)
00110       {
00111          MSG("FillAlg", Msg::kWarning) << "Failed to find NtpSREvent at "<< ievent << endl;
00112          continue; 
00113       }
00114       
00115       if(ntpevt -> index != ievent)
00116       {
00117          MSG("FillAlg", Msg::kError) << "Mismatched event index" << endl;
00118          continue;
00119       }
00120 
00121       const double pid             = fModuleA -> CalcPID(ntpevt, ntprec);
00122       const double costheta        = fModuleA -> GetRecoCosTheta(ntpevt,ntprec);
00123       const double eventenergy     = fModuleA -> GetRecoE(ntpevt,ntprec);
00124       const double trackcharge     = fModuleA -> GetTrackEMCharge(ntpevt,ntprec);
00125       const double trackenergy     = fModuleA -> GetTrackPlanes(ntpevt,ntprec);
00126       const double tracklikeplanes = fModuleA -> GetTrackLikePlanes(ntpevt,ntprec);
00127       const double trackphfrac     = fModuleA -> GetTrackPHfrac(ntpevt,ntprec);
00128       const double trackphmean     = fModuleA -> GetTrackPHmean(ntpevt,ntprec);
00129       const double trackqpsigmaqp  = fModuleA -> GetTrackQPsigmaQP(ntpevt,ntprec);
00130       const double recoy           = fModuleA -> GetRecoY(ntpevt,ntprec);
00131 
00132       for(EventIterator event = record.EventBegIterator(); event != record.EventEndIterator(); ++event)
00133       {
00134          if(event -> EventIndex() != ievent)
00135          {
00136             continue;
00137          }
00138 
00139          if(fKeyBaseA > 0)
00140          {
00141             event -> Add(fKeyBaseA + 0, pid);
00142             
00143             if(fFillAll)
00144             {
00145                event -> Add(fKeyBaseA + 1, costheta);
00146                event -> Add(fKeyBaseA + 2, eventenergy);
00147                event -> Add(fKeyBaseA + 3, trackcharge);
00148                event -> Add(fKeyBaseA + 4, trackenergy);
00149                event -> Add(fKeyBaseA + 5, tracklikeplanes);
00150                event -> Add(fKeyBaseA + 6, trackphfrac);
00151                event -> Add(fKeyBaseA + 7, trackphmean);
00152                event -> Add(fKeyBaseA + 8, trackqpsigmaqp);
00153                event -> Add(fKeyBaseA + 9, recoy);
00154             }
00155          }
00156 
00157          if(fKeyBaseD > 0)
00158          {
00159             const BeamType::BeamType_t beam = BeamType::TagToEnum(record.GetHeader().GetBEAMTYPE().c_str());
00160             if(beam == BeamType::kUnknown)
00161             {
00162                MSG("FillAlg", Msg::kWarning) << "Unknown beam type" << endl;
00163                continue;
00164             }
00165 
00166             Detector::Detector_t det = Detector::kUnknown;
00167             if     (record.GetHeader().IsNear()) det = Detector::kNear;
00168             else if(record.GetHeader().IsFar())  det = Detector::kFar;
00169             else continue;
00170             
00171             const TrackIter itrack = Anp::LongestTrack(*event, record);
00172             if(itrack == record.TrackEnd())
00173             {
00174                continue;
00175             }
00176 
00177             const map<int, NtpSRTrack *>::const_iterator strack = tmap.find(itrack->TrackIndex());
00178             if(strack == tmap.end())
00179             {
00180                MSG("FillAlg", Msg::kError) << "Failed to find track" << endl;
00181                continue;               
00182             }
00183 
00184             if(fModuleD -> ChoosePDFs(det,
00185                                       beam,
00186                                       ReleaseType::AsString(record.GetHeader().Release()),
00187                                       "daikon"))
00188             {
00189                const float dpid = fModuleD -> CalcPID(strack->second, ntpevt, &(ntprec->evthdr), det, 0);
00190                event -> Add(fKeyBaseA + 0, dpid);
00191             }
00192          }
00193       }
00194    }
00195 
00196    return true;
00197 }
00198 
00199 //-----------------------------------------------------------------------------
00200 void Anp::FillMadId::Config(const Registry &reg)
00201 {
00202    reg.Get("FillMadIdKeyBaseA", fKeyBaseA);
00203    reg.Get("FillMadIdKeyBaseA", fKeyBaseD);
00204    
00205    Anp::Read(reg, "FillMadIdAll", fFillAll);
00206 
00207    if(reg.KeyExists("PrintConfig"))
00208    {
00209       cout << "FillMadId::Config" << endl
00210            << "   FillAll = " << fFillAll << endl
00211            << "   KeyBaseA = " << fKeyBaseA << endl
00212            << "   KeyBaseD = " << fKeyBaseD << endl;
00213    }
00214 }
00215 
00216 //-----------------------------------------------------------------------------
00217 bool Anp::FillMadId::Reset(const Record &record)
00218 {
00219    const BeamType::BeamType_t beam = BeamType::TagToEnum(record.GetHeader().GetBEAMTYPE().c_str());
00220    if(beam == BeamType::kUnknown)
00221    {
00222       MSG("FillAlg", Msg::kWarning) << "Unknown beam type" << endl;
00223       return false;
00224    }
00225 
00226    if(beam == fBeam && fInit)
00227    {
00228       return true;
00229    }
00230 
00231    fInit = false;
00232 
00233    string wfile;
00234    if(record.GetHeader().IsFar())
00235    {
00236       if(beam == BeamType::kL010z185i)
00237       {
00238          wfile = "ab_pdf_far_le_cedar_daikon.root";
00239       }
00240       else if(beam == BeamType::kL250z200i)
00241       {
00242          wfile = "ab_pdf_far_phe_cedar_daikon.root";
00243       }
00244    }
00245    else if(record.GetHeader().IsNear())
00246    {
00247       if(beam == BeamType::kL010z000i)
00248       {
00249          wfile = "ab_pdf_near_le0_cedar_daikon.root";
00250       }
00251       else if(beam == BeamType::kL010z170i)
00252       {
00253          wfile = "ab_pdf_near_le170_cedar_daikon.root";
00254       }
00255       else if(beam == BeamType::kL010z200i)
00256       {
00257          wfile = "ab_pdf_near_le200_cedar_daikon.root";
00258       }
00259       else if(beam == BeamType::kL010z185i)
00260       {
00261          wfile = "ab_pdf_near_le_cedar_daikon.root";
00262       }
00263       else if(beam == BeamType::kL250z200i)
00264       {
00265          wfile = "ab_pdf_near_phe_cedar_daikon.root";
00266       }
00267       else if(beam == BeamType::kL100z200i)
00268       {
00269          wfile = "ab_pdf_near_pme_cedar_daikon.root";
00270       }
00271       else if(beam == BeamType::kL150z200i)
00272       {
00273          wfile = "ab_pdf_near_pmhe_cedar_daikon.root";
00274       }
00275    }
00276 
00277    if(wfile.empty())
00278    {
00279       MSG("FillAlg", Msg::kWarning) << "Unknown weight file name" << endl;
00280       return false;
00281    }
00282 
00283    string mpath = std::getenv("SRT_PRIVATE_CONTEXT");
00284    
00285    std::ifstream ufile((mpath + "/Mad/data/" + wfile).c_str());
00286    if(!ufile || !ufile.is_open())
00287    {
00288       mpath = std::getenv("SRT_PUBLIC_CONTEXT");
00289 
00290       std::ifstream ifile((mpath + "/Mad/data/" + wfile).c_str());
00291       if(!ifile || !ifile.is_open())
00292       {
00293          MSG("FillAlg", Msg::kWarning) << "FillMadId::Init - weight file does not exist "<< endl;
00294          return false;
00295       }
00296    }
00297 
00298    fInit = true;
00299    fBeam = beam;
00300 
00301    fModuleA -> ReadPDFs((mpath + "/Mad/data/" + wfile).c_str());
00302 
00303    MSG("FillAlg", Msg::kDebug) << "Read new weight file for " << BeamType::AsString(beam) << endl 
00304                                << mpath + "/Mad/data/" + wfile << endl;
00305    
00306    return true;
00307 }