CompareToReroot.cxx

Go to the documentation of this file.

// $Id: CompareToReroot.cxx,v 1.8 2006/12/01 20:12:11 rhatcher Exp $
//
// Test module. Makes a tree of ScintHits and the PEs that they create
// by way of either PhotonTransport or GMINOS. Used to compare the two.
//
// n.tagg1@physics.ox.ac.uk
#include "CompareToReroot.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include <map>
#include "RerootExodus/RerootExodus.h"
#include "Plex/PlexHandle.h"
#include "MINF_Classes/MINFast.h"
#include "REROOT_Classes/REROOT_GeomMisc.h"
#include "REROOT_Classes/REROOT_PlanePos.h"
#include "REROOT_Classes/REROOT_PlaneSpec.h"
#include "REROOT_Classes/REROOT_CellPos.h"
#include "Digitization/DigiList.h"
#include "Digitization/DigiPE.h"
#include "Digitization/DigiScintHit.h"
#include "UgliGeometry/UgliGeomHandle.h"

#include "LogCounter.h"

ClassImp(HitCompare)

JOBMODULE(CompareToReroot, "CompareToReroot",
          "CompareToReroot");
CVSID("$Id: CompareToReroot.cxx,v 1.8 2006/12/01 20:12:11 rhatcher Exp $");
//......................................................................

CompareToReroot::CompareToReroot()
{
  // Open a file/tree
  fFile = new TFile("reroot_compare.root","RECREATE");
  fTree = new TTree("comp","comp");
  fHitCompare = new HitCompare;
  fTree->Branch("branch","HitCompare",&fHitCompare,32000,1);
}

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

CompareToReroot::~CompareToReroot()
{
  fFile->cd();
  fTree->Write("comp");
  fFile->Close();
  delete fFile;
}

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

JobCResult CompareToReroot::Ana(const MomNavigator* mom)
{
//======================================================================
// Fill entries in a tree.
//======================================================================

  // 
  SimSnarlRecord* simsnarl = 0;
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();
  PlexStripEndId seid;
  Float_t adc;
  Int_t   isAB;
  static int event=0;
  event++;

  // Get the simsnarl.
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }
  const SimSnarlHeader* simHeader = simsnarl->GetSimSnarlHeader();
  if(simHeader ==0){
    MSG("Photon",Msg::kError) << "Cannot find SimSnarlHeader in SimSnarl." << endl;
    return JobCResult::kFailed;
  }
  
  // We need a context to do table lookups
  VldContext simContext = simHeader->GetVldContext();
  PlexHandle plex(simContext);
  
  // Get the DigiScintHits.
  // Get the scint hit array.
  const TObjArray* hitArray = 
    dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiScintHits"));
  if(hitArray==0) {
    cout << "Can't find scint hit array.\n"; 
    return JobCResult::kPassed;  
  };
  
  
  // Build a map of stripend.
  std::map<int,HitCompare> hmap;
  std::map<int,HitCompare>::iterator hmapit;

  // put scint hits in.
  TIter hitarrayIter(hitArray);
  while( (tobj = hitarrayIter.Next()) ) {
    const DigiScintHit* scinthit = dynamic_cast<DigiScintHit*>(tobj);
    if(scinthit) {
      seid = scinthit->StripEndId();
      seid.SetEnd(StripEnd::kPositive);
      int index1 = seid.GetEncoded();
      seid.SetEnd(StripEnd::kNegative);
      int index2 = seid.GetEncoded();
      // put it on both ends.
      hmap[index1].nhits++;
      hmap[index1].hit_de += scinthit->DE();
      hmap[index1].hit_ds += scinthit->DS();
      hmap[index1].hit_x  = scinthit->X1();
      hmap[index1].hit_t  = scinthit->T1();

      hmap[index2].nhits++;
      hmap[index2].hit_de += scinthit->DE();
      hmap[index2].hit_ds += scinthit->DS();
      hmap[index2].hit_x  = scinthit->X1();
      hmap[index2].hit_t  = scinthit->T1();
    }
  }

  // put in the reroot stuff.
  // Prepare to loop over REROOT digits
  const TClonesArray *flsdigits = RerootExodus::GetFLSDigitList();
  REROOT_FLSDigit    *flsdigit  = 0;
  TIter diter(flsdigits);


  Float_t ADAMO_RNULL =  699050. * (Float_t) (16<<26); // 699050*16.0**26
  while ( ( flsdigit = (REROOT_FLSDigit*) diter.Next() ) ) {
    // process A side
    adc  = flsdigit->RawA();
    isAB = 0;
    if ( adc > 0.0 && adc < ADAMO_RNULL ) {
      seid = RerootExodus::PECAB2SEId(flsdigit->IPln(),flsdigit->IExtr(),
                                      flsdigit->ICell(),isAB);

      if(flsdigit->SignalPEA()>0) {
        int index = seid.GetEncoded();
        hmap[index].adc_reroot += adc;
        hmap[index].tdc_reroot = flsdigit->TDCA();
        hmap[index].pe_reroot += flsdigit->SignalPEA();
      }
    }
    
    // process B side
    adc  = flsdigit->RawB();
    isAB = 1;
    if ( adc > 0.0 && adc < ADAMO_RNULL ) {
      seid = RerootExodus::PECAB2SEId(flsdigit->IPln(),flsdigit->IExtr(),
                                      flsdigit->ICell(),isAB);      
      if(flsdigit->SignalPEB()>0) {
        int index = seid.GetEncoded();
        hmap[index].adc_reroot += adc;
        hmap[index].tdc_reroot = flsdigit->TDCB();
        hmap[index].pe_reroot += flsdigit->SignalPEB();
      }
    }
  }



  // Get the DigiPes.
  const TObjArray* peListObj = 
    dynamic_cast<const TObjArray*>(simsnarl->FindTemporary(0,"DigiListPe"));
  if(peListObj == 0) {
    MSG("DetSim",Msg::kError) << "Cannot find DigiPE list in the SimSnarl. Aborting." << endl;
    return JobCResult::kWarning;
  }
  
  TObjArrayIter itr(peListObj);
  DigiPE* digipe;
  while( (digipe = dynamic_cast<DigiPE*>(itr.Next()))!=0 ) {
    int index = plex.GetStripEndId(digipe->GetPixelSpotId()).GetEncoded();
    float t = digipe->GetTime()*1e9;
    if( t < (hmap[index].tfirst_photon)) hmap[index].tfirst_photon = t;
    hmap[index].tmean_photon += t;
    hmap[index].pe_photon +=1;
  };


  UgliGeomHandle ugli(simContext);

  // Finish up 
  for(hmapit = hmap.begin(); hmapit!=hmap.end(); hmapit++){
    seid = PlexStripEndId(hmapit->first);
    hmapit->second.event = event;
    hmapit->second.plane = seid.GetPlane();
    hmapit->second.strip = seid.GetStrip();
    hmapit->second.end   = seid.GetEnd();
    
    if(seid.IsValid()) {
      UgliStripHandle ustrip = ugli.GetStripHandle(seid); 
      if(ustrip.IsValid()) {
        hmapit->second.halflength = ustrip.GetHalfLength();
        hmapit->second.pigtail = ustrip.WlsPigtail(seid.GetEnd());
        hmapit->second.clear   = ustrip.ClearFiber(seid.GetEnd());
      }
    }
    if(hmapit->second.pe_photon >0)
      hmapit->second.tmean_photon /= (float) hmapit->second.pe_photon; 
  }
  
  // And fill.

  double tot_photon = 0;
  double tot_reroot = 0;
  for(hmapit = hmap.begin(); hmapit!=hmap.end(); hmapit++){
    seid = PlexStripEndId(hmapit->first);
    // fill other end.
    seid.SetEnd( (seid.GetEnd()==StripEnd::kPositive)?StripEnd::kNegative:StripEnd::kPositive);
    int index = seid.GetEncoded();
    hmapit->second.pe_reroot_2 = hmap[index].pe_reroot;
    hmapit->second.adc_reroot_2 = hmap[index].adc_reroot;
    hmapit->second.tdc_reroot_2 = hmap[index].tdc_reroot;
    hmapit->second.pe_photon_2 = hmap[index].pe_photon;
    hmapit->second.tfirst_photon_2 = hmap[index].tfirst_photon;
    hmapit->second.tmean_photon_2 = hmap[index].tmean_photon;
    
    tot_photon += hmapit->second.pe_photon;
    tot_reroot += hmapit->second.pe_reroot;
    fHitCompare = &(hmapit->second); 
    fTree->Fill();
  }

  static LogCounter counter;
  if(counter.Increment()) {
    fTree->AutoSave();
    cout << counter.fCount << endl;
  }
  fFile->Write();

  static double tot_photon_running = 0;
  static double tot_reroot_running = 0;
  tot_photon_running += tot_photon;
  tot_reroot_running += tot_reroot;
  cout << "PhotonTotal: " << tot_photon << "\t RunningTotal:" << tot_photon_running << endl;
  cout << "RerootTotal: " << tot_reroot << "\t RunningTotal:" << tot_reroot_running << endl;
  cout << "Ratio:       " << tot_photon/tot_reroot << "\t RunningTotal:" << tot_photon_running/tot_reroot_running << endl;

  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

---