loon: PlotVtx.cxx Source File

00001 // $Id: PlotVtx.cxx,v 1.13 2008/09/29 18:51:33 rustem Exp $
00002 
00003 // C++
00004 #include <cmath>
00005 #include <iostream>
00006 
00007 // ROOT
00008 #include "TDirectory.h"
00009 #include "TH1.h"
00010 #include "TH2.h"
00011 
00012 // MINOS
00013 #include "Registry/Registry.h"
00014 
00015 // Local
00016 #include "PhysicsNtuple/Default.h"
00017 #include "PhysicsNtuple/Factory.h"
00018 #include "PhysicsNtuple/Hist/HistMan.h"
00019 #include "PhysicsNtuple/Vertex.h"
00020 
00021 // Local
00022 #include "PlotVtx.h"
00023 
00024 REGISTER_ANP_OBJECT(AlgEvent,PlotVtx)
00025 REGISTER_ANP_OBJECT(AlgSnarl,PlotVtx)
00026 
00027 using namespace std;
00028 
00029 //----------------------------------------------------------------------
00030 Anp::HistVtx::HistVtx()
00031    :fPlot(false),
00032     fNMiss(0),
00033     fCenterX(1.4828),
00034     fCenterY(0.2384),
00035     fDir(0),
00036     fX(0),
00037     fY(0),
00038     fU(0),
00039     fV(0),
00040     fZ(0),
00041     fR2(0),
00042     fCosX(0),
00043     fCosY(0),
00044     fCosZ(0),
00045     fYX(0),
00046     fXZ(0),
00047     fYZ(0),
00048     fDiffX(0),
00049     fDiffY(0),
00050     fDiffZ(0),
00051     fDiffR2d(0),
00052     fDiffR3d(0),
00053     fProjRad(0),
00054     fProjAng(0),
00055     fDiffR2dForYX(0),
00056     fDiffR2dForXZ(0),
00057     fDiffR2dForYZ(0),
00058     fProjRadVsRad(0)
00059 {
00060 }
00061 
00062 //----------------------------------------------------------------------
00063 Anp::HistVtx::~HistVtx()
00064 { 
00065 }
00066 
00067 //----------------------------------------------------------------------
00068 bool Anp::HistVtx::Fill(const Vertex &vertex, const double weight)
00069 {
00070    if(!fPlot || !fDir) return false;
00071    
00072    fX -> Fill(vertex.X(), weight);
00073    fY -> Fill(vertex.Y(), weight);
00074    fU -> Fill(vertex.U(), weight);
00075    fV -> Fill(vertex.V(), weight);
00076    fZ -> Fill(vertex.Z(), weight);
00077 
00078    fR2 -> Fill(vertex.X()*vertex.X() + vertex.Y()*vertex.Y(), weight);
00079    
00080    fCosX -> Fill(vertex.CosX(), weight);
00081    fCosY -> Fill(vertex.CosY(), weight);
00082    fCosZ -> Fill(vertex.CosZ(), weight);
00083    
00084    if(fYX) fYX -> Fill(vertex.X(), vertex.Y(), weight);
00085    if(fXZ) fXZ -> Fill(vertex.Z(), vertex.X(), weight);
00086    if(fYZ) fYZ -> Fill(vertex.Z(), vertex.Y(), weight);
00087 
00088    return true;
00089 }
00090 
00091 
00092 //----------------------------------------------------------------------
00093 bool Anp::HistVtx::Fill(const Vertex &lhs, const Vertex &rhs, const double weight)
00094 {
00095    if(!fPlot || !fDir) return false;
00096   
00097    const double diffX = lhs.X() - rhs.X();
00098    const double diffY = lhs.Y() - rhs.Y();
00099    const double diffZ = lhs.Z() - rhs.Z();
00100    const double diffR2d = std::sqrt(diffX*diffX + diffY*diffY);
00101    const double diffR3d = std::sqrt(diffX*diffX + diffY*diffY + diffZ*diffZ);
00102 
00103    fDiffX -> Fill(diffX, weight);
00104    fDiffY -> Fill(diffY, weight);
00105    fDiffZ -> Fill(diffZ, weight);
00106 
00107    fDiffR2d -> Fill(diffR2d, weight);
00108    fDiffR3d -> Fill(diffR3d, weight);
00109 
00110    //
00111    // Compute 2d (xy) vector between lhs vertex and fiducial center
00112    //
00113    const double coordX = lhs.X() - fCenterX;
00114    const double coordY = lhs.Y() - fCenterY;
00115    const double coordR = std::sqrt(coordX*coordX + coordY*coordY);
00116    
00117    if(coordR > 0.0)
00118    {
00119       fProjRad -> Fill((+coordX*diffX + coordY*diffY)/coordR, weight);
00120       fProjAng -> Fill((-coordY*diffX + coordX*diffY)/coordR, weight);
00121       
00122       if(fProjRadVsRad) fProjRadVsRad -> Fill(coordR, (+coordX*diffX + coordY*diffY)/coordR, weight);
00123    }
00124    else
00125    {
00126       fProjRad -> Fill(0.0, weight);
00127       fProjAng -> Fill(0.0, weight);      
00128    }
00129    
00130    if(fDiffR2dForYX) fDiffR2dForYX -> Fill(lhs.X(), lhs.Y(), diffR2d);
00131    if(fDiffR2dForXZ) fDiffR2dForXZ -> Fill(lhs.Z(), lhs.X(), diffR2d);
00132    if(fDiffR2dForYZ) fDiffR2dForYZ -> Fill(lhs.Z(), lhs.Y(), diffR2d);
00133 
00134    return true;
00135 }
00136 
00137 //----------------------------------------------------------------------
00138 bool Anp::HistVtx::Make(TDirectory *dir, const string &option)
00139 {
00140    if(!dir) return false;
00141 
00142    fDir  = dir;
00143    fPlot = false;
00144 
00145    fX    = HistVtx::GetTH1("x");
00146    fY    = HistVtx::GetTH1("y");
00147    fZ    = HistVtx::GetTH1("z");
00148    fU    = HistVtx::GetTH1("u");
00149    fV    = HistVtx::GetTH1("v");   
00150    fR2   = HistVtx::GetTH1("r2");   
00151    fCosX = HistVtx::GetTH1("cosx");
00152    fCosY = HistVtx::GetTH1("cosy");
00153    fCosZ = HistVtx::GetTH1("cosz");
00154 
00155    if(option.find("plot_2d") != string::npos)
00156    {
00157       fYX = HistVtx::GetTH2("yx");
00158       fXZ = HistVtx::GetTH2("xz");
00159       fYZ = HistVtx::GetTH2("yz");
00160    }
00161 
00162    if(option.find("plot_diff") != string::npos)
00163    {
00164       fDiffX   = HistVtx::GetTH1("diff_x");
00165       fDiffY   = HistVtx::GetTH1("diff_y");
00166       fDiffZ   = HistVtx::GetTH1("diff_z");      
00167       fDiffR2d = HistVtx::GetTH1("diff_r2d");
00168       fDiffR3d = HistVtx::GetTH1("diff_r3d");      
00169       fProjRad = HistVtx::GetTH1("proj_rad");
00170       fProjAng = HistVtx::GetTH1("proj_ang");
00171 
00172       if(option.find("plot_2d") != string::npos)
00173       {
00174          fDiffR2dForYX = HistVtx::GetTH2("diff_r2d_yx");
00175          fDiffR2dForXZ = HistVtx::GetTH2("diff_r2d_xz");
00176          fDiffR2dForYZ = HistVtx::GetTH2("diff_r2d_yz");
00177          fProjRadVsRad = HistVtx::GetTH2("proj_rad_vs_rad");
00178       }
00179    }
00180 
00181    if(fNMiss == 0)
00182    {
00183       fPlot = true;
00184    }
00185    else
00186    {
00187       cerr << "HistVtx::Make - missed " << fNMiss << " histograms" << endl;
00188    }
00189 
00190    return fPlot;
00191 }
00192 
00193 //----------------------------------------------------------------------
00194 void Anp::HistVtx::SetCenter(const double x, const double y)
00195 {
00196    fCenterX = x;
00197    fCenterY = y;
00198 }
00199 
00200 //-----------------------------------------------------------------------------
00201 TH1* Anp::HistVtx::GetTH1(const std::string &key, const std::string &name)
00202 {
00203    TH1 *h = HistMan::Instance().CreateTH1(key, "vertex");
00204    if(h)
00205    {
00206       Anp::SetDir(h, fDir, name);
00207    }
00208    else
00209    {
00210       ++fNMiss;
00211    }
00212 
00213    return h;
00214 }
00215 
00216 //-----------------------------------------------------------------------------
00217 TH2* Anp::HistVtx::GetTH2(const std::string &key, const std::string &name)
00218 {
00219    TH2 *h = HistMan::Instance().CreateTH2(key, "vertex");
00220    if(h)
00221    {
00222       Anp::SetDir(h, fDir, name);
00223    }
00224    else
00225    {
00226       ++fNMiss;
00227    }
00228 
00229    return h;
00230 }
00231 
00232 //----------------------------------------------------------------------
00233 Anp::PlotVtx::PlotVtx()
00234    :fDirName("vertex"),
00235     fInput("all"),
00236     fDir(0),
00237     fInit(false),
00238     fPlot2d(true),
00239     fPlotDiff(false),
00240     fCenterX(1.4885),
00241     fCenterY(0.1397)
00242 {
00243 }
00244 
00245 //----------------------------------------------------------------------
00246 Anp::PlotVtx::~PlotVtx()
00247 {
00248 }
00249 
00250 //----------------------------------------------------------------------
00251 bool Anp::PlotVtx::Run(Record &record)
00252 {
00253    //
00254    // Plot vertex information for events, tracks and showers
00255    //
00256    if(!fDir) return true;
00257 
00258    for(EventIterator ievent = record.EventBegIterator(); ievent != record.EventEndIterator(); ++ievent)
00259    {      
00260       PlotVtx::Run(*ievent, record, true);
00261    }
00262 
00263    for(TrackIter itrack = record.TrackBeg(); itrack != record.TrackEnd(); ++itrack)
00264    {
00265       PlotVtx::Fill("tracks_beg", itrack -> GetBegVtx(), itrack -> Weight());
00266       PlotVtx::Fill("tracks_end", itrack -> GetEndVtx(), itrack -> Weight());
00267    }
00268    
00269    for(TruthIter itruth = record.TruthBeg(); itruth != record.TruthEnd(); ++itruth)
00270    {
00271       PlotVtx::Fill("truths", itruth -> GetVtx(), itruth -> Weight());
00272    }
00273 
00274    if(fPlotDiff && !record.GetHeader().IsData())
00275    {
00276       for(TrackIter itrack = record.TrackBeg(); itrack != record.TrackEnd(); ++itrack)
00277       {
00278          const TruthIter itruth = record.FindTruth(*itrack);
00279          if(itruth != record.TruthEnd())
00280          {
00281             PlotVtx::Fill("truth_minus_track", itruth->GetVtx(), itrack->Weight());
00282             PlotVtx::Fill("truth_minus_track", itruth->GetVtx(), itrack->GetBegVtx(),itrack->Weight());
00283          }
00284       }
00285    }
00286 
00287    return true;
00288 }
00289 
00290 //----------------------------------------------------------------------
00291 bool Anp::PlotVtx::Run(Event &event, const Record &record, bool pass)
00292 {
00293    //
00294    // Plot event vertex 
00295    //
00296    if(!pass || !fDir) return true;
00297 
00298    PlotVtx::Fill("event", event.GetVtx(), event.Weight());
00299 
00300    const TrackIter  itrack  = Anp::LongestTrack (event, record);
00301    const ShowerIter ishower = Anp::PrimaryShower(event, record);
00302 
00303    if(ishower != record.ShowerEnd())
00304    {
00305       PlotVtx::Fill("event_shower_prm", ishower -> GetVtx(), event.Weight());
00306    }
00307 
00308    if(itrack != record.TrackEnd())
00309    {
00310       PlotVtx::Fill("event_track_beg", itrack -> GetBegVtx(), event.Weight());
00311       PlotVtx::Fill("event_track_end", itrack -> GetEndVtx(), event.Weight());
00312 
00313       if(ishower != record.ShowerEnd())
00314       {
00315          PlotVtx::Fill("event_shower_prm", itrack -> GetEndVtx(), ishower -> GetVtx(), event.Weight());
00316       }
00317    }
00318    
00319    return true;
00320 }
00321 
00322 //----------------------------------------------------------------------
00323 void Anp::PlotVtx::Config(const Registry &reg)
00324 {
00325    const char *value_char = 0;
00326    if(reg.Get("PlotVtxDirName", value_char) && value_char)
00327    {
00328       fDirName = value_char;
00329    }
00330 
00331    value_char = 0;
00332    if(reg.Get("PlotVtxInput", value_char) && value_char)
00333    {
00334       fInput = value_char;
00335    }
00336    
00337    Anp::Read(reg, "PlotVtx2d",   fPlot2d);
00338    Anp::Read(reg, "PlotVtxDiff", fPlotDiff);
00339 
00340    reg.Get("PlotVtxCenterX", fCenterX);
00341    reg.Get("PlotVtxCenterY", fCenterY);
00342 
00343    if(reg.KeyExists("PrintConfig"))
00344    {
00345       cout << "PlotVtx::Config" << endl
00346            << "   DirName = " << fDirName << endl
00347            << "   Input = " << fInput << endl
00348            << "   Plot2d = " << fPlot2d << endl
00349            << "   PlotDiff = " << fPlotDiff << endl
00350            << "   CenterX = " << fCenterX << endl
00351            << "   CenterY = " << fCenterY << endl;
00352    }
00353 }
00354 
00355 //----------------------------------------------------------------------
00356 void Anp::PlotVtx::Set(TDirectory *dir)
00357 {
00358    fInit = false;
00359 
00360    if(!dir)
00361    {
00362       return;
00363    }
00364 
00365    fDir = Anp::GetDir(fDirName, dir);
00366 }
00367 
00368 //----------------------------------------------------------------------
00369 void Anp::PlotVtx::End(const DataBlock &)
00370 {
00371    fMap.clear();
00372    fDir = 0;
00373 }
00374 
00375 //----------------------------------------------------------------------
00376 bool Anp::PlotVtx::Fill(const string &input, const Vertex &vtx, const double weight)
00377 {
00378    if(fInput != "all" && fInput != input)
00379    {
00380       return false;
00381    }
00382 
00383    Handle<HistVtx> hist = PlotVtx::GetPlot(input);
00384    if(hist.valid())
00385    {
00386       return hist -> Fill(vtx, weight);
00387    }
00388 
00389    return false;
00390 }
00391 
00392 
00393 //----------------------------------------------------------------------
00394 bool Anp::PlotVtx::Fill(const string &input, const Vertex &lhs, const Vertex &rhs, const double weight)
00395 {
00396    if(!fPlotDiff)
00397    {
00398       return false;
00399    }
00400 
00401    if(fInput != "all" && fInput != input)
00402    {
00403       return false;
00404    }
00405 
00406    Handle<HistVtx> hist = PlotVtx::GetPlot(input);
00407    if(hist.valid())
00408    {
00409       return hist -> Fill(lhs, rhs, weight);
00410    }
00411 
00412    return false;
00413 }
00414 
00415 //----------------------------------------------------------------------
00416 Anp::Handle<Anp::HistVtx> Anp::PlotVtx::GetPlot(const string &input)
00417 {
00418    //
00419    // Find, create first, histogram object
00420    //
00421 
00422    PlotMap::iterator hit = fMap.find(input);
00423    if(hit == fMap.end())
00424    {
00425       Handle<HistVtx> hist(new HistVtx());
00426       TDirectory *dir = Anp::GetDir(input, fDir);
00427 
00428       if(dir)
00429       {
00430          hist -> SetCenter(fCenterX, fCenterY);
00431          
00432          string option;
00433          if(fPlot2d)   option += "plot_2d ";
00434          if(fPlotDiff) option += "plot_diff ";
00435          
00436          //
00437          // Attempt to create histograms
00438          //
00439          if(!(hist -> Make(dir, option))) hist.release();
00440       }
00441       else
00442       {
00443          hist.release();
00444       }
00445 
00446       hit = fMap.insert(PlotMap::value_type(input, hist)).first;
00447    }
00448 
00449    return hit -> second;
00450 }