SimpleTransport.cxx

Go to the documentation of this file.

// $Id: SimpleTransport.cxx,v 1.19 2009/01/14 19:15:27 brebel Exp $
//
// From Caius' code: "This is a really silly photon transporter that 
// just uses some silly fudge factor."
//
// jpaley@indiana.edu

#include "PhotonTransport/SimpleTransport.h"
#include "PhotonTransport/PhotonSignal.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/TruthFilter.h"

#include "JobControl/Exception.h"
#include "Geometry/CellUniqueId.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"
#include "RawData/RawDigit.h"
#include "RawData/DAQHeader.h"

//#include "ConnectionMap/ConnectionMap.h"

#include <cmath>
#include <iostream>

using namespace photrans;
using namespace std;

MODULE_DECL(SimpleTransport);

SimpleTransport::SimpleTransport(const char* version) : 
  PhotonTransporter(), 
  jobc::Module("SimpleTransport"),
  fMeVToNPhot(0), 
  fCollectionEff(0),
  fAttenLength(0)
{
  this->SetCfgVersion(version);
  
}

//............................................................

SimpleTransport::~SimpleTransport()
{
}

//............................................................

void SimpleTransport::Update(const cfg::Config& c)
{
  c("MeVToNPhot").   Get(fMeVToNPhot);
  c("CollectionEff").Get(fCollectionEff);
  c("AttenLength").  Get(fAttenLength);
  c("MessageLevel"). Get(msglevel);
}

//............................................................

jobc::Result SimpleTransport::Reco(edm::EventHandle& evt) 
{
  
  std::vector<const rawdata::DAQHeader*> header;
  short int det = rawdata::kFar;
  try{ evt.DAQ().Get("./",header); }
  catch(edm::Exception e){
    std::cerr << "Error retrieving daq header, while looking "<<
      "in SimpleTransport::Reco(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();
  geo::Geometry& fGeo = geo::Geometry::Instance(evt.Header().Run(), det);
  
  vector<const sim::FLSHitList*> hitlist(0);
  sim::TruthFilter::GetFilteredHits(evt,hitlist);
  vector<photrans::PhotonSignal> photlist(0);
  
  int nhits = 0;
  
  int planeId = -1;
  int cellId  = -1;

  float fGeVToMeV = 1000.;
  
  //geo::CellUniqueId lastuid;
  
  for (unsigned int i=0; i<hitlist.size(); ++i) {
    
    //  geo::CellUniqueId uid;
    
    float edep = 0.;
    double time = 0.;

    float dist=0.;

    unsigned short trkid = 0;    
     
    for (unsigned int j=0; j< hitlist[i]->fHits.size(); ++j) {
      ++nhits;
      
      //      uid = hitlist[i]->fHits[j].fId;
      
      //if(!(lastuid==uid))
      if(planeId!=hitlist[i]->fHits[j].fPlaneId || cellId!=hitlist[i]->fHits[j].fCellId)
        try{     
          
          planeId=hitlist[i]->fHits[j].fPlaneId;
          cellId=hitlist[i]->fHits[j].fCellId;
          //cout <<"plane = "<<planeId<<endl;
          //cout << "("<<planeId<<","<<cellId<<")    (" << hitlist[i]->fHits[j].fPlaneId<<","<<hitlist[i]->fHits[j].fCellId<<")"<<endl;
          //  if (&cell == 0) continue;
          //    cout << "cell = " << cellId << endl;
          //    cellHalfLength = cell.HalfL();

          //  cout << "cell Half length = " << cellHalfLength << endl;
          //    lastuid=uid;
        }catch(...){printf("bad photon geometry\n");continue;}
      
      
      edep = hitlist[i]->fHits[j].fEdep;
      time = hitlist[i]->fHits[j].fT;
      //      double pos1[3] = {0.,0.,hitlist[i]->fHits[j].fEntryZ()};
      //      double pos2[3];
      //      cell.WorldToLocal(pos1,pos2);

      dist = fGeo.
        Plane((unsigned int)planeId).
        Cell((unsigned int)cellId).
        DistToReadOut(0.5*(hitlist[i]->fHits[j].fEntryZ + 
                           hitlist[i]->fHits[j].fExitZ));
      trkid =  hitlist[i]->fHits[j].fTrackId;

      photrans::PhotonSignal 
        photon(time,int(edep*fMeVToNPhot*fGeVToMeV*fCollectionEff*exp(-dist/fAttenLength)));
      photon.SetPlaneCell(planeId, cellId);
      photon.SetTrackId(trkid);
      
      if(photon.NPhoton() > 0.) photlist.push_back(photon);
      
    }
    
  }

  if(msglevel>=INFO)cout << "Read " << nhits << " FLS hits" << endl;
  if(msglevel>=INFO)cout << "Created " << photlist.size() << " photon signals" << endl;
  
  for (unsigned int i=0; i<photlist.size(); ++i)
    evt.Raw().Put(photlist[i],"");
  
  photlist.clear();

  return jobc::kPassed;

}

---