Majorana Simulation: io/src/MJOutputLANLClover.cc Source File

00001 //---------------------------------------------------------------------------//
00002 //bb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nu//
00003 //                                                                           //
00004 //                         MAJORANA Simulation                               //
00005 //                                                                           //
00006 //      This code implementation is the intellectual property of the         //
00007 //      MAJORANA Collaboration. It is based on Geant4, an intellectual       //
00008 //      property of the RD44 GEANT4 collaboration.                           //
00009 //                                                                           //
00010 //                        *********************                              //
00011 //                                                                           //
00012 //    Neither the authors of this software system, nor their employing       //
00013 //    institutes, nor the agencies providing financial support for this      //
00014 //    work  make  any representation or  warranty, express or implied,       //
00015 //    regarding this software system or assume any liability for its use.    //
00016 //    By copying, distributing or modifying the Program (or any work based   //
00017 //    on on the Program) you indicate your acceptance of this statement,     //
00018 //    and all its terms.                                                     //
00019 //                                                                           //
00020 //bb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nubb0nu//
00021 //---------------------------------------------------------------------------//
00022 //                                                          
00023 // $Id: MJOutputLANLClover.cc,v 1.2 2004/11/09 13:42:39 xliu Exp $ 
00024 //      
00025 // CLASS IMPLEMENTATION:  MJOutputLANLClover.cc
00026 //
00027 //---------------------------------------------------------------------------//
00033 // 
00034 //---------------------------------------------------------------------------//
00051 //---------------------------------------------------------------------------//
00052 //
00053 
00054 #include <CLHEP/Random/RandGauss.h>
00055 #include <CLHEP/Units/PhysicalConstants.h>
00056 
00057 #include "globals.hh"
00058 #include "G4PrimaryVertex.hh"
00059 #include "G4PrimaryParticle.hh"
00060 #include "G4Material.hh"
00061 #include "G4Step.hh"
00062 #include "G4Track.hh"
00063 #include "G4Vector3D.hh"
00064 
00065 #include "database/MJDatabase.hh"
00066 #include "database/MJDatabaseCloverCrystal.hh"
00067 #include "database/MJDatabaseCloverDetector.hh"
00068 #include "io/MJLogger.hh"
00069 #include "io/MJOutputRoot.hh"
00070 
00071 //---------------------------------------------------------------------------//
00072 
00073 #include "io/MJOutputLANLClover.hh"
00074 
00075 //---------------------------------------------------------------------------//
00076 //---------------------------------------------------------------------------//
00077 
00078 MJOutputLANLClover::MJOutputLANLClover(G4String serNum, G4bool isMother):
00079   MJOutputRoot(isMother), fSerialNumber(serNum), fTreePointer(0)
00080 {
00081   SetSchemaDefined(false);
00082   SetWaveformsSaved(false);
00083 }
00084 
00085 //---------------------------------------------------------------------------//
00086 
00087 MJOutputLANLClover::~MJOutputLANLClover()
00088 {
00089   delete fTreePointer;
00090 }
00091 
00092 //--------------------------------------------------------------------------//
00093 
00094 void MJOutputLANLClover::BeginOfEventAction(const G4Event *event)
00095 {
00096   MJOutputLANLCloverDataNoPS   *tP = fTreePointer;
00097 
00098   fE1_MC = fE2_MC = fE3_MC = fE4_MC = fPl_MC = fPm_MC = fPr_MC = 0;
00099 
00100   tP->T1 = tP->A1 = tP->M1 = tP->TriggerT1 = 
00101   tP->T2 = tP->A2 = tP->M2 = tP->TriggerT2 =
00102   tP->T3 = tP->A3 = tP->M3 = tP->TriggerT3 = 
00103   tP->T4 = tP->A4 = tP->M4 = tP->TriggerT4 = 
00104   tP->Tl = tP->Al = tP->Ml = tP->TriggerTl = 
00105   tP->Tm = tP->Am = tP->Mm = tP->TriggerTm =
00106   tP->Tr = tP->Ar = tP->Mr = tP->TriggerTr = 
00107   tP->ECsI = tP->TECsI = tP->AECsI = tP->MECsI = tP->TriggerTECsI = 0.0;
00108 
00109   tP->EventNumber = G4ToRoot(event->GetEventID());
00110 }
00111 
00112 //---------------------------------------------------------------------------//
00113 
00114 void MJOutputLANLClover::BeginOfRunAction()
00115 {
00116   // Get sensitive material.
00117   fSensitiveMaterialIndex = 
00118     G4Material::GetMaterial("Germanium-Nat")->GetIndex();
00119 
00120   // Get x coordinates of yz planes separating LM and MR segments.
00121   // This is required to determine which segment records energy.
00122   MJDatabaseCloverDetector *theDBDetector = 
00123     MJDatabase::GetCloverDetector(fSerialNumber);
00124   string *crystalSerialNumbers = theDBDetector->GetCrystalSerialNumbers();
00125   G4double spacerWidth = theDBDetector->GetSpacerWidth();
00126   MJDatabaseCloverCrystal *leftDBCrystal = 
00127     MJDatabase::GetCloverCrystal(crystalSerialNumbers[0]);
00128   MJDatabaseCloverCrystal *rightDBCrystal = 
00129     MJDatabase::GetCloverCrystal(crystalSerialNumbers[1]);
00130   fLMPlaneX = -(leftDBCrystal->GetCrystalRadius() - 
00131                 leftDBCrystal->GetRightDeficit() + spacerWidth);
00132   fMRPlaneX = rightDBCrystal->GetCrystalRadius() - 
00133               rightDBCrystal->GetLeftDeficit() + spacerWidth;
00134 
00135   // Load calibration data from database.
00136   // *** Not Implemented in dbase yet. Do by hand for now. ***
00137 
00138   fInverseCalibrationA[0] = 1.0 / 0.053209;
00139   fInverseCalibrationA[1] = 1.0 / 0.050614; 
00140   fInverseCalibrationA[2] = 1.0 / 0.051829;
00141   fInverseCalibrationA[3] = 1.0 / 0.051893;
00142   fInverseCalibrationA[4] = 1.0 / 0.057337;
00143   fInverseCalibrationA[5] = 1.0 / 0.064544;
00144   fInverseCalibrationA[6] = 1.0 / 0.06;  // Guesstimate
00145 
00146   fCalibrationB[0] = 0.24;
00147   fCalibrationB[1] = 0.23;
00148   fCalibrationB[2] = 0.21;
00149   fCalibrationB[3] = 0.31;
00150   fCalibrationB[4] = -0.86;
00151   fCalibrationB[5] = 1.1;
00152   fCalibrationB[6] = 0.0;  //Guesstimate
00153 
00154   fStandardDeviationA[0] = 1.0028;
00155   fStandardDeviationA[1] = 1.0615;
00156   fStandardDeviationA[2] = 1.0842;
00157   fStandardDeviationA[3] = 0.97019;
00158   fStandardDeviationA[4] = 3.9196;
00159   fStandardDeviationA[5] = 4.0077;
00160   fStandardDeviationA[6] = 4.0;   // Guesstimate
00161 
00162   fStandardDeviationB[0] = 0.063079;
00163   fStandardDeviationB[1] = 0.057654;
00164   fStandardDeviationB[2] = 0.057929;
00165   fStandardDeviationB[3] = 0.075027;
00166   fStandardDeviationB[4] = 0.10728;
00167   fStandardDeviationB[5] = 0.22727;
00168   fStandardDeviationB[6] = 0.2;  // Guesstimate
00169 
00170   // Convert FWHM to Standard Deviation for RandGauss generator.
00171   for(Int_t i = 0; i < 7; i++) {
00172     fStandardDeviationA[i] /= 2.355;
00173     fStandardDeviationB[i] /= 2.355;
00174   }
00175 
00176   // Root init.
00177   if(IsMother()) 
00178     OpenRootFile(GetFileName());
00179   DefineSchema();
00180 
00181 }
00182 
00183 //--------------------------------------------------------------------------//
00184 
00185 Float_t MJOutputLANLClover::ConvertEnergyToBin(Float_t edep, const Int_t ch)
00186 {
00187   Float_t channel;
00188   // Account for finite energy resolution by adding Gaussian fluctuations
00189   // of form sigma = a + b * sqrt(E). a, b from calibration data.
00190   edep += RandGauss::shoot(0.0, fStandardDeviationA[ch] + 
00191                            fStandardDeviationB[ch] * TMath::Sqrt(edep));
00192   
00193   // "Undo" calibration for each crystal and convert to 14 bit ADC channel
00194   channel = (edep - fCalibrationB[ch]) * fInverseCalibrationA[ch];
00195   return (channel < 0.0) ? 0.0 : channel;
00196 }
00197 
00198 //--------------------------------------------------------------------------//
00199 
00200 void MJOutputLANLClover::DefineSchema()
00201 {
00202 
00203   if(!SchemaDefined()) {
00204     TTree *nT;
00205 
00206     // Create output Tree if it has not been assigned.
00207     if(fTree)
00208       nT = fTree;
00209     else {
00210       if(!IsMother())
00211         MJLog(warning) << "No tree assigned to child output class." << endlog;
00212       nT = fTree = new TTree("fTree", "LANL Clover Tree");
00213     }
00214 
00215     // Create Tree schema from object.
00216     MJOutputLANLCloverDataNoPS   *tP;
00217     tP = fTreePointer = new MJOutputLANLCloverDataNoPS;
00218 
00219     // MC Header info.
00220     nT->Branch("EventNumber", &tP->EventNumber, "EventNumber/i");//Event Number
00221     nT->Branch("NumberOfParticles", &tP->NumberOfParticles,
00222                "NumberOfParticles/I");
00223     // Num. of particles generated at beginning of event.
00224     nT->Branch("Particle", tP->Particle, "Particle[NumberOfParticles]/I");
00225     // Enumerated particle type (PDG code). 
00226     nT->Branch("MomentumX", tP->MomentumX, "MomentumX[NumberOfParticles]/F");
00227     nT->Branch("MomentumY", tP->MomentumY, "MomentumY[NumberOfParticles]/F");
00228     nT->Branch("MomentumZ", tP->MomentumZ, "MomentumZ[NumberOfParticles]/F");
00229     // Initial momentum vector of each particle (keV/c).
00230     nT->Branch("PositionX", tP->PositionX, "PositionX[NumberOfParticles]/F");
00231     nT->Branch("PositionY", tP->PositionY, "PositionY[NumberOfParticles]/F");
00232     nT->Branch("PositionZ", tP->PositionZ, "PositionZ[NumberOfParticles]/F");
00233     // Initial position of each particle (cm).
00234 
00235     nT->Branch("E1_MC", &tP->E1_MC, "E1_MC/F"); // Energy deposit (MC truth)
00236     nT->Branch("E2_MC", &tP->E2_MC, "E2_MC/F"); // in keV
00237     nT->Branch("E3_MC", &tP->E3_MC, "E3_MC/F"); 
00238     nT->Branch("E4_MC", &tP->E4_MC, "E4_MC/F"); 
00239     nT->Branch("Pl_MC", &tP->Pl_MC, "Pl_MC/F"); 
00240     nT->Branch("Pm_MC", &tP->Pm_MC, "Pm_MC/F"); 
00241     nT->Branch("Pr_MC", &tP->Pr_MC, "Pr_MC/F"); 
00242 
00243   
00244     //Energy channel-1
00245     nT->Branch("E1",&tP->E1,"E1/F"); // Energy output, channels
00246     nT->Branch("T1",&tP->T1,"T1/F"); // Pulse Width
00247     nT->Branch("A1",&tP->A1,"A1/F"); // Pulse Asymmetry
00248     nT->Branch("M1",&tP->M1,"M1/F"); // Normalized moment
00249     nT->Branch("TriggerT1",&tP->TriggerT1,"TriggerT1/F"); // Trigger Time
00250 
00251     //Energy channel-2
00252     nT->Branch("E2",&tP->E2,"E2/F");
00253     nT->Branch("T2",&tP->T2,"T2/F");
00254     nT->Branch("A2",&tP->A2,"A2/F");
00255     nT->Branch("M2",&tP->M2,"M2/F");
00256     nT->Branch("TriggerT2",&tP->TriggerT2,"TriggerT2/F");
00257 
00258     //Energy channel-3
00259     nT->Branch("E3",&tP->E3,"E3/F");
00260     nT->Branch("T3",&tP->T3,"T3/F");
00261     nT->Branch("A3",&tP->A3,"A3/F");
00262     nT->Branch("M3",&tP->M3,"M3/F");
00263     nT->Branch("TriggerT3",&tP->TriggerT3,"TriggerT3/F");
00264 
00265     //Energy channel-4
00266     nT->Branch("E4",&tP->E4,"E4/F");
00267     nT->Branch("T4",&tP->T4,"T4/F");
00268     nT->Branch("A4",&tP->A4,"A4/F");
00269     nT->Branch("M4",&tP->M4,"M4/F");
00270     nT->Branch("TriggerT4",&tP->TriggerT4,"TriggerT4/F");
00271 
00272     //Position channel-1
00273     nT->Branch("Pl",&tP->Pl,"PL/F");
00274     nT->Branch("Tl",&tP->Tl,"Tl/F");
00275     nT->Branch("Al",&tP->Al,"Al/F");
00276     nT->Branch("Ml",&tP->Ml,"Ml/F");
00277     nT->Branch("TriggerTl",&tP->TriggerTl,"TriggerTl/F");
00278 
00279     //Position channel-2
00280     nT->Branch("Pm",&tP->Pm,"Pm/F");
00281     nT->Branch("Tm",&tP->Tm,"Tm/F");
00282     nT->Branch("Am",&tP->Am,"Am/F");
00283     nT->Branch("Mm",&tP->Mm,"Mm/F");
00284     nT->Branch("TriggerTm",&tP->TriggerTm,"TriggerTm/F");
00285 
00286     //Position channel-3
00287     nT->Branch("Pr",&tP->Pr,"Pr/F");
00288     nT->Branch("Tr",&tP->Tr,"Tr/F");
00289     nT->Branch("Ar",&tP->Ar,"Ar/F");
00290     nT->Branch("Mr",&tP->Mr,"Mr/F");
00291     nT->Branch("TriggerTr",&tP->TriggerTr,"TriggerTr/F");
00292 
00293     //ECsI detector
00294     nT->Branch("ECsI",&tP->ECsI,"ECsI/F");
00295     nT->Branch("TECsI",&tP->TECsI,"TECsI/F");
00296     nT->Branch("AECsI",&tP->AECsI,"AECsI/F");
00297     nT->Branch("MECsI",&tP->MECsI,"MECsI/F");
00298     nT->Branch("TriggerTECsI",&tP->TriggerTECsI,"TriggerTECsI/F");
00299 
00300     SetSchemaDefined(true);
00301   }
00302 }
00303 
00304 
00305 //---------------------------------------------------------------------------//
00306 
00307 void MJOutputLANLClover::EndOfEventAction(const G4Event *event)
00308 {
00309   // Save information about primary particle(s)
00310   G4PrimaryParticle *primaryParticle;
00311   G4PrimaryVertex   *primaryVertex;
00312   MJOutputLANLCloverDataNoPS  *tP = fTreePointer;
00313   G4int numberOfParticles = event->GetNumberOfPrimaryVertex();
00314 
00315   if(numberOfParticles > MAX_N_PARTICLES) {
00316     MJLog(warning) << numberOfParticles << "particles generated. Currently"
00317                    << " only " << MAX_N_PARTICLES << " are supported. "
00318                    << " EventID : " << event->GetEventID() << endlog;
00319     numberOfParticles = MAX_N_PARTICLES;
00320   }
00321 
00322   for(G4int i = 0; i < numberOfParticles; i++) {
00323     primaryVertex = event->GetPrimaryVertex(i);
00324     primaryParticle = primaryVertex->GetPrimary();
00325     tP->Particle[i] = G4ToRoot(primaryParticle->GetPDGcode());
00326     tP->MomentumX[i] = 
00327       G4ToRoot(primaryParticle->GetPx() / keV);
00328     tP->MomentumY[i] = 
00329       G4ToRoot(primaryParticle->GetPy() / keV);
00330     tP->MomentumZ[i] = 
00331       G4ToRoot(primaryParticle->GetPz() / keV);
00332     tP->PositionX[i] = G4ToRoot(primaryVertex->GetX0()/cm);
00333     tP->PositionY[i] = G4ToRoot(primaryVertex->GetY0()/cm);
00334     tP->PositionZ[i] = G4ToRoot(primaryVertex->GetZ0()/cm);
00335   }
00336 
00337   tP->NumberOfParticles = numberOfParticles;
00338 
00339   // Convert and save accumulated energies in keV.
00340   
00341   tP->E1_MC = (Float_t) G4ToRoot(fE1_MC / keV);
00342   tP->E2_MC = (Float_t) G4ToRoot(fE2_MC / keV);
00343   tP->E3_MC = (Float_t) G4ToRoot(fE3_MC / keV);
00344   tP->E4_MC = (Float_t) G4ToRoot(fE4_MC / keV);
00345   tP->Pl_MC = (Float_t) G4ToRoot(fPl_MC / keV);
00346   tP->Pm_MC = (Float_t) G4ToRoot(fPr_MC / keV);
00347   tP->Pr_MC = (Float_t) G4ToRoot(fPm_MC / keV);
00348   
00349   tP->E1 = ConvertEnergyToBin(tP->E1_MC, 0);
00350   tP->E2 = ConvertEnergyToBin(tP->E2_MC, 1);
00351   tP->E3 = ConvertEnergyToBin(tP->E3_MC, 2);
00352   tP->E4 = ConvertEnergyToBin(tP->E4_MC, 3);
00353   tP->Pl = ConvertEnergyToBin(tP->Pl_MC, 4);
00354   tP->Pm = ConvertEnergyToBin(tP->Pm_MC, 5);
00355   tP->Pr = ConvertEnergyToBin(tP->Pr_MC, 6);
00356  
00357   if(IsMother())
00358     FillTree();
00359 }
00360 
00361 //---------------------------------------------------------------------------//
00362 
00363 void MJOutputLANLClover::EndOfRunAction()
00364 {
00365   if(IsMother())
00366     CloseRootFile();
00367 }
00368 
00369 //---------------------------------------------------------------------------//
00370 
00371 void MJOutputLANLClover::RootSteppingAction(const G4Step* step)
00372 {
00373   G4Track *track = step->GetTrack();
00374   G4Material *material = track->GetMaterial();
00375   
00376   // Determine in which crystal and segment the step occurred. 
00377   if(material->GetIndex() == fSensitiveMaterialIndex){
00378     G4ThreeVector position = track->GetPosition();
00379     G4double x = position.x();
00380     G4double deltaE = step->GetTotalEnergyDeposit();
00381     if(x > 0.0) {
00382       if(position.y() > 0.0){
00383         fE2_MC += deltaE;
00384         if(x > fMRPlaneX){
00385           fPr_MC += deltaE;
00386         } else {
00387           fPm_MC += deltaE;
00388         }
00389       } else {
00390         fE3_MC += deltaE;
00391         if(x > fMRPlaneX){
00392           fPr_MC += deltaE;
00393         } else {
00394           fPm_MC += deltaE;
00395         }       
00396       }
00397     } else {
00398       if(position.y() > 0.0){
00399         fE1_MC += deltaE;
00400         if(x < fLMPlaneX){
00401           fPl_MC += deltaE;
00402         } else {
00403           fPm_MC += deltaE;
00404         }
00405       } else {
00406         fE4_MC += deltaE;
00407         if(x < fLMPlaneX){
00408           fPl_MC += deltaE;
00409         } else {
00410           fPm_MC += deltaE;
00411         }
00412       }
00413     }
00414   }
00415 }
00416 
00417 //---------------------------------------------------------------------------//
00418 //---------------------------------------------------------------------------//
00419 
00420 
00421 
00422 
00423 
00424