AlcapDAQ: rootana/src/Al50/TME_Al50

00001 #include "TME_Al50_EvdE.h"
00002 #include "RegisterModule.inc"
00003 #include "TGlobalData.h"
00004 #include "TSetupData.h"
00005 #include "ModulesOptions.h"
00006 #include "definitions.h"
00007 #include "TMuonEvent.h"
00008 #include "debug_tools.h"
00009 
00010 #include <TH1F.h>
00011 #include <TH2F.h>
00012 #include <TCanvas.h>
00013 #include <TApplication.h>
00014 
00015 #include <iostream>
00016 #include <cmath>
00017 using std::cout;
00018 using std::endl;
00019 
00020 extern MuonEventList gMuonEvents;
00021 
00022 TME_Al50_EvdE::TME_Al50_EvdE(modules::options* opts):
00023   BaseModule("TME_Al50_EvdE",opts),fNullCount(0),fTdpCount(0){
00024 
00025   fStoppedProtonCut=opts->GetBool("stopped_proton_cut", false);
00026 }
00027 
00028 TME_Al50_EvdE::~TME_Al50_EvdE(){
00029 }
00030 
00031 int TME_Al50_EvdE::BeforeFirstEntry(TGlobalData* gData,const TSetupData *setup){
00032 
00033   fLeftArm.detname = "SiL";
00034   fLeftArmQuad1.detname = "SiL_Quad1";
00035   fLeftArmQuad2.detname = "SiL_Quad2";
00036   fLeftArmQuad3.detname = "SiL_Quad3";
00037   fLeftArmQuad4.detname = "SiL_Quad4";
00038   fRightArm.detname = "SiR";
00039   fRightArmQuad1.detname = "SiR_Quad1";
00040   fRightArmQuad2.detname = "SiR_Quad2";
00041   fRightArmQuad3.detname = "SiR_Quad3";
00042   fRightArmQuad4.detname = "SiR_Quad4";
00043 
00044 
00045   using namespace IDs;
00046   fSiL1.push_back(IDs::channel (kSiL1_1 , kNotApplicable ));
00047   fSiL1.push_back(IDs::channel (kSiL1_2 , kNotApplicable ));
00048   fSiL1.push_back(IDs::channel (kSiL1_3 , kNotApplicable ));
00049   fSiL1.push_back(IDs::channel (kSiL1_4 , kNotApplicable ));
00050   fSiL2 = new IDs::channel (kSiL2   , kNotApplicable );
00051   fSiR1.push_back(IDs::channel (kSiR1_1 , kNotApplicable ));
00052   fSiR1.push_back(IDs::channel (kSiR1_2 , kNotApplicable ));
00053   fSiR1.push_back(IDs::channel (kSiR1_3 , kNotApplicable ));
00054   fSiR1.push_back(IDs::channel (kSiR1_4 , kNotApplicable ));
00055   fSiR2 = new IDs::channel (kSiR2   , kNotApplicable );
00056 
00057   fLeftArm.thin = fSiL1;
00058   fLeftArm.thick = fSiL2;
00059   DetectorList left_quad_1; left_quad_1.push_back(IDs::channel (kSiL1_1 , kNotApplicable));
00060   fLeftArmQuad1.thin = left_quad_1;
00061   fLeftArmQuad1.thick = fSiL2;
00062   DetectorList left_quad_2; left_quad_2.push_back(IDs::channel (kSiL1_2 , kNotApplicable));
00063   fLeftArmQuad2.thin = left_quad_2;
00064   fLeftArmQuad2.thick = fSiL2;
00065   DetectorList left_quad_3; left_quad_3.push_back(IDs::channel (kSiL1_3 , kNotApplicable));
00066   fLeftArmQuad3.thin = left_quad_3;
00067   fLeftArmQuad3.thick = fSiL2;
00068   DetectorList left_quad_4; left_quad_4.push_back(IDs::channel (kSiL1_4 , kNotApplicable));
00069   fLeftArmQuad4.thin = left_quad_4;
00070   fLeftArmQuad4.thick = fSiL2;
00071 
00072   fRightArm.thin = fSiR1;
00073   fRightArm.thick = fSiR2;
00074   DetectorList right_quad_1; right_quad_1.push_back(IDs::channel (kSiR1_1 , kNotApplicable));
00075   fRightArmQuad1.thin = right_quad_1;
00076   fRightArmQuad1.thick = fSiR2;
00077   DetectorList right_quad_2; right_quad_2.push_back(IDs::channel (kSiR1_2 , kNotApplicable));
00078   fRightArmQuad2.thin = right_quad_2;
00079   fRightArmQuad2.thick = fSiR2;
00080   DetectorList right_quad_3; right_quad_3.push_back(IDs::channel (kSiR1_3 , kNotApplicable));
00081   fRightArmQuad3.thin = right_quad_3;
00082   fRightArmQuad3.thick = fSiR2;
00083   DetectorList right_quad_4; right_quad_4.push_back(IDs::channel (kSiR1_4 , kNotApplicable));
00084   fRightArmQuad4.thin = right_quad_4;
00085   fRightArmQuad4.thick = fSiR2;
00086 
00087 
00088   // Hard-coded for the time being
00089   fLeftArm.lower_time_cut = 0;
00090   fLeftArm.upper_time_cut = 10000;
00091   fLeftArmQuad1.lower_time_cut = 0;
00092   fLeftArmQuad1.upper_time_cut = 999999999;
00093   fLeftArmQuad2.lower_time_cut = 0;
00094   fLeftArmQuad2.upper_time_cut = 999999999;
00095   fLeftArmQuad3.lower_time_cut = 0;
00096   fLeftArmQuad3.upper_time_cut = 999999999;
00097   fLeftArmQuad4.lower_time_cut = 0;
00098   fLeftArmQuad4.upper_time_cut = 999999999;
00099 
00100   fRightArm.lower_time_cut = 0;
00101   fRightArm.upper_time_cut = 10000;
00102   fRightArmQuad1.lower_time_cut = 0;
00103   fRightArmQuad1.upper_time_cut = 999999999;
00104   fRightArmQuad2.lower_time_cut = 0;
00105   fRightArmQuad2.upper_time_cut = 999999999;
00106   fRightArmQuad3.lower_time_cut = 0;
00107   fRightArmQuad3.upper_time_cut = 999999999;
00108   fRightArmQuad4.lower_time_cut = 0;
00109   fRightArmQuad4.upper_time_cut = 999999999;
00110 
00111   fArms.push_back(fLeftArm);
00112   fArms.push_back(fLeftArmQuad1);
00113   fArms.push_back(fLeftArmQuad2);
00114   fArms.push_back(fLeftArmQuad3);
00115   fArms.push_back(fLeftArmQuad4);
00116   fArms.push_back(fRightArm);
00117   fArms.push_back(fRightArmQuad1);
00118   fArms.push_back(fRightArmQuad2);
00119   fArms.push_back(fRightArmQuad3);
00120   fArms.push_back(fRightArmQuad4);
00121 
00122   // Create the histograms for each arm
00123   for (std::vector<Arm>::iterator i_arm = fArms.begin(); i_arm != fArms.end(); ++i_arm) {
00124     std::string evde_histname = i_arm->detname + "_EvdE";
00125     std::string evde_histtitle = "The E v dE plot for " + i_arm->detname;
00126     i_arm->h_EvdE = new TH2F(evde_histname.c_str(), evde_histtitle.c_str(), 250,0,25000, 100,0,10000);
00127     i_arm->h_EvdE->SetXTitle("E + dE [keV]");
00128     i_arm->h_EvdE->SetYTitle("dE [keV]");
00129 
00130     std::string time_histname = i_arm->detname + "_Time";
00131     std::string time_histtitle = "The time of hits in " + i_arm->detname;
00132     i_arm->h_Time = new TH1F(time_histname.c_str(), time_histtitle.c_str(), 2500,0,10000);
00133     i_arm->h_Time->SetXTitle("Time [ns]");
00134     i_arm->h_Time->SetYTitle("Count");    
00135 
00136     if (fStoppedProtonCut) {
00137       TTree* pid_cuts_tree = new TTree();
00138       std::string pid_cuts_filename = "/gpfs/home/edmonds_a/AlcapDAQ/analyzer/rootana/src/Al50/pid-cuts-" + i_arm->detname + ".txt";
00139       pid_cuts_tree->ReadFile(pid_cuts_filename.c_str());
00140       pid_cuts_tree->Print();
00141 
00142       TBranch* energy_branch = (TBranch*) pid_cuts_tree->GetBranch("energy");
00143       double energy;
00144       energy_branch->SetAddress(&energy);
00145 
00146       TBranch* d_energy_branch = (TBranch*) pid_cuts_tree->GetBranch("dEnergy");
00147       double d_energy;
00148       d_energy_branch->SetAddress(&d_energy);
00149 
00150       TBranch* stopped_proton_prob_branch = (TBranch*) pid_cuts_tree->GetBranch("p_stop_prob");
00151       double stopped_proton_prob;
00152       stopped_proton_prob_branch->SetAddress(&stopped_proton_prob);
00153 
00154       std::string prob_histname = i_arm->detname + "_StoppedProtonProb";
00155       std::string prob_histtitle = "Probability that at a given (E, dE) that the hit is a proton in " + i_arm->detname;
00156       i_arm->h_stopped_proton_prob = new TH2F(prob_histname.c_str(), prob_histtitle.c_str(), 100,0,10000, 100,0,10000);
00157       i_arm->h_stopped_proton_prob->SetXTitle("E [keV]");
00158       i_arm->h_stopped_proton_prob->SetYTitle("dE [keV]");
00159 
00160       for (int i_entry = 0; i_entry < pid_cuts_tree->GetEntries(); ++i_entry) {
00161         pid_cuts_tree->GetEntry(i_entry);
00162         i_arm->h_stopped_proton_prob->Fill(energy, d_energy, stopped_proton_prob);
00163         //      std::cout << i_arm->detname << ": Filling @ (" << energy << ", " << d_energy << "): " << stopped_proton_prob << std::endl;
00164       }
00165       delete pid_cuts_tree;
00166     }
00167     else {
00168       i_arm->h_stopped_proton_prob = NULL;
00169     }
00170 
00171   }
00172 
00173   return 0;
00174 }
00175 
00176 int TME_Al50_EvdE::ProcessEntry(TGlobalData* gData,const TSetupData *setup){
00177 
00178     for(MuonEventList::const_iterator i_tme=gMuonEvents.begin();
00179             i_tme!=gMuonEvents.end(); ++i_tme){
00180 
00181       // First, check for pile-up in the muSc
00182       if ( (*i_tme)->HasMuonPileup()) {
00183         continue;
00184       }
00185 
00186       double tme_time= (*i_tme)->GetTime(); // this is the same as the muSc time
00187 
00188       // Loop through the arms and plots
00189       for (std::vector<Arm>::const_iterator i_arm = fArms.begin(); i_arm != fArms.end(); ++i_arm) {
00190 
00191         // Check that there are pulses in both the thin and the thick silicon detectors
00192         DetectorList si_thin = (*i_arm).thin;
00193         IDs::channel* si_thick = (*i_arm).thick;
00194 
00195         // Get the thick detector
00196         int si_thick_source_index=(*i_tme)->GetSourceIndex(*si_thick);
00197 
00198         // While there are thick detector sources still around
00199         while (si_thick_source_index>-1) {
00200           const IDs::source& si_thick_source=(*i_tme)->GetSource(si_thick_source_index);
00201           int n_si_thick = (*i_tme)->NumPulses(si_thick_source);
00202 
00203           // Loop through the pulses in the thick detector
00204           for (int i_thick_pulse=0; i_thick_pulse<n_si_thick; ++i_thick_pulse) {
00205             const TDetectorPulse* tdp_si_thick=(*i_tme)->GetPulse(si_thick_source,i_thick_pulse);
00206 
00207             double thick_energy = tdp_si_thick->GetTAP(TDetectorPulse::kFast)->GetEnergy();
00208             // Skip to next pulse if energy <100 keV
00209             if (thick_energy < 100) {
00210               continue;
00211             }
00212 
00213             double thick_amplitude = tdp_si_thick->GetAmplitude();
00214             double thick_time = tdp_si_thick->GetTime();
00215 
00216             // Loop through the thin quadrants
00217             for(DetectorList::const_iterator i_det=si_thin.begin();
00218                 i_det!=si_thin.end(); ++i_det){
00219 
00220               // Get the thin quadrant and loop through its sources
00221               int si_thin_source_index=(*i_tme)->GetSourceIndex(*i_det);
00222               while(si_thin_source_index>-1){
00223                 const IDs::source& si_thin_source=(*i_tme)->GetSource(si_thin_source_index);
00224                 int n_si_thin = (*i_tme)->NumPulses(si_thin_source);
00225 
00226                 // Loop through the pulses in this quadrant
00227                 for(int i_thin_pulse=0; i_thin_pulse<n_si_thin; ++i_thin_pulse){
00228                   const TDetectorPulse* tdp_si_thin=(*i_tme)->GetPulse(si_thin_source,i_thin_pulse);
00229                   double thin_energy = tdp_si_thin->GetTAP(TDetectorPulse::kFast)->GetEnergy();
00230 
00231                   if (thin_energy < 100) {
00232                     continue;
00233                   }
00234                   double thin_amplitude = tdp_si_thin->GetAmplitude();
00235                   double thin_time = tdp_si_thin->GetTime();
00236                   
00237                   bool passes_cuts = false;
00238 
00239                   double arrival_time = thick_time - tme_time;
00240                   if ( arrival_time > i_arm->lower_time_cut && arrival_time < i_arm->upper_time_cut ) { 
00241                     // Now check if this passes our proton cut
00242                     if (fStoppedProtonCut) {
00243                       int bin = i_arm->h_stopped_proton_prob->FindBin(thick_energy+thin_energy, thin_energy);
00244                       double probability = i_arm->h_stopped_proton_prob->GetBinContent(bin);
00245                       //                      std::cout << "(E+dE, dE) = (" << thick_energy+thin_energy << ", " << thin_energy << "): Prob = " << probability << std::endl;
00246                       if (probability > 0.99) {
00247                         passes_cuts = true;
00248                       }
00249                     }
00250                     else { // everything passes if we're not cutting on protons
00251                       passes_cuts = true;
00252                     }
00253                   }
00254 
00255                   if (std::abs(thin_time - thick_time) > 500) {
00256                     //              std::cout << "abs(t_thin - t_thick) = " << std::fabs(thin_time - thick_time) << std::endl;
00257                     passes_cuts=false;
00258                   }
00259                   if (fStoppedProtonCut && (thick_energy+thin_energy < 1500 || thick_energy+thin_energy > 8000)) {
00260                     passes_cuts=false; // only accept proton if E is between 1.5 and 8 MeV (like Nam's cut)
00261                   }
00262                   if (passes_cuts) {
00263                     //              std::cout << "Amplitude --> Energy (thick): " << thick_amplitude << " --> " << thick_energy << std::endl;
00264                     //              std::cout << "Amplitude --> Energy (thin): " << thin_amplitude << " --> " << thin_energy << std::endl;
00265                     //              std::cout << "Plotting: " << thick_energy+thin_energy << ", " << thin_energy << std::endl;
00266                     i_arm->h_EvdE->Fill(thick_energy+thin_energy, thin_energy);
00267                     i_arm->h_Time->Fill(arrival_time);
00268                   }
00269 
00270                 }
00271                 si_thin_source_index=(*i_tme)->GetSourceIndex(*i_det,si_thin_source_index+1);
00272               }
00273             }
00274           }
00275           si_thick_source_index=(*i_tme)->GetSourceIndex(*si_thick,si_thick_source_index+1);
00276         }
00277       }
00278     }
00279     
00280     return 0;
00281 }
00282 
00283 int TME_Al50_EvdE::AfterLastEntry(TGlobalData* gData,const TSetupData *setup){
00284 
00285   return 0;
00286 }
00287 
00288 ALCAP_REGISTER_MODULE(TME_Al50_EvdE, stopped_proton_cut);
rootana/src/Al50/TME_Al50_EvdE.cpp
