loon: ScintHitTruthModule class Reference

00035    : fOut(0)
00036 {
00037 //======================================================================
00038 // FILL_IN: [Document your code!!]
00039 //======================================================================
00040 }

00045 {
00046 //======================================================================
00047 // FILL_IN: [Document your code!!]
00048 //======================================================================
00049 }

00116 {
00117 //======================================================================
00118 // FILL_IN: [Document your code!!]
00119 //======================================================================
00120    
00121    // get the digi scint hits from mom
00122    SimSnarlRecord* simsnarl = 0;
00123    TObject* tobj;
00124    TIter    fragiter = mom->FragmentIter();
00125   
00126    // Get the simsnarl.
00127    while( ( tobj = fragiter.Next() ) ) {
00128       simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
00129       if(simsnarl) break;
00130    }
00131   
00132    // Verify there IS a simsnarl.
00133    if(!simsnarl) {
00134       MSG("CalDetDST",Msg::kError) << "No SimSnarl found. You must run RerootToTruthModule()!" << endl;
00135       return JobCResult::kFailed;
00136    }
00137 
00138 
00139    const SimSnarlHeader* simHeader = simsnarl->GetSimSnarlHeader();
00140    if(simHeader ==0){
00141       MSG("CalDetDST",Msg::kError) << "Cannot find SimSnarlHeader in SimSnarl." << endl;
00142       return JobCResult::kFailed;
00143    }
00144 
00145    const TObjArray* hits = 
00146       dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiScintHits"));
00147      if(hits==0) {
00148         MSG("CalDetDST",Msg::kError) << "Can't find scint hit array.\n"; 
00149         return JobCResult::kFailed;  
00150      };
00151 
00152 
00153      static const int pidvec[7]={11,13,22,211,111,2212,2112};
00154 
00155      TObject* obj;
00156      TIter hitarrayIter(hits);
00157      double etot=0.0;
00158      double edep[60]={0.0};
00159      for(int i=0; i<60; i++){edep[i]=0.0;}
00160      while( (obj = hitarrayIter.Next()) ) {
00161         DigiScintHit* hit = dynamic_cast<DigiScintHit*>(obj);
00162         if(hit) {
00163            int plane= hit->Plane();
00164            int strip= hit->Strip();
00165            if( (plane<0) || (plane>59) ) continue;
00166            etot += hit->DE();      
00167            edep[plane]=hit->DE();
00168            h_bvp->Fill(plane,hit->Z1()*100.0,hit->DE()*1000.0);
00169            h_evp->Fill(plane,hit->Z2()*100.0,hit->DE()*1000.0);
00170            h_evds->Fill(hit->DS()*1000.0,hit->DE()*1000.0);
00171            int pid=abs(hit->ParticleId());
00172            int ipid=-1;
00173 //         cout<<"pid: "<<pid<<endl;
00174            for(int i=0; i<7; i++){            
00175               if(pidvec[i]==pid) {ipid=i; break;}
00176            }
00177            if(ipid>=0){
00178               
00179               if(hit->DS()>0) h_dedsvpid->Fill(ipid,hit->DE()*10.0/hit->DS());
00180               float ke=hit->ParticleKineticEnergy();          
00181               if(ke>0) h_kevpid->Fill(ipid,log10(ke));
00182               if(hit->DS()>0 && ke>0){
00183                  if(ipid==0) 
00184                     h_dedsvke_e->Fill(log10(ke),hit->DE()*10.0/hit->DS());
00185                  if(ipid==1) 
00186                     h_dedsvke_mu->Fill(log10(ke),hit->DE()*10.0/hit->DS());
00187               }
00188               
00189               if(hit->DE()>0){
00190                  h_edep_ps[ipid]->Fill(plane,strip,hit->DE());
00191                  if(ke>0) 
00192                     h_edep_pske[ipid]->Fill(plane,strip,log10(ke),hit->DE());
00193               }
00194            }
00195            
00196         }
00197      }  
00198      h_etot->Fill(etot);
00199 
00200      
00201 
00202 
00203   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
00204 }

00054 {
00055 //======================================================================
00056 // FILL_IN: [Document your code!!]
00057 //======================================================================
00058    TDirectory* dsave = gDirectory;
00059    fOut = new TFile("scinttruth.root", "recreate");
00060    h_bvp = new TH2F("h_bvp", 
00061                     "z1 in scint plane vs. plane wtd. by dE;plane;z1 (cm)", 
00062                     60, -0.5, 59.5, 110, -0.55,0.55);
00063    h_evp = new TH2F("h_evp", 
00064                     "z2 in scint plane vs. plane wtd. by dE;plane,z2 (cm)", 
00065                     60, -0.5, 59.5, 110, -0.55,0.55);
00066    h_evds = new TH2F("h_evds","dE vs. dS;dS (mm);dE (MeV)", 500, 0, 20.0, 400,0,5.0);
00067    h_etot = new TH1F("h_etot", "Total Energy;Etot (GeV)", 700, 0, 0.7);
00068 
00069    h_dedsvpid = new TH2F("h_dedsvpid", "dE/dS v. pid;pid; dE/ds (MeV/cm)", 7,0,7, 500, 0, 100.0 );
00070    h_kevpid = new TH2F("h_kevpid", "log10(kE) v. pid; pid; log10(kE) (GeV)", 7,0,7, 700, -6, 1); 
00071    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);
00072    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);
00073 
00074 
00075    const char* pid_names[7]={"e","#mu","#gamma","#pi^{#pm}",
00076                              "#pi^{0}","p","n"};
00077    const char* pid_names2[7]={"e","mu","g","pi",
00078                              "pz","p","n"};
00079 
00080    for(int i=0; i<7; i++){
00081       h_dedsvpid->GetXaxis()->SetBinLabel(i+1,pid_names[i]);
00082       h_kevpid->GetXaxis()->SetBinLabel(i+1,pid_names[i]);
00083       TString name="h_edep_ps_"; name+=pid_names2[i];
00084       TString name2="h_edep_pske_"; name2+=pid_names2[i];
00085 
00086       TString title="energy deposited by "; title+=pid_names[i];
00087       title+=";plane;strip";
00088 
00089       TString title2="energy deposited by "; title2+=pid_names[i];
00090       title2+=";plane;strip;log10(KE)(GeV)";
00091 
00092       h_edep_ps[i]=new TH2F(name.Data(),title.Data(),
00093                             60,-0.5,59.5,24,-0.5,23.5);
00094       h_edep_pske[i]=new TH3F(name2.Data(), title2.Data(),
00095                               60,-0.5,59.5,24,-0.5,23.5,70,-6,1);
00096    }
00097 
00098    if(dsave) dsave->cd();
00099 
00100 }

00231 {
00232 //======================================================================
00233 // Configure the module given the Registry r
00234 //======================================================================
00235 //  int    tmpi;
00236 
00237 //  if (r.Get("Mode",tmpi)) { fData = tmpi; }
00238 }

00209 {
00210 //======================================================================
00211 // Supply the default configuration for the module
00212 //======================================================================
00213   static Registry r; // Default configuration for module
00214 
00215   // Set name of config
00216   std::string name = this->GetName();
00217   name += ".config.default";
00218   r.SetName(name.c_str());
00219 
00220   // Set values in configuration
00221   r.UnLockValues();
00222   r.Set("Mode", 1);
00223   r.LockValues();
00224 
00225   return r;
00226 }

00105 {
00106 //======================================================================
00107 // FILL_IN: [Document your code!!]
00108 //======================================================================
00109    fOut->Write();
00110 
00111 }


Public Member Functions
	ScintHitTruthModule ()
	~ScintHitTruthModule ()
void	BeginJob ()
void	EndJob ()
JobCResult	Ana (const MomNavigator *mom)
const Registry &	DefaultConfig () const
void	Config (const Registry &r)
Private Attributes
TFile *	fOut
TH2 *	h_bvp
TH2 *	h_evp
TH2 *	h_evds
TH1 *	h_etot
TH2 *	h_dedsvpid
TH2 *	h_kevpid
TH2 *	h_dedsvke_e
TH2 *	h_dedsvke_mu
TH2 *	h_edep_ps [7]
TH3 *	h_edep_pske [7]

ScintHitTruthModule Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation