ScintHitTruthModule.cxx

Go to the documentation of this file.

// $Id: ScintHitTruthModule.cxx,v 1.4 2007/03/01 17:37:56 rhatcher Exp $
//
// FILL_IN: [Document your code!!]
//
// kordosky@hep.utexas.edu
#include <cmath>

#include "TFile.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TObjArray.h"
#include "TObject.h"
#include "TIterator.h"

#include "Digitization/DigiScintHit.h"
#include "CalDetDST/ScintHitTruthModule.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"

using namespace std;

JOBMODULE(ScintHitTruthModule, "ScintHitTruthModule",
          "analyzes digiscinthits");
CVSID("$Id: ScintHitTruthModule.cxx,v 1.4 2007/03/01 17:37:56 rhatcher Exp $");
//......................................................................

ScintHitTruthModule::ScintHitTruthModule() 
   : fOut(0)
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

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

ScintHitTruthModule::~ScintHitTruthModule()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
}

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

void ScintHitTruthModule::BeginJob()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
   TDirectory* dsave = gDirectory;
   fOut = new TFile("scinttruth.root", "recreate");
   h_bvp = new TH2F("h_bvp", 
                    "z1 in scint plane vs. plane wtd. by dE;plane;z1 (cm)", 
                    60, -0.5, 59.5, 110, -0.55,0.55);
   h_evp = new TH2F("h_evp", 
                    "z2 in scint plane vs. plane wtd. by dE;plane,z2 (cm)", 
                    60, -0.5, 59.5, 110, -0.55,0.55);
   h_evds = new TH2F("h_evds","dE vs. dS;dS (mm);dE (MeV)", 500, 0, 20.0, 400,0,5.0);
   h_etot = new TH1F("h_etot", "Total Energy;Etot (GeV)", 700, 0, 0.7);

   h_dedsvpid = new TH2F("h_dedsvpid", "dE/dS v. pid;pid; dE/ds (MeV/cm)", 7,0,7, 500, 0, 100.0 );
   h_kevpid = new TH2F("h_kevpid", "log10(kE) v. pid; pid; log10(kE) (GeV)", 7,0,7, 700, -6, 1); 
   h_dedsvke_e = new TH2F("h_dedsvke_e", "e dE/dS v. log10(KE); log10(KE) (GeV); dE/dS (MeV/cm)", 700,-6,1, 500,0,100);
   h_dedsvke_mu = new TH2F("h_dedsvke_mu", "#mu dE/dS v. log10(KE); log10(KE) (GeV); dE/dS (MeV/cm)", 700,-6,1, 500,0,100);


   const char* pid_names[7]={"e","#mu","#gamma","#pi^{#pm}",
                             "#pi^{0}","p","n"};
   const char* pid_names2[7]={"e","mu","g","pi",
                             "pz","p","n"};

   for(int i=0; i<7; i++){
      h_dedsvpid->GetXaxis()->SetBinLabel(i+1,pid_names[i]);
      h_kevpid->GetXaxis()->SetBinLabel(i+1,pid_names[i]);
      TString name="h_edep_ps_"; name+=pid_names2[i];
      TString name2="h_edep_pske_"; name2+=pid_names2[i];

      TString title="energy deposited by "; title+=pid_names[i];
      title+=";plane;strip";

      TString title2="energy deposited by "; title2+=pid_names[i];
      title2+=";plane;strip;log10(KE)(GeV)";

      h_edep_ps[i]=new TH2F(name.Data(),title.Data(),
                            60,-0.5,59.5,24,-0.5,23.5);
      h_edep_pske[i]=new TH3F(name2.Data(), title2.Data(),
                              60,-0.5,59.5,24,-0.5,23.5,70,-6,1);
   }

   if(dsave) dsave->cd();

}

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

void ScintHitTruthModule::EndJob()
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
   fOut->Write();

}

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

JobCResult ScintHitTruthModule::Ana(const MomNavigator* mom)
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
   
   // get the digi scint hits from mom
   SimSnarlRecord* simsnarl = 0;
   TObject* tobj;
   TIter    fragiter = mom->FragmentIter();
  
   // Get the simsnarl.
   while( ( tobj = fragiter.Next() ) ) {
      simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
      if(simsnarl) break;
   }
  
   // Verify there IS a simsnarl.
   if(!simsnarl) {
      MSG("CalDetDST",Msg::kError) << "No SimSnarl found. You must run RerootToTruthModule()!" << endl;
      return JobCResult::kFailed;
   }


   const SimSnarlHeader* simHeader = simsnarl->GetSimSnarlHeader();
   if(simHeader ==0){
      MSG("CalDetDST",Msg::kError) << "Cannot find SimSnarlHeader in SimSnarl." << endl;
      return JobCResult::kFailed;
   }

   const TObjArray* hits = 
      dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiScintHits"));
     if(hits==0) {
        MSG("CalDetDST",Msg::kError) << "Can't find scint hit array.\n"; 
        return JobCResult::kFailed;  
     };


     static const int pidvec[7]={11,13,22,211,111,2212,2112};

     TObject* obj;
     TIter hitarrayIter(hits);
     double etot=0.0;
     double edep[60]={0.0};
     for(int i=0; i<60; i++){edep[i]=0.0;}
     while( (obj = hitarrayIter.Next()) ) {
        DigiScintHit* hit = dynamic_cast<DigiScintHit*>(obj);
        if(hit) {
           int plane= hit->Plane();
           int strip= hit->Strip();
           if( (plane<0) || (plane>59) ) continue;
           etot += hit->DE();      
           edep[plane]=hit->DE();
           h_bvp->Fill(plane,hit->Z1()*100.0,hit->DE()*1000.0);
           h_evp->Fill(plane,hit->Z2()*100.0,hit->DE()*1000.0);
           h_evds->Fill(hit->DS()*1000.0,hit->DE()*1000.0);
           int pid=abs(hit->ParticleId());
           int ipid=-1;
//         cout<<"pid: "<<pid<<endl;
           for(int i=0; i<7; i++){            
              if(pidvec[i]==pid) {ipid=i; break;}
           }
           if(ipid>=0){
              
              if(hit->DS()>0) h_dedsvpid->Fill(ipid,hit->DE()*10.0/hit->DS());
              float ke=hit->ParticleKineticEnergy();          
              if(ke>0) h_kevpid->Fill(ipid,log10(ke));
              if(hit->DS()>0 && ke>0){
                 if(ipid==0) 
                    h_dedsvke_e->Fill(log10(ke),hit->DE()*10.0/hit->DS());
                 if(ipid==1) 
                    h_dedsvke_mu->Fill(log10(ke),hit->DE()*10.0/hit->DS());
              }
              
              if(hit->DE()>0){
                 h_edep_ps[ipid]->Fill(plane,strip,hit->DE());
                 if(ke>0) 
                    h_edep_pske[ipid]->Fill(plane,strip,log10(ke),hit->DE());
              }
           }
           
        }
     }  
     h_etot->Fill(etot);

     


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

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

const Registry& ScintHitTruthModule::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("Mode", 1);
  r.LockValues();

  return r;
}

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

void ScintHitTruthModule::Config(const Registry& /* r */)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
//  int    tmpi;

//  if (r.Get("Mode",tmpi)) { fData = tmpi; }
}

---