// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // $Id: HistoManager.cc,v 1.1 2010-09-21 17:55:58 maire Exp $ // GEANT4 tag $Name: not supported by cvs2svn $ // //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... #include "HistoManager.hh" #include "HistoMessenger.hh" #include "G4UnitsTable.hh" #include "G4SystemOfUnits.hh" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" #include "G4Threading.hh" #include "G4UIcommand.hh" // line to check svn //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //G4ThreadLocal HistoManager* HistoManager::fHistoManager = 0; HistoManager* HistoManager::GetHistoManager() { if(fHistoManager) return fHistoManager; else{ fHistoManager = new HistoManager(); return fHistoManager;} } HistoManager::HistoManager() :factoryOn(false) { fileName[0] = "iThembaPassiveBeamline"; fileType = "root"; fileOption = "export=root"; histoMessenger = new HistoMessenger(this); //threadID = G4Threading::G4GetThreadId(); if(fHistoManager) { G4Exception("HistoManager::HistoManager", "Run9999", FatalException, "G4HistoManager constructed twice."); } fHistoManager = this; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... HistoManager::~HistoManager() { //delete histoMessenger; //delete fHistoManager; //delete file; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::book() { std::cout << "BOOKY BOOKY" << std::endl; // Creating a tree mapped to an hbook file. // on a master thread threadID=-1 and so we dont add threadID G4int threadID = G4Threading::G4GetThreadId(); if(threadID<0) fileName[1] = fileName[0] + "." + fileType; else fileName[1] = fileName[0] + "_threadID"+ G4UIcommand::ConvertToString(threadID) + "." + fileType; file = new TFile(fileName[1], "RECREATE"); // create selected histograms for (G4int k=0; k Branch( "rndmSeed", &treeInfo[k].rndmSeed); treeInfo[k].tree -> Branch( "evtNumber" , &treeInfo[k].evtNumber ); treeInfo[k].tree -> Branch( "x", &treeInfo[k].x); treeInfo[k].tree -> Branch( "y", &treeInfo[k].y); treeInfo[k].tree -> Branch( "z", &treeInfo[k].z); treeInfo[k].tree -> Branch( "px", &treeInfo[k].px); treeInfo[k].tree -> Branch( "py", &treeInfo[k].py); treeInfo[k].tree -> Branch( "pz", &treeInfo[k].pz); treeInfo[k].tree -> Branch( "trackID", &treeInfo[k].trackID); treeInfo[k].tree -> Branch( "parentID", &treeInfo[k].parentID); treeInfo[k].tree -> Branch( "isIncident", &treeInfo[k].isIncident); treeInfo[k].tree -> Branch( "isExit", &treeInfo[k].isExit); treeInfo[k].tree -> Branch( "kineticEnergy", &treeInfo[k].kineticEnergy); treeInfo[k].tree -> Branch( "edep", &treeInfo[k].edep ); treeInfo[k].tree -> Branch( "niel", &treeInfo[k].niel); treeInfo[k].tree -> Branch( "pdgid", &treeInfo[k].pdgid); treeInfo[k].tree -> Branch( "charge", &treeInfo[k].charge); treeInfo[k].tree -> Branch( "processSubType", &treeInfo[k].processType); factoryOn = true; } } //if (factoryOn) // G4cout << "\n----> Histogram Tree is opened in " << file->GetName() << G4endl; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... #include void HistoManager::save() { // G4int threadID = G4Threading::G4GetThreadId(); // if(threadID<0) // fileName[1] = fileName[0] + "." + fileType; // else // fileName[1] = fileName[0] + "_threadID"+ G4UIcommand::ConvertToString(threadID) + "." + fileType; // file = new TFile(fileName[1], "RECREATE"); if (factoryOn && file->IsOpen()) { // print out the RT to a txt file format for ME // call here incase do not reach 1000 before end of run HistoManager::StoreRTInfo(); binout.close(); //saveAscii(); // Write ascii file, if any file->cd(); /*for (G4int k=0; kWrite(); } } // create selected 2D histograms for (G4int k=0; kWrite(); } } // create selectred TProfiles for(G4int k=0; kWrite(); } } for(G4int k=0; kWrite(); } } */ rtTree->Write(); trackerTree->Write(); stripRTTree->Write(); file->ls(); file->Write(); G4cout << "\n----> Histogram Tree is saved in " << file->GetName() << G4endl; file->Close(); //delete file; //factoryOn = false; } } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::FillHisto(const G4String& label, G4double e, G4double weight) { G4int ih = GetHistoAddress(label);//+"_threadID"+G4UIcommand::ConvertToString(G4Threading::G4GetThreadId())); if (ih < 0) { G4cout << "---> warning from HistoManager::FillHisto() : histo " << ih << "does not exist; e= " << e << " w= " << weight << G4endl; return; } if(histoInfo[ih].exist) histoInfo[ih].histo->Fill(e/histoInfo[ih].Unit, weight); } void HistoManager::FillHisto2D(const G4String& label, G4double x, G4double y, G4double weight) { G4int ih = GetHisto2DAddress(label); if (ih < 0) { G4cout << "---> warning from HistoManager::FillHisto2D() : histo2D " << ih << "does not exist; x= " << x << " y= " << y << " w= " << weight << G4endl; return; } if(histo2DInfo[ih].exist) histo2DInfo[ih].histo->Fill(x/histo2DInfo[ih].UnitX, y/histo2DInfo[ih].UnitY, weight); } void HistoManager::FillProfile(const G4String& label, G4double e, G4double weight) { G4int ih = GetProfileAddress(label); if (ih < 0) { G4cout << "---> warning from HistoManager::FillProfile() : profile " << ih << "does not exist; e= " << e << " w= " << weight << G4endl; return; } if(profileInfo[ih].exist) { profileInfo[ih].profile->Fill(e/profileInfo[ih].Unit, weight); } } void HistoManager::FillTree(const G4String& label, G4double rndmSeed, G4double evtNumber, G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, G4double trackID, G4double parentID, bool isIncident, bool isExit, G4double kineticEnergy, G4double edep, G4double niel, G4double pdgid, G4double charge, G4int processSubType) { G4int ih = GetTreeAddress(label); if (ih < 0) { G4cout << "---> warning from HistoManager::FillTree() : tree " << ih << "does not exist;" << G4endl; return; } if(treeInfo[ih].exist) { treeInfo[ih].rndmSeed = rndmSeed; treeInfo[ih].evtNumber = evtNumber; treeInfo[ih].x = x; treeInfo[ih].y = y; treeInfo[ih].z = z; treeInfo[ih].px = px; treeInfo[ih].py = py; treeInfo[ih].pz = pz; treeInfo[ih].trackID = trackID; treeInfo[ih].parentID = parentID; treeInfo[ih].isIncident = isIncident; treeInfo[ih].isExit = isExit; treeInfo[ih].kineticEnergy = kineticEnergy; treeInfo[ih].edep = edep; treeInfo[ih].niel = niel; treeInfo[ih].pdgid = pdgid; treeInfo[ih].charge = charge; treeInfo[ih].processType = double(processSubType); } treeInfo[ih].tree->Fill(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::SetHisto(const G4String& label, const G4String& title, G4int nbins, G4double valmin, G4double valmax, const G4String& unit) { G4String lab = label;// + "_threadID" + G4UIcommand::ConvertToString(G4Threading::G4GetThreadId()); G4cout << lab << G4endl; for(unsigned i=0; i warning from HistoManager::SetHisto() : histo " << histoInfo[i].Label << "already exists" << G4endl; return; } } HistoInfo hi; hi.Label = lab; histoInfo.push_back(hi); G4int ih = GetHistoAddress(lab); G4String titl = title; G4double vmin = valmin, vmax = valmax; histoInfo[ih].Unit = 1.; if (unit != "none") { titl = title + " (" + unit + ")"; histoInfo[ih].Unit = G4UnitDefinition::GetValueOf(unit); vmin = valmin/histoInfo[ih].Unit; vmax = valmax/histoInfo[ih].Unit; } histoInfo[ih].exist = true; histoInfo[ih].Label = lab; histoInfo[ih].Title = titl; histoInfo[ih].Nbins = nbins; histoInfo[ih].Vmin = vmin; histoInfo[ih].Vmax = vmax; histoInfo[ih].Width = (valmax-valmin)/nbins; G4cout << "----> SetHisto " << ih << ": " << lab << "; " << titl << "; " << nbins << " bins from " << vmin << " " << unit << " to " << vmax << " " << unit << G4endl; } void HistoManager::SetProfile(const G4String& label, const G4String& title, G4int nbins, G4double valmin, G4double valmax, const G4String& unit) { for(unsigned i=0; i warning from HistoManager::SetProfile() : profile " << profileInfo[i].Label << "already exists" << G4endl; return; } } ProfileInfo pi; pi.Label = label; G4String titl = title; G4double vmin = valmin, vmax = valmax; pi.Unit = 1.; if (unit != "none") { titl = title + " (" + unit + ")"; pi.Unit = G4UnitDefinition::GetValueOf(unit); vmin = valmin/pi.Unit; vmax = valmax/pi.Unit; } pi.exist = true; pi.Label = label; pi.Title = titl; pi.Nbins = nbins; pi.Vmin = vmin; pi.Vmax = vmax; pi.Width = (valmax-valmin)/nbins; profileInfo.push_back(pi); G4int ih = GetProfileAddress(label); G4cout << "----> SetProfile " << ih << ": " << label << "; "<< titl << "; " << nbins << " bins from " << vmin << " " << unit << " to " << vmax << " " << unit << G4endl; } void HistoManager::SetHisto2D(const G4String& label, const G4String& title, G4int nbinsx, G4double vminx, G4double vmaxx, G4int nbinsy, G4double vminy, G4double vmaxy, const G4String& unitx, const G4String& unity) { for(unsigned i=0; i warning from HistoManager::SetHisto() : histo " << histoInfo[i].Label << "already exists" << G4endl; return; } } Histo2DInfo hi; hi.Label = label; histo2DInfo.push_back(hi); G4int ih = GetHisto2DAddress(label); G4String titl = title; histo2DInfo[ih].UnitX = 1.; if (unitx != "none") { histo2DInfo[ih].UnitX = G4UnitDefinition::GetValueOf(unitx); vminx = vminx/histo2DInfo[ih].UnitX; vmaxx = vmaxx/histo2DInfo[ih].UnitX; } histo2DInfo[ih].UnitY = 1.; if (unity != "none") { histo2DInfo[ih].UnitY = G4UnitDefinition::GetValueOf(unity); vminy = vminy/histo2DInfo[ih].UnitY; vmaxy = vmaxy/histo2DInfo[ih].UnitY; } histo2DInfo[ih].exist = true; histo2DInfo[ih].Label = label; histo2DInfo[ih].Title = titl; histo2DInfo[ih].NbinsX = nbinsx; histo2DInfo[ih].NbinsY = nbinsy; histo2DInfo[ih].Xmin = vminx; histo2DInfo[ih].Xmax = vmaxx; histo2DInfo[ih].Ymin = vminy; histo2DInfo[ih].Ymax = vmaxy; G4cout << "----> SetHisto2D " << ih << ": " << label << "; " << titl << "; " << nbinsx << " x-bins from " << vminx << " " << unitx << " to " << vmaxx << " " << unitx << "; " << nbinsy << " y-bins from " << vminy << " " << unity << " to " << vmaxy << " " << unity << G4endl; } void HistoManager::SetTree(const G4String& label, const G4String& title) { for(unsigned i=0; i warning from HistoManager::SetTree() : tree " << treeInfo[i].Label << "already exists" << G4endl; return; } } TreeInfo ti; ti.Label = label; ti.Title = title; ti.exist = true; treeInfo.push_back(ti); } G4int HistoManager::GetHistoAddress(const G4String& label) { G4int address = -1; for(unsigned i=0; i=0); return address; } G4int HistoManager::GetHisto2DAddress(const G4String& label) { G4int address = -1; for(unsigned i=0; i=0); return address; } G4int HistoManager::GetProfileAddress(const G4String& label) { G4int address = -1; for(unsigned i=0; i=0); return address; } G4int HistoManager::GetTreeAddress(const G4String& label) { G4int address = -1; for(unsigned i=0; i=0); return address; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::Scale(const G4String& label, G4double fac) { G4int ih = GetHistoAddress(label); if (ih > MaxHisto) { G4cout << "---> warning from HistoManager::Scale() : histo " << ih << "does not exist (fac = " << fac << ")" << G4endl; return; } //if(histoInfo[ih].exist) histoInfo[ih].histo->scale(fac); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::RemoveHisto(const G4String& label) { G4int ih = GetHistoAddress(label); if (ih > MaxHisto) { G4cout << "---> warning from HistoManager::RemoveHisto() : histo " << ih << "does not exist" << G4endl; return; } histoInfo[ih].histo = 0; histoInfo[ih].exist = false; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::PrintHisto(const G4String& label) { G4int ih = GetHistoAddress(label); if (ih < MaxHisto) { histoInfo[ih].ascii = true; histoInfo[0].ascii = true; } else G4cout << "---> warning from HistoManager::PrintHisto() : histo " << ih << "does not exist" << G4endl; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... #include void HistoManager::saveAscii() { /* if (!ascii[0]) return; G4String name = fileName[0] + ".ascii"; std::ofstream File(name, std::ios::out); File.setf( std::ios::scientific, std::ios::floatfield ); //write selected histograms for (G4int ih=0; ihaxis().binLowerEdge(iBin) + histo[ih]->axis().binUpperEdge(iBin)) << "\t" << histo[ih]->binHeight(iBin) << G4endl; } } }*/ } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void HistoManager::AddBranch( TTree* input_tree, G4String title, std::vector data ) { if( (input_tree->GetEntries() != (int)data.size()) && (input_tree->GetEntries() != 0) ) { G4cout << "ERROR: CANNOT ADD BRANCH " << title << " DATA SIZE DID NOT MATCH TREE." << G4endl; return; } unsigned int size = (unsigned)data.size();//input_tree->GetEntries(); double thisValue=0.; TString name=title.c_str(); name.Append("/D"); TBranch* output_branch_obj = input_tree->Branch( title.c_str(), &thisValue, name ); input_tree->SetEntries( (int)size ); for( unsigned int i=0; i< size; ++i ) { thisValue = data[i]; output_branch_obj->Fill(); } } void HistoManager::AddBranch( TTree* input_tree, G4String title, std::vector data ) { if( (input_tree->GetEntries() != (int)data.size()) && (input_tree->GetEntries() != 0) ) { G4cout << "ERROR: CANNOT ADD BRANCH " << title << " DATA SIZE DID NOT MATCH TREE." << G4endl; return; } unsigned int size = (unsigned)data.size();//input_tree->GetEntries(); TString thisValue=""; TString name=title.c_str(); name.Append("/C"); TBranch* output_branch_obj = input_tree->Branch( title.c_str(), &thisValue, name ); input_tree->SetEntries( (int)size ); for( unsigned int i=0; i< size; ++i ) { thisValue = data[i] ; output_branch_obj->Fill(); } } void HistoManager::AddBranch( TTree* input_tree, G4String title, std::vector data ) { if( (input_tree->GetEntries() != (int)data.size()) && (input_tree->GetEntries() != 0) ) { G4cout << "ERROR: CANNOT ADD BRANCH " << title << " DATA SIZE DID NOT MATCH TREE." << G4endl; return; } unsigned int size = (unsigned)data.size();//input_tree->GetEntries(); bool thisValue=false; TString name=title.c_str(); name.Append("/O"); TBranch* output_branch_obj = input_tree->Branch( title.c_str(), &thisValue, name ); input_tree->SetEntries( (int)size ); for( unsigned int i=0; i< size; ++i ) { thisValue = data[i] ; output_branch_obj->Fill(); } } void HistoManager::CreateRTTree(const std::string& treeName) { rtTree = new TTree( treeName.data() , treeName.data() ); //Event variables rtTree->Branch( "evtNumber" , &rtEvtNumber ); rtTree->Branch( "pdgid", &rtPDGID ); rtTree->Branch( "layerNumber", &rtLayerNumber ); rtTree->Branch( "x", &rtX ); rtTree->Branch( "y", &rtY ); rtTree->Branch( "nElectrons", &rtElectrons ); // setup the binary files G4String name = fileName[0] + ".txt"; binout.open(name, std::ios::out); } void HistoManager::AddRTEvent(const G4int event, const G4int pdgid, const G4int layerNumber, const G4int x, const G4int y, const G4int nElectrons) { rtEvtNumber = event; rtPDGID = pdgid; rtLayerNumber = layerNumber; rtX = x; rtY = y; rtElectrons = nElectrons; if(rtTree == 0) return; else rtTree->Fill(); evtsnumbers.push_back(rtEvtNumber); pdgids.push_back(rtPDGID); layernumbers.push_back(rtLayerNumber); xs.push_back(rtX); ys.push_back(rtY); electrons.push_back(rtElectrons); if(xs.size()>10000) HistoManager::StoreRTInfo(); } void HistoManager::StoreRTInfo() { int nevts = xs.size(); for(int i=0; iAutoSave(); trackerTree->AutoSave(); // clear the vectors and starts again evtsnumbers.clear(); pdgids.clear(); layernumbers.clear(); xs.clear(); ys.clear(); electrons.clear(); } void HistoManager::CreateTrackerTree(const std::string& treeName, const int n_planes, const int n_strips ) { trackerTree = new TTree( treeName.data() , treeName.data() ); nStrips = n_strips; // used to set size of strip signal arrays Signal_x1 = new Float_t[nStrips]; Signal_u1 = new Float_t[nStrips]; Signal_v1 = new Float_t[nStrips]; Signal_x2 = new Float_t[nStrips]; Signal_u2 = new Float_t[nStrips]; Signal_v2 = new Float_t[nStrips]; Signal_TP1 = new Float_t[nStrips]; Signal_TP2 = new Float_t[nStrips]; Signal_TP3 = new Float_t[nStrips]; Signal_TP4 = new Float_t[nStrips]; Signal_x3 = new Float_t[nStrips]; Signal_u3 = new Float_t[nStrips]; Signal_v3 = new Float_t[nStrips]; Signal_x4 = new Float_t[nStrips]; Signal_u4 = new Float_t[nStrips]; Signal_v4 = new Float_t[nStrips]; //for ( Int_t plane = 0 ; plane < n_planes; ++plane ) {Det_rotations[plane] = det_angles[plane];} //for ( Int_t geom = 0 ; geom < 5; ++geom ) {Geom_params[geom] = geom_pos[geom];} for ( Int_t strip = 0 ; strip < nStrips ; ++strip ) { Signal_x1[strip] = 0; Signal_u1[strip] = 0; Signal_v1[strip] = 0; Signal_x2[strip] = 0; Signal_u2[strip] = 0; Signal_v2[strip] = 0; Signal_TP1[strip] = 0; Signal_TP2[strip] = 0; Signal_TP3[strip] = 0; Signal_TP4[strip] = 0; Signal_x3[strip] = 0; Signal_u3[strip] = 0; Signal_v3[strip] = 0; Signal_x4[strip] = 0; Signal_u4[strip] = 0; Signal_v4[strip] = 0; } char branch[50]; //Event variables trackerTree->Branch( "event_no" , &Event_no ); //contains parameters for position and angles that can be checked //against reconstruction config file //sprintf(branch, "det_rotations[%i]/F", n_planes); //trackerTree->Branch( "det_rotations", Det_rotations , branch ); //can't prevent from being written each event! //sprintf(branch, "geom_params[%i]/F", 5); //trackerTree->Branch( "geom_params", Geom_params , branch ); //can't prevent from being written each event! trackerTree->Branch( "ke_in" , &KE_in ); trackerTree->Branch( "ke_out" , &KE_out ); trackerTree->Branch( "ke_TP1" , &KE_TP1 ); trackerTree->Branch( "ke_TP2" , &KE_TP2 ); trackerTree->Branch( "ke_TP3" , &KE_TP3 ); trackerTree->Branch( "ke_TP4" , &KE_TP4 ); trackerTree->Branch( "truth_x_pos" , &Truth_x_pos ); trackerTree->Branch( "truth_y_pos" , &Truth_y_pos ); trackerTree->Branch( "truth_z_pos" , &Truth_z_pos ); trackerTree->Branch( "truthTheta_x" , &TruthTheta_x ); trackerTree->Branch( "truthTheta_y" , &TruthTheta_y ); trackerTree->Branch( "det_mult" , &Det_mult ); trackerTree->Branch( "hit_mult_x1" , &Hit_mult_x1 ); trackerTree->Branch( "hit_mult_u1" , &Hit_mult_u1 ); trackerTree->Branch( "hit_mult_v1" , &Hit_mult_v1 ); trackerTree->Branch( "hit_mult_x2" , &Hit_mult_x2 ); trackerTree->Branch( "hit_mult_u2" , &Hit_mult_u2 ); trackerTree->Branch( "hit_mult_v2" , &Hit_mult_v2 ); trackerTree->Branch( "hit_mult_TP1" , &Hit_mult_TP1 ); trackerTree->Branch( "hit_mult_TP2" , &Hit_mult_TP2 ); trackerTree->Branch( "hit_mult_TP3" , &Hit_mult_TP3 ); trackerTree->Branch( "hit_mult_TP4" , &Hit_mult_TP4 ); trackerTree->Branch( "hit_mult_x3" , &Hit_mult_x3 ); trackerTree->Branch( "hit_mult_u3" , &Hit_mult_u3 ); trackerTree->Branch( "hit_mult_v3" , &Hit_mult_v3 ); trackerTree->Branch( "hit_mult_x4" , &Hit_mult_x4 ); trackerTree->Branch( "hit_mult_u4" , &Hit_mult_u4 ); trackerTree->Branch( "hit_mult_v4" , &Hit_mult_v4 ); trackerTree->Branch( "clusterSize_x1" , &ClusterSize_x1 ); trackerTree->Branch( "clusterSize_u1" , &ClusterSize_u1 ); trackerTree->Branch( "clusterSize_v1" , &ClusterSize_v1 ); trackerTree->Branch( "clusterSize_x2" , &ClusterSize_x2 ); trackerTree->Branch( "clusterSize_u2" , &ClusterSize_u2 ); trackerTree->Branch( "clusterSize_v2" , &ClusterSize_v2 ); trackerTree->Branch( "clusterSize_TP1" , &ClusterSize_TP1 ); trackerTree->Branch( "clusterSize_TP2" , &ClusterSize_TP2 ); trackerTree->Branch( "clusterSize_TP3" , &ClusterSize_TP3 ); trackerTree->Branch( "clusterSize_TP4" , &ClusterSize_TP4 ); trackerTree->Branch( "clusterSize_x3" , &ClusterSize_x3 ); trackerTree->Branch( "clusterSize_u3" , &ClusterSize_u3 ); trackerTree->Branch( "clusterSize_v3" , &ClusterSize_v3 ); trackerTree->Branch( "clusterSize_x4" , &ClusterSize_x4 ); trackerTree->Branch( "clusterSize_u4" , &ClusterSize_u4 ); trackerTree->Branch( "clusterSize_v4" , &ClusterSize_v4 ); //Digit variables sprintf(branch, "signal_x1[%i]/F", nStrips); trackerTree->Branch( "signal_x1", Signal_x1 , branch ); sprintf(branch, "signal_u1[%i]/F", nStrips); trackerTree->Branch( "signal_u1", Signal_u1 , branch ); sprintf(branch, "signal_v1[%i]/F", nStrips); trackerTree->Branch( "signal_v1", Signal_v1 , branch ); sprintf(branch, "signal_x2[%i]/F", nStrips); trackerTree->Branch( "signal_x2", Signal_x2 , branch ); sprintf(branch, "signal_u2[%i]/F", nStrips); trackerTree->Branch( "signal_u2", Signal_u2 , branch ); sprintf(branch, "signal_v2[%i]/F", nStrips); trackerTree->Branch( "signal_v2", Signal_v2 , branch ); sprintf(branch, "signal_TP1[%i]/F", nStrips); trackerTree->Branch( "signal_TP1", Signal_TP1 , branch ); sprintf(branch, "signal_TP2[%i]/F", nStrips); trackerTree->Branch( "signal_TP2", Signal_TP2 , branch ); sprintf(branch, "signal_TP3[%i]/F", nStrips); trackerTree->Branch( "signal_TP3", Signal_TP3 , branch ); sprintf(branch, "signal_TP4[%i]/F", nStrips); trackerTree->Branch( "signal_TP4", Signal_TP4 , branch ); sprintf(branch, "signal_x3[%i]/F", nStrips); trackerTree->Branch( "signal_x3", Signal_x3 , branch ); sprintf(branch, "signal_u3[%i]/F", nStrips); trackerTree->Branch( "signal_u3", Signal_u3 , branch ); sprintf(branch, "signal_v3[%i]/F", nStrips); trackerTree->Branch( "signal_v3", Signal_v3 , branch ); sprintf(branch, "signal_x4[%i]/F", nStrips); trackerTree->Branch( "signal_x4", Signal_x4 , branch ); sprintf(branch, "signal_u4[%i]/F", nStrips); trackerTree->Branch( "signal_u4", Signal_u4 , branch ); sprintf(branch, "signal_v4[%i]/F", nStrips); trackerTree->Branch( "signal_v4", Signal_v4 , branch ); trackerTree->Branch( "strip_no_x1" , &Strip_no_x1 ); trackerTree->Branch( "strip_no_u1" , &Strip_no_u1 ); trackerTree->Branch( "strip_no_v1" , &Strip_no_v1 ); trackerTree->Branch( "strip_no_x2" , &Strip_no_x2 ); trackerTree->Branch( "strip_no_u2" , &Strip_no_u2 ); trackerTree->Branch( "strip_no_v2" , &Strip_no_v2 ); trackerTree->Branch( "strip_no_TP1" , &Strip_no_TP1 ); trackerTree->Branch( "strip_no_TP2" , &Strip_no_TP2 ); trackerTree->Branch( "strip_no_TP3" , &Strip_no_TP3 ); trackerTree->Branch( "strip_no_TP4" , &Strip_no_TP4 ); trackerTree->Branch( "strip_no_x3" , &Strip_no_x3 ); trackerTree->Branch( "strip_no_u3" , &Strip_no_u3 ); trackerTree->Branch( "strip_no_v3" , &Strip_no_v3 ); trackerTree->Branch( "strip_no_x4" , &Strip_no_x4 ); trackerTree->Branch( "strip_no_u4" , &Strip_no_u4 ); trackerTree->Branch( "strip_no_v4" , &Strip_no_v4 ); //Hit variables trackerTree->Branch( "edep_x1" , &Edep_x1 ); trackerTree->Branch( "edep_u1" , &Edep_u1 ); trackerTree->Branch( "edep_v1" , &Edep_v1 ); trackerTree->Branch( "edep_x2" , &Edep_x2 ); trackerTree->Branch( "edep_u2" , &Edep_u2 ); trackerTree->Branch( "edep_v2" , &Edep_v2 ); trackerTree->Branch( "edep_TP1" , &Edep_TP1 ); trackerTree->Branch( "edep_TP2" , &Edep_TP2 ); trackerTree->Branch( "edep_TP3" , &Edep_TP3 ); trackerTree->Branch( "edep_TP4" , &Edep_TP4 ); trackerTree->Branch( "edep_x3" , &Edep_x3 ); trackerTree->Branch( "edep_u3" , &Edep_u3 ); trackerTree->Branch( "edep_v3" , &Edep_v3 ); trackerTree->Branch( "edep_x4" , &Edep_x4 ); trackerTree->Branch( "edep_u4" , &Edep_u4 ); trackerTree->Branch( "edep_v4" , &Edep_v4 ); trackerTree->Branch( "particleZ_x1" , &ParticleZ_x1 ); trackerTree->Branch( "particleZ_u1" , &ParticleZ_u1 ); trackerTree->Branch( "particleZ_v1" , &ParticleZ_v1 ); trackerTree->Branch( "particleZ_x2" , &ParticleZ_x2 ); trackerTree->Branch( "particleZ_u2" , &ParticleZ_u2 ); trackerTree->Branch( "particleZ_v2" , &ParticleZ_v2 ); trackerTree->Branch( "particleZ_TP1" , &ParticleZ_TP1 ); trackerTree->Branch( "particleZ_TP2" , &ParticleZ_TP2 ); trackerTree->Branch( "particleZ_TP3" , &ParticleZ_TP3 ); trackerTree->Branch( "particleZ_TP4" , &ParticleZ_TP4 ); trackerTree->Branch( "particleZ_x3" , &ParticleZ_x3 ); trackerTree->Branch( "particleZ_u3" , &ParticleZ_u3 ); trackerTree->Branch( "particleZ_v3" , &ParticleZ_v3 ); trackerTree->Branch( "particleZ_x4" , &ParticleZ_x4 ); trackerTree->Branch( "particleZ_u4" , &ParticleZ_u4 ); trackerTree->Branch( "particleZ_v4" , &ParticleZ_v4 ); trackerTree->Branch( "isPrimaryParticle_x1" , &IsPrimaryParticle_x1 ); trackerTree->Branch( "isPrimaryParticle_u1" , &IsPrimaryParticle_u1 ); trackerTree->Branch( "isPrimaryParticle_v1" , &IsPrimaryParticle_v1 ); trackerTree->Branch( "isPrimaryParticle_x2" , &IsPrimaryParticle_x2 ); trackerTree->Branch( "isPrimaryParticle_u2" , &IsPrimaryParticle_u2 ); trackerTree->Branch( "isPrimaryParticle_v2" , &IsPrimaryParticle_v2 ); trackerTree->Branch( "isPrimaryParticle_TP1" , &IsPrimaryParticle_TP1 ); trackerTree->Branch( "isPrimaryParticle_TP2" , &IsPrimaryParticle_TP2 ); trackerTree->Branch( "isPrimaryParticle_TP3" , &IsPrimaryParticle_TP3 ); trackerTree->Branch( "isPrimaryParticle_TP4" , &IsPrimaryParticle_TP4 ); trackerTree->Branch( "isPrimaryParticle_x3" , &IsPrimaryParticle_x3 ); trackerTree->Branch( "isPrimaryParticle_u3" , &IsPrimaryParticle_u3 ); trackerTree->Branch( "isPrimaryParticle_v3" , &IsPrimaryParticle_v3 ); trackerTree->Branch( "isPrimaryParticle_x4" , &IsPrimaryParticle_x4 ); trackerTree->Branch( "isPrimaryParticle_u4" , &IsPrimaryParticle_u4 ); trackerTree->Branch( "isPrimaryParticle_v4" , &IsPrimaryParticle_v4 ); trackerTree->Branch( "isUp_x1" , &IsUp_x1 ); trackerTree->Branch( "isUp_u1" , &IsUp_u1 ); trackerTree->Branch( "isUp_v1" , &IsUp_v1 ); trackerTree->Branch( "isUp_x2" , &IsUp_x2 ); trackerTree->Branch( "isUp_u2" , &IsUp_u2 ); trackerTree->Branch( "isUp_v2" , &IsUp_v2 ); trackerTree->Branch( "isUp_TP1" , &IsUp_TP1 ); trackerTree->Branch( "isUp_TP2" , &IsUp_TP2 ); trackerTree->Branch( "isUp_TP3" , &IsUp_TP3 ); trackerTree->Branch( "isUp_TP4" , &IsUp_TP4 ); trackerTree->Branch( "isUp_x3" , &IsUp_x3 ); trackerTree->Branch( "isUp_u3" , &IsUp_u3 ); trackerTree->Branch( "isUp_v3" , &IsUp_v3 ); trackerTree->Branch( "isUp_x4" , &IsUp_x4 ); trackerTree->Branch( "isUp_u4" , &IsUp_u4 ); trackerTree->Branch( "isUp_v4" , &IsUp_v4 ); trackerTree->Branch( "trackNumber_x1" , &Track_no_x1 ); trackerTree->Branch( "trackNumber_u1" , &Track_no_u1 ); trackerTree->Branch( "trackNumber_v1" , &Track_no_v1 ); trackerTree->Branch( "trackNumber_x2" , &Track_no_x2 ); trackerTree->Branch( "trackNumber_u2" , &Track_no_u2 ); trackerTree->Branch( "trackNumber_v2" , &Track_no_v2 ); trackerTree->Branch( "trackNumber_TP1" , &Track_no_TP1 ); trackerTree->Branch( "trackNumber_TP2" , &Track_no_TP2 ); trackerTree->Branch( "trackNumber_TP3" , &Track_no_TP3 ); trackerTree->Branch( "trackNumber_TP4" , &Track_no_TP4 ); trackerTree->Branch( "trackNumber_x3" , &Track_no_x3 ); trackerTree->Branch( "trackNumber_u3" , &Track_no_u3 ); trackerTree->Branch( "trackNumber_v3" , &Track_no_v3 ); trackerTree->Branch( "trackNumber_x4" , &Track_no_x4 ); trackerTree->Branch( "trackNumber_u4" , &Track_no_u4 ); trackerTree->Branch( "trackNumber_v4" , &Track_no_v4 ); trackerTree->Branch( "particleName_x1" , &ParticleName_x1 ); trackerTree->Branch( "particleName_u1" , &ParticleName_u1 ); trackerTree->Branch( "particleName_v1" , &ParticleName_v1 ); trackerTree->Branch( "particleName_x2" , &ParticleName_x2 ); trackerTree->Branch( "particleName_u2" , &ParticleName_u2 ); trackerTree->Branch( "particleName_v2" , &ParticleName_v2 ); trackerTree->Branch( "particleName_TP1" , &ParticleName_TP1 ); trackerTree->Branch( "particleName_TP2" , &ParticleName_TP2 ); trackerTree->Branch( "particleName_TP3" , &ParticleName_TP3 ); trackerTree->Branch( "particleName_TP4" , &ParticleName_TP4 ); trackerTree->Branch( "particleName_x3" , &ParticleName_x3 ); trackerTree->Branch( "particleName_u3" , &ParticleName_u3 ); trackerTree->Branch( "particleName_v3" , &ParticleName_v3 ); trackerTree->Branch( "particleName_x4" , &ParticleName_x4 ); trackerTree->Branch( "particleName_u4" , &ParticleName_u4 ); trackerTree->Branch( "particleName_v4" , &ParticleName_v4 ); trackerTree->Branch( "x_pos_x1" , &X_pos_x1 ); trackerTree->Branch( "x_pos_u1" , &X_pos_u1 ); trackerTree->Branch( "x_pos_v1" , &X_pos_v1 ); trackerTree->Branch( "y_pos_x1" , &Y_pos_x1 ); trackerTree->Branch( "y_pos_u1" , &Y_pos_u1 ); trackerTree->Branch( "y_pos_v1" , &Y_pos_v1 ); trackerTree->Branch( "z_pos_x1" , &Z_pos_x1 ); trackerTree->Branch( "z_pos_u1" , &Z_pos_u1 ); trackerTree->Branch( "z_pos_v1" , &Z_pos_v1 ); trackerTree->Branch( "x_pos_x2" , &X_pos_x2 ); trackerTree->Branch( "x_pos_u2" , &X_pos_u2 ); trackerTree->Branch( "x_pos_v2" , &X_pos_v2 ); trackerTree->Branch( "y_pos_x2" , &Y_pos_x2 ); trackerTree->Branch( "y_pos_u2" , &Y_pos_u2 ); trackerTree->Branch( "y_pos_v2" , &Y_pos_v2 ); trackerTree->Branch( "z_pos_x2" , &Z_pos_x2 ); trackerTree->Branch( "z_pos_u2" , &Z_pos_u2 ); trackerTree->Branch( "z_pos_v2" , &Z_pos_v2 ); trackerTree->Branch( "x_pos_TP1" , &X_pos_TP1 ); trackerTree->Branch( "x_pos_TP2" , &X_pos_TP2 ); trackerTree->Branch( "x_pos_TP3" , &X_pos_TP3 ); trackerTree->Branch( "x_pos_TP4" , &X_pos_TP4 ); trackerTree->Branch( "y_pos_TP1" , &Y_pos_TP1 ); trackerTree->Branch( "y_pos_TP2" , &Y_pos_TP2 ); trackerTree->Branch( "y_pos_TP3" , &Y_pos_TP3 ); trackerTree->Branch( "y_pos_TP4" , &Y_pos_TP4 ); trackerTree->Branch( "z_pos_TP1" , &Z_pos_TP1 ); trackerTree->Branch( "z_pos_TP2" , &Z_pos_TP2 ); trackerTree->Branch( "z_pos_TP3" , &Z_pos_TP3 ); trackerTree->Branch( "z_pos_TP4" , &Z_pos_TP4 ); trackerTree->Branch( "x_pos_x3" , &X_pos_x3 ); trackerTree->Branch( "x_pos_u3" , &X_pos_u3 ); trackerTree->Branch( "x_pos_v3" , &X_pos_v3 ); trackerTree->Branch( "y_pos_x3" , &Y_pos_x3 ); trackerTree->Branch( "y_pos_u3" , &Y_pos_u3 ); trackerTree->Branch( "y_pos_v3" , &Y_pos_v3 ); trackerTree->Branch( "z_pos_x3" , &Z_pos_x3 ); trackerTree->Branch( "z_pos_u3" , &Z_pos_u3 ); trackerTree->Branch( "z_pos_v3" , &Z_pos_v3 ); trackerTree->Branch( "x_pos_x4" , &X_pos_x4 ); trackerTree->Branch( "x_pos_u4" , &X_pos_u4 ); trackerTree->Branch( "x_pos_v4" , &X_pos_v4 ); trackerTree->Branch( "y_pos_x4" , &Y_pos_x4 ); trackerTree->Branch( "y_pos_u4" , &Y_pos_u4 ); trackerTree->Branch( "y_pos_v4" , &Y_pos_v4 ); trackerTree->Branch( "z_pos_x4" , &Z_pos_x4 ); trackerTree->Branch( "z_pos_u4" , &Z_pos_u4 ); trackerTree->Branch( "z_pos_v4" , &Z_pos_v4 ); trackerTree->Branch( "theta_x1" , &Theta_x1 ); trackerTree->Branch( "theta_u1" , &Theta_u1 ); trackerTree->Branch( "theta_v1" , &Theta_v1 ); trackerTree->Branch( "theta_x2" , &Theta_x2 ); trackerTree->Branch( "theta_u2" , &Theta_u2 ); trackerTree->Branch( "theta_v2" , &Theta_v2 ); trackerTree->Branch( "theta_TP1" , &Theta_TP1 ); trackerTree->Branch( "theta_TP2" , &Theta_TP2 ); trackerTree->Branch( "theta_TP3" , &Theta_TP3 ); trackerTree->Branch( "theta_TP4" , &Theta_TP4 ); trackerTree->Branch( "theta_x3" , &Theta_x3 ); trackerTree->Branch( "theta_u3" , &Theta_u3 ); trackerTree->Branch( "theta_v3" , &Theta_v3 ); trackerTree->Branch( "theta_x4" , &Theta_x4 ); trackerTree->Branch( "theta_u4" , &Theta_u4 ); trackerTree->Branch( "theta_v4" , &Theta_v4 ); } void HistoManager::AddTrackerEvent(const G4int event, const StripHitCollection* const hits_x1, const StripHitCollection* const hits_u1, const StripHitCollection* const hits_v1, const StripHitCollection* const hits_x2, const StripHitCollection* const hits_u2, const StripHitCollection* const hits_v2, const StripHitCollection* const hits_TP1, const StripHitCollection* const hits_TP2, const StripHitCollection* const hits_TP3, const StripHitCollection* const hits_TP4, const StripHitCollection* const hits_x3, const StripHitCollection* const hits_u3, const StripHitCollection* const hits_v3, const StripHitCollection* const hits_x4, const StripHitCollection* const hits_u4, const StripHitCollection* const hits_v4, const StripDigiCollection* const digits ,const G4ThreeVector& primPos, const G4ThreeVector& primMom, const G4float K_E_in, const G4float K_E_out) { //If root TTree is not created ends if ( trackerTree == 0 ) {return;} //Initialise variables Event_no = event; Det_mult = 0; KE_in = K_E_in; KE_out = K_E_out; //set at the end with truth variables ClusterSize_x1 = 0; ClusterSize_u1 = 0; ClusterSize_v1 = 0; ClusterSize_x2 = 0; ClusterSize_u2 = 0; ClusterSize_v2 = 0; ClusterSize_TP1 = 0; ClusterSize_TP2 = 0; ClusterSize_TP3 = 0; ClusterSize_TP4 = 0; ClusterSize_x3 = 0; ClusterSize_u3 = 0; ClusterSize_v3 = 0; ClusterSize_x4 = 0; ClusterSize_u4 = 0; ClusterSize_v4 = 0; Strip_no_x1.clear(); Strip_no_u1.clear(); Strip_no_v1.clear(); Strip_no_x2.clear(); Strip_no_u2.clear(); Strip_no_v2.clear(); Strip_no_TP1.clear(); Strip_no_TP2.clear(); Strip_no_TP3.clear(); Strip_no_TP4.clear(); Strip_no_x3.clear(); Strip_no_u3.clear(); Strip_no_v3.clear(); Strip_no_x4.clear(); Strip_no_u4.clear(); Strip_no_v4.clear(); //Store Digits information if ( digits ) { G4int nDigits = digits->entries(); G4int is_det_x1, is_det_u1, is_det_v1 ,is_det_x2, is_det_u2, is_det_v2, is_det_x3, is_det_u3, is_det_v3, is_det_x4, is_det_u4, is_det_v4; is_det_x1 = 0; is_det_u1 = 0; is_det_v1 = 0; is_det_x2 = 0; is_det_u2 = 0; is_det_v2 = 0; is_det_x3 = 0; is_det_u3 = 0; is_det_v3 = 0; is_det_x4 = 0; is_det_u4 = 0; is_det_v4 = 0; for ( G4int d = 0 ; d( digits->GetDigi( d ) ); G4int stripNum = digi->GetStripNumber(); //Safety check if ( stripNum >= nStrips ) { G4cerr << "Digi Error: Strip number "<< stripNum << " expected max value:" << nStrips << G4endl; continue;//Go to next digit } G4int planeNum = digi->GetPlaneNumber(); if ( planeNum == 0 ) { Signal_x1[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_x1[ stripNum ] > 0) {is_det_x1 = 1; ClusterSize_x1++; Strip_no_x1.push_back(stripNum);} } else if ( planeNum == 1 ) { Signal_u1[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_u1[ stripNum ] > 0) {is_det_u1 = 1; ClusterSize_u1++; Strip_no_u1.push_back(stripNum);} } else if ( planeNum == 2 ) { Signal_v1[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_v1[ stripNum ] > 0) {is_det_v1 = 1; ClusterSize_v1++; Strip_no_v1.push_back(stripNum);} } else if ( planeNum == 3 ) { Signal_x2[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_x2[ stripNum ] > 0) {is_det_x2 = 1; ClusterSize_x2++; Strip_no_x2.push_back(stripNum);} } else if ( planeNum == 4 ) { Signal_u2[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_u2[ stripNum ] > 0) {is_det_u2 = 1; ClusterSize_u2++; Strip_no_u2.push_back(stripNum);} } else if ( planeNum == 5 ) { Signal_v2[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_v2[ stripNum ] > 0) {is_det_v2 = 1; ClusterSize_v2++; Strip_no_v2.push_back(stripNum);} } else if ( planeNum == 6 ) { Signal_TP1[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_TP1[ stripNum ] > 0) {ClusterSize_TP1++; Strip_no_TP1.push_back(stripNum);} } else if ( planeNum == 7 ) { Signal_TP2[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_TP2[ stripNum ] > 0) {ClusterSize_TP2++; Strip_no_TP2.push_back(stripNum);} } else if ( planeNum == 8 ) { Signal_TP3[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_TP3[ stripNum ] > 0) {ClusterSize_TP3++; Strip_no_TP3.push_back(stripNum);} } else if ( planeNum == 9 ) { Signal_TP4[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_TP4[ stripNum ] > 0) {ClusterSize_TP4++; Strip_no_TP4.push_back(stripNum);} } else if ( planeNum == 10 ) { Signal_x3[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_x3[ stripNum ] > 0) {is_det_x3 = 1; ClusterSize_x3++; Strip_no_x3.push_back(stripNum);} } else if ( planeNum == 11 ) { Signal_u3[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_u3[ stripNum ] > 0) {is_det_u3 = 1; ClusterSize_u3++; Strip_no_u3.push_back(stripNum);} } else if ( planeNum == 12 ) { Signal_v3[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_v3[ stripNum ] > 0) {is_det_v3 = 1; ClusterSize_v3++; Strip_no_v3.push_back(stripNum);} } else if ( planeNum == 13 ) { Signal_x4[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_x4[ stripNum ] > 0) {is_det_x4 = 1; ClusterSize_x4++; Strip_no_x4.push_back(stripNum);} } else if ( planeNum == 14 ) { Signal_u4[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_u4[ stripNum ] > 0) {is_det_u4 = 1; ClusterSize_u4++; Strip_no_u4.push_back(stripNum);} } else if ( planeNum == 15 ) { Signal_v4[ stripNum ] = static_cast(digi->GetCharge()); if(Signal_v4[ stripNum ] > 0) {is_det_v4 = 1; ClusterSize_v4++; Strip_no_v4.push_back(stripNum);} } else { G4cerr << "Digi Error: Plane number not set correctly in DetectorConstruction.cc" << G4endl; } } //Generate detector multiplicity, truth planes not included. Det_mult = is_det_x1 + is_det_u1 + is_det_v1 + is_det_x2 + is_det_u2 + is_det_v2 + is_det_x3 + is_det_u3 + is_det_v3 + is_det_x4 + is_det_u4 + is_det_v4; } else { G4cerr << "Error: No digi collection passed to RootSaver" << G4endl; } //Store Hits information if ( hits_x1 || hits_u1 || hits_v1 || hits_x2 || hits_u2 || hits_v2 || hits_x3 || hits_u3 || hits_v3 || hits_x4 || hits_u4 || hits_v4 ) { //Set defaults Hit_mult_x1 = 0; Hit_mult_u1 = 0; Hit_mult_v1 = 0; Hit_mult_x2 = 0; Hit_mult_u2 = 0; Hit_mult_v2 = 0; Hit_mult_TP1 = 0; Hit_mult_TP2 = 0; Hit_mult_TP3 = 0; Hit_mult_TP4 = 0; Hit_mult_x3 = 0; Hit_mult_u3 = 0; Hit_mult_v3 = 0; Hit_mult_x4 = 0; Hit_mult_u4 = 0; Hit_mult_v4 = 0; Edep_x1.clear(); Edep_u1.clear(); Edep_v1.clear(); Edep_x2.clear(); Edep_u2.clear(); Edep_v2.clear(); Edep_TP1.clear(); Edep_TP2.clear(); Edep_TP3.clear(); Edep_TP4.clear(); Edep_x3.clear(); Edep_u3.clear(); Edep_v3.clear(); Edep_x4.clear(); Edep_u4.clear(); Edep_v4.clear(); G4float th_x1 = 0; G4float th_u1 = 0; G4float th_v1 = 0; G4float th_x2 = 0; G4float th_u2 = 0; G4float th_v2 = 0; G4float th_TP1 = 0; G4float th_TP2 = 0; G4float th_TP3 = 0; G4float th_TP4 = 0; G4float th_x3 = 0; G4float th_u3 = 0; G4float th_v3 = 0; G4float th_x4 = 0; G4float th_u4 = 0; G4float th_v4 = 0; G4float x0 = 0; G4float y0 = 0; G4float z0 = 0; G4float x1 = 0; G4float y1 = 0; G4float z1 = 0; IsPrimaryParticle_x1.clear(); IsPrimaryParticle_u1.clear(); IsPrimaryParticle_v1.clear(); IsPrimaryParticle_x2.clear(); IsPrimaryParticle_u2.clear(); IsPrimaryParticle_v2.clear(); IsPrimaryParticle_TP1.clear(); IsPrimaryParticle_TP2.clear(); IsPrimaryParticle_TP3.clear(); IsPrimaryParticle_TP4.clear(); IsPrimaryParticle_x3.clear(); IsPrimaryParticle_u3.clear(); IsPrimaryParticle_v3.clear(); IsPrimaryParticle_x4.clear(); IsPrimaryParticle_u4.clear(); IsPrimaryParticle_v4.clear(); IsUp_x1.clear(); IsUp_u1.clear(); IsUp_v1.clear(); IsUp_x2.clear(); IsUp_u2.clear(); IsUp_v2.clear(); IsUp_TP1.clear(); IsUp_TP2.clear(); IsUp_TP3.clear(); IsUp_TP4.clear(); IsUp_x3.clear(); IsUp_u3.clear(); IsUp_v3.clear(); IsUp_x4.clear(); IsUp_u4.clear(); IsUp_v4.clear(); Track_no_x1.clear(); Track_no_u1.clear(); Track_no_v1.clear(); Track_no_x2.clear(); Track_no_u2.clear(); Track_no_v2.clear(); Track_no_TP1.clear(); Track_no_TP2.clear(); Track_no_TP3.clear(); Track_no_TP4.clear(); Track_no_x3.clear(); Track_no_u3.clear(); Track_no_v3.clear(); Track_no_x4.clear(); Track_no_u4.clear(); Track_no_v4.clear(); ParticleName_x1.clear(); ParticleName_u1.clear(); ParticleName_v1.clear(); ParticleName_x2.clear(); ParticleName_u2.clear(); ParticleName_v2.clear(); ParticleName_TP1.clear(); ParticleName_TP2.clear(); ParticleName_TP3.clear(); ParticleName_TP4.clear(); ParticleName_x3.clear(); ParticleName_u3.clear(); ParticleName_v3.clear(); ParticleName_x4.clear(); ParticleName_u4.clear(); ParticleName_v4.clear(); X_pos_x1.clear(); X_pos_u1.clear(); X_pos_v1.clear(); Y_pos_x1.clear(); Y_pos_u1.clear(); Y_pos_v1.clear(); Z_pos_x1.clear(); Z_pos_u1.clear(); Z_pos_v1.clear(); X_pos_x2.clear(); X_pos_u2.clear(); X_pos_v2.clear(); Y_pos_x2.clear(); Y_pos_u2.clear(); Y_pos_v2.clear(); Z_pos_x2.clear(); Z_pos_u2.clear(); Z_pos_v2.clear(); X_pos_TP1.clear(); X_pos_TP2.clear(); X_pos_TP3.clear(); X_pos_TP4.clear(); Y_pos_TP1.clear(); Y_pos_TP2.clear(); Y_pos_TP3.clear(); Y_pos_TP4.clear(); Z_pos_TP1.clear(); Z_pos_TP2.clear(); Z_pos_TP3.clear(); Z_pos_TP4.clear(); X_pos_x3.clear(); X_pos_u3.clear(); X_pos_v3.clear(); Y_pos_x3.clear(); Y_pos_u3.clear(); Y_pos_v3.clear(); Z_pos_x3.clear(); Z_pos_u3.clear(); Z_pos_v3.clear(); X_pos_x4.clear(); X_pos_u4.clear(); X_pos_v4.clear(); Y_pos_x4.clear(); Y_pos_u4.clear(); Y_pos_v4.clear(); Z_pos_x4.clear(); Z_pos_u4.clear(); Z_pos_v4.clear(); Theta_x1.clear(); Theta_u1.clear(); Theta_v1.clear(); Theta_x2.clear(); Theta_u2.clear(); Theta_v2.clear(); Theta_TP1.clear(); Theta_TP2.clear(); Theta_TP3.clear(); Theta_TP4.clear(); Theta_x3.clear(); Theta_u3.clear(); Theta_v3.clear(); Theta_x4.clear(); Theta_u4.clear(); Theta_v4.clear(); ParticleZ_x1.clear(); ParticleZ_u1.clear(); ParticleZ_v1.clear(); ParticleZ_x2.clear(); ParticleZ_u2.clear(); ParticleZ_v2.clear(); ParticleZ_TP1.clear(); ParticleZ_TP2.clear(); ParticleZ_TP3.clear(); ParticleZ_TP4.clear(); ParticleZ_x3.clear(); ParticleZ_u3.clear(); ParticleZ_v3.clear(); ParticleZ_x4.clear(); ParticleZ_u4.clear(); ParticleZ_v4.clear(); //Loop on all hits, to obtain energy, pos and angle //Position is weighted average of hit x(), see SensitiveDetector.cc //********************************** FIRST TRACKER MODULE **********************************// G4int nHits = hits_x1->entries(); for ( G4int h = 0 ; (h( hits_x1->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; th_x1 = std::atan( (x1-x0) / (z1-z0) ); th_x1 /= mrad; //only one element per track (particle) per sensitive detector is recorded with the index of each element in the vector being equal to the track no. //this could be changed such that multiple hits per track (particle) per sensitive detector are recorded. //hit->GetTrackNumber()) == (Track_no_x1.back()) ensures the same element is overwritten each time if the hit belongs to the same track (particle), //except for energy which is not overwritten but added up. if( Track_no_x1.size()>0 && (hit->GetTrackNumber()) == (Track_no_x1.back()) ) //ensures each element represents a track, element [0] = track 1 { //Use printouts below to check that all arrays are the same size. /*if( Track_no_x1.size() != ParticleName_x1.size() ) G4cout << "Track_no_x1.size() != ParticleName_x1.size()" << G4endl; if( Track_no_x1.size() != IsPrimaryParticle_x1.size() ) G4cout << "Track_no_x1.size() != IsPrimaryParticle_x1.size()" << G4endl; if( Track_no_x1.size() != X_pos_x1.size() ) G4cout << "Track_no_x1.size() != X_pos_x1.size()" << G4endl; if( Track_no_x1.size() != Y_pos_x1.size() ) G4cout << "Track_no_x1.size() != Y_pos_x1.size()" << G4endl; if( Track_no_x1.size() != Z_pos_x1.size() ) G4cout << "Track_no_x1.size() != Z_pos_x1.size()" << G4endl; if( Track_no_x1.size() != Theta_x1.size() ) G4cout << "Track_no_x1.size() != Theta_x1.size()" << G4endl; if( Track_no_x1.size() != Edep_x1.size() ) G4cout << "Track_no_x1.size() != Edep_x1.size()" << G4endl; if( Track_no_x1.size() != ParticleZ_x1.size() ) G4cout << "Track_no_x1.size() != ParticleZ_x1.size()" << G4endl; if( Track_no_x1.size() != IsUp_x1.size() ) G4cout << "Track_no_x1.size() != IsUp_x1.size()" << G4endl;*/ /*ParticleName_x1[(ParticleName_x1.size())-1] = hit->GetParticleName(); Track_no_x1[(Track_no_x1.size())-1] = hit->GetTrackNumber(); IsPrimaryParticle_x1[(IsPrimaryParticle_x1.size())-1] = hit->GetIsPrimary(); X_pos_x1[(X_pos_x1.size())-1] = x1; Y_pos_x1[(Y_pos_x1.size())-1] = y1; Z_pos_x1[(Z_pos_x1.size())-1] = z1; Theta_x1[(Theta_x1.size())-1] = th_x1; Edep_x1[(Edep_x1.size())-1] += edep; ParticleZ_x1[(ParticleZ_x1.size())-1] = hit->GetParticleZ(); IsUp_x1[(IsUp_x1.size())-1] = (y1 > 0) ? true : false;*/ // Track_no_x1.size()-1 gives means elements from different vectors correspond to the // same particle when the same element no. is chosen. For each events all vectors are same // size, this can be verified with print out statement above. ParticleName_x1[Track_no_x1.size()-1] = hit->GetParticleName(); Track_no_x1[Track_no_x1.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_x1[Track_no_x1.size()-1] = hit->GetIsPrimary(); X_pos_x1[Track_no_x1.size()-1] = x1; Y_pos_x1[Track_no_x1.size()-1] = y1; Z_pos_x1[Track_no_x1.size()-1] = z1; Theta_x1[Track_no_x1.size()-1] = th_x1; Edep_x1[Track_no_x1.size()-1] += edep; ParticleZ_x1[Track_no_x1.size()-1] = hit->GetParticleZ(); IsUp_x1[Track_no_x1.size()-1] = (y1 > 0) ? true : false; } else { //G4cout << "Track_no_x1 = " << hit->GetTrackNumber() << G4endl; ParticleName_x1.push_back(hit->GetParticleName()); Track_no_x1.push_back(hit->GetTrackNumber()); IsPrimaryParticle_x1.push_back(hit->GetIsPrimary()); X_pos_x1.push_back(x1); Y_pos_x1.push_back(y1); Z_pos_x1.push_back(z1); Theta_x1.push_back(th_x1); Edep_x1.push_back(edep); ParticleZ_x1.push_back(hit->GetParticleZ()); IsUp_x1.push_back((y1 > 0) ? true : false); } Hit_mult_x1 = X_pos_x1.size(); //G4int planeNum = hit->GetPlaneNumber(); /* if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det x1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det x1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_u1->entries(); Edep_u1.clear(); for ( G4int h = 0 ; (h( hits_u1->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; //Edep_u1 += edep; th_u1 = std::atan( (x1-x0) / (z1-z0) ); th_u1 /= mrad; //if( (hit->GetTrackNumber()) <= (ParticleName_u1.size()) ) //ensures each element represents a track, element [0] = track 1 if( Track_no_u1.size()>0 && (hit->GetTrackNumber()) == (Track_no_u1.back()) ) //ensures each element represents a track, element [0] = track 1 { /*ParticleName_u1[(ParticleName_u1.size())-1] = hit->GetParticleName(); Track_no_u1[(Track_no_u1.size())-1] = hit->GetTrackNumber(); IsPrimaryParticle_u1[(IsPrimaryParticle_u1.size())-1] = hit->GetIsPrimary(); X_pos_u1[(X_pos_u1.size())-1] = x1; Y_pos_u1[(Y_pos_u1.size())-1] = y1; Z_pos_u1[(Z_pos_u1.size())-1] = z1; Theta_u1[(Theta_u1.size())-1] = th_u1; Edep_u1[(Edep_u1.size())-1] += edep; ParticleZ_u1[(ParticleZ_u1.size())-1] = hit->GetParticleZ(); IsUp_u1[(IsUp_u1.size())-1] = (y1 > 0) ? true : false; //Needs to be modified as plane is at stereo angle != 0*/ ParticleName_u1[Track_no_u1.size()-1] = hit->GetParticleName(); Track_no_u1[Track_no_u1.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_u1[Track_no_u1.size()-1] = hit->GetIsPrimary(); X_pos_u1[Track_no_u1.size()-1] = x1; Y_pos_u1[Track_no_u1.size()-1] = y1; Z_pos_u1[Track_no_u1.size()-1] = z1; Theta_u1[Track_no_u1.size()-1] = th_u1; Edep_u1[Track_no_u1.size()-1] += edep; ParticleZ_u1[Track_no_u1.size()-1] = hit->GetParticleZ(); IsUp_u1[Track_no_u1.size()-1] = (y1 > 0) ? true : false; //Needs to be modified as plane is at stereo angle != 0 } else { ParticleName_u1.push_back(hit->GetParticleName()); Track_no_u1.push_back(hit->GetTrackNumber()); IsPrimaryParticle_u1.push_back(hit->GetIsPrimary()); X_pos_u1.push_back(x1); Y_pos_u1.push_back(y1); Z_pos_u1.push_back(z1); Theta_u1.push_back(th_u1); Edep_u1.push_back(edep); ParticleZ_u1.push_back(hit->GetParticleZ()); IsUp_u1.push_back((y1 > 0) ? true : false); //Needs to be modified as plane is at stereo angle != 0 } Hit_mult_u1 = X_pos_u1.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det u1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det u1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_v1->entries(); Edep_v1.clear(); for ( G4int h = 0 ; (h( hits_v1->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_v1 += edep; th_v1 = std::atan( (x1-x0) / (z1-z0) ); th_v1 /= mrad; if( Track_no_v1.size()>0 && (hit->GetTrackNumber()) == (Track_no_v1.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_v1[Track_no_v1.size()-1] = hit->GetParticleName(); Track_no_v1[Track_no_v1.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_v1[Track_no_v1.size()-1] = hit->GetIsPrimary(); X_pos_v1[Track_no_v1.size()-1] = x1; Y_pos_v1[Track_no_v1.size()-1] = y1; Z_pos_v1[Track_no_v1.size()-1] = z1; Theta_v1[Track_no_v1.size()-1] = th_v1; Edep_v1[Track_no_v1.size()-1] += edep; ParticleZ_v1[Track_no_v1.size()-1] = hit->GetParticleZ(); IsUp_v1[Track_no_v1.size()-1] = (y1 > 0) ? true : false; //Needs to be modified as plane is at stereo angle != 0 //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v1 for particle: " << hit->GetParticleName() << ", Edep_v1.back() = " << Edep_v1.back() << G4endl; } else { ParticleName_v1.push_back(hit->GetParticleName()); Track_no_v1.push_back(hit->GetTrackNumber()); IsPrimaryParticle_v1.push_back(hit->GetIsPrimary()); X_pos_v1.push_back(x1); Y_pos_v1.push_back(y1); Z_pos_v1.push_back(z1); Theta_v1.push_back(th_v1); Edep_v1.push_back(edep); ParticleZ_v1.push_back(hit->GetParticleZ()); IsUp_v1.push_back((y1 > 0) ? true : false); //Needs to be modified as plane is at stereo angle != 0 //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v1 for particle: " << hit->GetParticleName() << ", Edep_v1[(hit->GetTrackNumber())-1] = " << Edep_v1.back() << G4endl; } Hit_mult_v1 = X_pos_v1.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det v1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det v1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** SECOND TRACKER MODULE **********************************// nHits = hits_x2->entries(); for ( G4int h = 0 ; (h( hits_x2->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; th_x2 = std::atan( (x1-x0) / (z1-z0) ); th_x2 /= mrad; if( Track_no_x2.size()>0 && (hit->GetTrackNumber()) == (Track_no_x2.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_x2[Track_no_x2.size()-1] = hit->GetParticleName(); Track_no_x2[Track_no_x2.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_x2[Track_no_x2.size()-1] = hit->GetIsPrimary(); X_pos_x2[Track_no_x2.size()-1] = x1; Y_pos_x2[Track_no_x2.size()-1] = y1; Z_pos_x2[Track_no_x2.size()-1] = z1; Theta_x2[Track_no_x2.size()-1] = th_x2; Edep_x2[Track_no_x2.size()-1] += edep; ParticleZ_x2[Track_no_x2.size()-1] = hit->GetParticleZ(); IsUp_x2[Track_no_x2.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_x2.push_back(hit->GetParticleName()); Track_no_x2.push_back(hit->GetTrackNumber()); IsPrimaryParticle_x2.push_back(hit->GetIsPrimary()); X_pos_x2.push_back(x1); Y_pos_x2.push_back(y1); Z_pos_x2.push_back(z1); Theta_x2.push_back(th_x2); Edep_x2.push_back(edep); ParticleZ_x2.push_back(hit->GetParticleZ()); IsUp_x2.push_back((y1 > 0) ? true : false); } Hit_mult_x2 = X_pos_x2.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det x2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det x2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_u2->entries(); Edep_u2.clear(); for ( G4int h = 0 ; (h( hits_u2->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; //Edep_u2 += edep; th_u2 = std::atan( (x1-x0) / (z1-z0) ); th_u2 /= mrad; //if( (hit->GetTrackNumber()) <= (ParticleName_u2.size()) ) //ensures each element represents a track, element [0] = track 1 if( Track_no_u2.size()>0 && (hit->GetTrackNumber()) == (Track_no_u2.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_u2[Track_no_u2.size()-1] = hit->GetParticleName(); Track_no_u2[Track_no_u2.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_u2[Track_no_u2.size()-1] = hit->GetIsPrimary(); X_pos_u2[Track_no_u2.size()-1] = x1; Y_pos_u2[Track_no_u2.size()-1] = y1; Z_pos_u2[Track_no_u2.size()-1] = z1; Theta_u2[Track_no_u2.size()-1] = th_u2; Edep_u2[Track_no_u2.size()-1] += edep; ParticleZ_u2[Track_no_u2.size()-1] = hit->GetParticleZ(); IsUp_u2[Track_no_u2.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_u2.push_back(hit->GetParticleName()); Track_no_u2.push_back(hit->GetTrackNumber()); IsPrimaryParticle_u2.push_back(hit->GetIsPrimary()); X_pos_u2.push_back(x1); Y_pos_u2.push_back(y1); Z_pos_u2.push_back(z1); Theta_u2.push_back(th_u2); Edep_u2.push_back(edep); ParticleZ_u2.push_back(hit->GetParticleZ()); IsUp_u2.push_back((y1 > 0) ? true : false); } Hit_mult_u2 = X_pos_u2.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det u2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det u2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_v2->entries(); Edep_v2.clear(); for ( G4int h = 0 ; (h( hits_v2->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_v2 += edep; th_v2 = std::atan( (x1-x0) / (z1-z0) ); th_v2 /= mrad; if( Track_no_v2.size()>0 && (hit->GetTrackNumber()) == (Track_no_v2.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_v2[Track_no_v2.size()-1] = hit->GetParticleName(); Track_no_v2[Track_no_v2.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_v2[(IsPrimaryParticle_v2.size())-1] = hit->GetIsPrimary(); X_pos_v2[Track_no_v2.size()-1] = x1; Y_pos_v2[Track_no_v2.size()-1] = y1; Z_pos_v2[Track_no_v2.size()-1] = z1; Theta_v2[Track_no_v2.size()-1] = th_v2; Edep_v2[Track_no_v2.size()-1] += edep; ParticleZ_v2[Track_no_v2.size()-1] = hit->GetParticleZ(); IsUp_v2[Track_no_v2.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v2 for particle: " << hit->GetParticleName() << ", Edep_v2.back() = " << Edep_v2.back() << G4endl; } else { ParticleName_v2.push_back(hit->GetParticleName()); Track_no_v2.push_back(hit->GetTrackNumber()); IsPrimaryParticle_v2.push_back(hit->GetIsPrimary()); X_pos_v2.push_back(x1); Y_pos_v2.push_back(y1); Z_pos_v2.push_back(z1); Theta_v2.push_back(th_v2); Edep_v2.push_back(edep); ParticleZ_v2.push_back(hit->GetParticleZ()); IsUp_v2.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v2 for particle: " << hit->GetParticleName() << ", Edep_v2[(hit->GetTrackNumber())-1] = " << Edep_v2.back() << G4endl; } Hit_mult_v2 = X_pos_v2.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det v2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det v2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** TRUTH PLANE 1 **********************************// nHits = hits_TP1->entries(); Edep_TP1.clear(); for ( G4int h = 0 ; (h( hits_TP1->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; //Use to obtain exit KE of primary proton after passing through tracker and phantom if ( hit->GetIsPrimary() == true ){KE_TP1 = static_cast(hit->GetKE());} G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_TP1 += edep; th_TP1 = std::atan( (x1-x0) / (z1-z0) ); th_TP1 /= mrad; if( Track_no_TP1.size()>0 && (hit->GetTrackNumber()) == (Track_no_TP1.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_TP1[Track_no_TP1.size()-1] = hit->GetParticleName(); Track_no_TP1[Track_no_TP1.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_TP1[Track_no_TP1.size()-1] = hit->GetIsPrimary(); X_pos_TP1[Track_no_TP1.size()-1] = x1; Y_pos_TP1[Track_no_TP1.size()-1] = y1; Z_pos_TP1[Track_no_TP1.size()-1] = z1; Theta_TP1[Track_no_TP1.size()-1] = th_TP1; Edep_TP1[Track_no_TP1.size()-1] += edep; ParticleZ_TP1[Track_no_TP1.size()-1] = hit->GetParticleZ(); IsUp_TP1[Track_no_TP1.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP1 for particle: " << hit->GetParticleName() << ", Edep_TP1.back() = " << Edep_TP1.back() << G4endl; } else { ParticleName_TP1.push_back(hit->GetParticleName()); Track_no_TP1.push_back(hit->GetTrackNumber()); IsPrimaryParticle_TP1.push_back(hit->GetIsPrimary()); X_pos_TP1.push_back(x1); Y_pos_TP1.push_back(y1); Z_pos_TP1.push_back(z1); Theta_TP1.push_back(th_TP1); Edep_TP1.push_back(edep); ParticleZ_TP1.push_back(hit->GetParticleZ()); IsUp_TP1.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP1 for particle: " << hit->GetParticleName() << ", Edep_TP1[(hit->GetTrackNumber())-1] = " << Edep_TP1.back() << G4endl; } Hit_mult_TP1 = X_pos_TP1.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det TP1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det TP1 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP1.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** TRUTH PLANE 2 **********************************// nHits = hits_TP2->entries(); Edep_TP2.clear(); for ( G4int h = 0 ; (h( hits_TP2->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; //Use to obtain exit KE of primary proton after passing through tracker and phantom if ( hit->GetIsPrimary() == true ){KE_TP2 = static_cast(hit->GetKE());} G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_TP2 += edep; th_TP2 = std::atan( (x1-x0) / (z1-z0) ); th_TP2 /= mrad; if( Track_no_TP2.size()>0 && (hit->GetTrackNumber()) == (Track_no_TP2.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_TP2[Track_no_TP2.size()-1] = hit->GetParticleName(); Track_no_TP2[Track_no_TP2.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_TP2[Track_no_TP2.size()-1] = hit->GetIsPrimary(); X_pos_TP2[Track_no_TP2.size()-1] = x1; Y_pos_TP2[Track_no_TP2.size()-1] = y1; Z_pos_TP2[Track_no_TP2.size()-1] = z1; Theta_TP2[Track_no_TP2.size()-1] = th_TP2; Edep_TP2[Track_no_TP2.size()-1] += edep; ParticleZ_TP2[Track_no_TP2.size()-1] = hit->GetParticleZ(); IsUp_TP2[Track_no_TP2.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP2 for particle: " << hit->GetParticleName() << ", Edep_TP2.back() = " << Edep_TP2.back() << G4endl; } else { ParticleName_TP2.push_back(hit->GetParticleName()); Track_no_TP2.push_back(hit->GetTrackNumber()); IsPrimaryParticle_TP2.push_back(hit->GetIsPrimary()); X_pos_TP2.push_back(x1); Y_pos_TP2.push_back(y1); Z_pos_TP2.push_back(z1); Theta_TP2.push_back(th_TP2); Edep_TP2.push_back(edep); ParticleZ_TP2.push_back(hit->GetParticleZ()); IsUp_TP2.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP2 for particle: " << hit->GetParticleName() << ", Edep_TP2[(hit->GetTrackNumber())-1] = " << Edep_TP2.back() << G4endl; } Hit_mult_TP2 = X_pos_TP2.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det TP2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det TP2 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP2.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** TRUTH PLANE 3 **********************************// nHits = hits_TP3->entries(); Edep_TP3.clear(); for ( G4int h = 0 ; (h( hits_TP3->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; //Use to obtain exit KE of primary proton after passing through tracker and phantom if ( hit->GetIsPrimary() == true ){KE_TP3 = static_cast(hit->GetKE());} G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_TP3 += edep; th_TP3 = std::atan( (x1-x0) / (z1-z0) ); th_TP3 /= mrad; if( Track_no_TP3.size()>0 && (hit->GetTrackNumber()) == (Track_no_TP3.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_TP3[Track_no_TP3.size()-1] = hit->GetParticleName(); Track_no_TP3[Track_no_TP3.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_TP3[Track_no_TP3.size()-1] = hit->GetIsPrimary(); X_pos_TP3[Track_no_TP3.size()-1] = x1; Y_pos_TP3[Track_no_TP3.size()-1] = y1; Z_pos_TP3[Track_no_TP3.size()-1] = z1; Theta_TP3[Track_no_TP3.size()-1] = th_TP3; Edep_TP3[Track_no_TP3.size()-1] += edep; ParticleZ_TP3[Track_no_TP3.size()-1] = hit->GetParticleZ(); IsUp_TP3[Track_no_TP3.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP3 for particle: " << hit->GetParticleName() << ", Edep_TP3.back() = " << Edep_TP3.back() << G4endl; } else { ParticleName_TP3.push_back(hit->GetParticleName()); Track_no_TP3.push_back(hit->GetTrackNumber()); IsPrimaryParticle_TP3.push_back(hit->GetIsPrimary()); X_pos_TP3.push_back(x1); Y_pos_TP3.push_back(y1); Z_pos_TP3.push_back(z1); Theta_TP3.push_back(th_TP3); Edep_TP3.push_back(edep); ParticleZ_TP3.push_back(hit->GetParticleZ()); IsUp_TP3.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP3 for particle: " << hit->GetParticleName() << ", Edep_TP3[(hit->GetTrackNumber())-1] = " << Edep_TP3.back() << G4endl; } Hit_mult_TP3 = X_pos_TP3.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det TP3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det TP3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** TRUTH PLANE 4 **********************************// nHits = hits_TP4->entries(); Edep_TP4.clear(); for ( G4int h = 0 ; (h( hits_TP4->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; //Use to obtain exit KE of primary proton after passing through tracker and phantom if ( hit->GetIsPrimary() == true ){KE_TP4 = static_cast(hit->GetKE());} G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_TP4 += edep; th_TP4 = std::atan( (x1-x0) / (z1-z0) ); th_TP4 /= mrad; if( Track_no_TP4.size()>0 && (hit->GetTrackNumber()) == (Track_no_TP4.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_TP4[Track_no_TP4.size()-1] = hit->GetParticleName(); Track_no_TP4[Track_no_TP4.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_TP4[Track_no_TP4.size()-1] = hit->GetIsPrimary(); X_pos_TP4[Track_no_TP4.size()-1] = x1; Y_pos_TP4[Track_no_TP4.size()-1] = y1; Z_pos_TP4[Track_no_TP4.size()-1] = z1; Theta_TP4[Track_no_TP4.size()-1] = th_TP4; Edep_TP4[Track_no_TP4.size()-1] += edep; ParticleZ_TP4[Track_no_TP4.size()-1] = hit->GetParticleZ(); IsUp_TP4[Track_no_TP4.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP4 for particle: " << hit->GetParticleName() << ", Edep_TP4.back() = " << Edep_TP4.back() << G4endl; } else { ParticleName_TP4.push_back(hit->GetParticleName()); Track_no_TP4.push_back(hit->GetTrackNumber()); IsPrimaryParticle_TP4.push_back(hit->GetIsPrimary()); X_pos_TP4.push_back(x1); Y_pos_TP4.push_back(y1); Z_pos_TP4.push_back(z1); Theta_TP4.push_back(th_TP4); Edep_TP4.push_back(edep); ParticleZ_TP4.push_back(hit->GetParticleZ()); IsUp_TP4.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. TP4 for particle: " << hit->GetParticleName() << ", Edep_TP4[(hit->GetTrackNumber())-1] = " << Edep_TP4.back() << G4endl; } Hit_mult_TP4 = X_pos_TP4.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det TP4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det TP4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_TP4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** THIRD TRACKER MODULE **********************************// nHits = hits_x3->entries(); for ( G4int h = 0 ; (h( hits_x3->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; th_x3 = std::atan( (x1-x0) / (z1-z0) ); th_x3 /= mrad; if( Track_no_x3.size()>0 && (hit->GetTrackNumber()) == (Track_no_x3.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_x3[Track_no_x3.size()-1] = hit->GetParticleName(); Track_no_x3[Track_no_x3.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_x3[Track_no_x3.size()-1] = hit->GetIsPrimary(); X_pos_x3[Track_no_x3.size()-1] = x1; Y_pos_x3[Track_no_x3.size()-1] = y1; Z_pos_x3[Track_no_x3.size()-1] = z1; Theta_x3[Track_no_x3.size()-1] = th_x3; Edep_x3[Track_no_x3.size()-1] += edep; ParticleZ_x3[Track_no_x3.size()-1] = hit->GetParticleZ(); IsUp_x3[Track_no_x3.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_x3.push_back(hit->GetParticleName()); Track_no_x3.push_back(hit->GetTrackNumber()); IsPrimaryParticle_x3.push_back(hit->GetIsPrimary()); X_pos_x3.push_back(x1); Y_pos_x3.push_back(y1); Z_pos_x3.push_back(z1); Theta_x3.push_back(th_x3); Edep_x3.push_back(edep); ParticleZ_x3.push_back(hit->GetParticleZ()); IsUp_x3.push_back((y1 > 0) ? true : false); } Hit_mult_x3 = X_pos_x3.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det x3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det x3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_u3->entries(); Edep_u3.clear(); for ( G4int h = 0 ; (h( hits_u3->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; //Edep_u3 += edep; th_u3 = std::atan( (x1-x0) / (z1-z0) ); th_u3 /= mrad; //if( (hit->GetTrackNumber()) <= (ParticleName_u3.size()) ) //ensures each element represents a track, element [0] = track 1 if( Track_no_u3.size()>0 && (hit->GetTrackNumber()) == (Track_no_u3.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_u3[Track_no_u3.size()-1] = hit->GetParticleName(); Track_no_u3[Track_no_u3.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_u3[Track_no_u3.size()-1] = hit->GetIsPrimary(); X_pos_u3[Track_no_u3.size()-1] = x1; Y_pos_u3[Track_no_u3.size()-1] = y1; Z_pos_u3[Track_no_u3.size()-1] = z1; Theta_u3[Track_no_u3.size()-1] = th_u3; Edep_u3[Track_no_u3.size()-1] += edep; ParticleZ_u3[Track_no_u3.size()-1] = hit->GetParticleZ(); IsUp_u3[Track_no_u3.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_u3.push_back(hit->GetParticleName()); Track_no_u3.push_back(hit->GetTrackNumber()); IsPrimaryParticle_u3.push_back(hit->GetIsPrimary()); X_pos_u3.push_back(x1); Y_pos_u3.push_back(y1); Z_pos_u3.push_back(z1); Theta_u3.push_back(th_u3); Edep_u3.push_back(edep); ParticleZ_u3.push_back(hit->GetParticleZ()); IsUp_u3.push_back((y1 > 0) ? true : false); } Hit_mult_u3 = X_pos_u3.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det u3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det u3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_v3->entries(); Edep_v3.clear(); for ( G4int h = 0 ; (h( hits_v3->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_v3 += edep; th_v3 = std::atan( (x1-x0) / (z1-z0) ); th_v3 /= mrad; if( Track_no_v3.size()>0 && (hit->GetTrackNumber()) == (Track_no_v3.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_v3[Track_no_v3.size()-1] = hit->GetParticleName(); Track_no_v3[Track_no_v3.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_v3[Track_no_v3.size()-1] = hit->GetIsPrimary(); X_pos_v3[Track_no_v3.size()-1] = x1; Y_pos_v3[Track_no_v3.size()-1] = y1; Z_pos_v3[Track_no_v3.size()-1] = z1; Theta_v3[Track_no_v3.size()-1] = th_v3; Edep_v3[Track_no_v3.size()-1] += edep; ParticleZ_v3[Track_no_v3.size()-1] = hit->GetParticleZ(); IsUp_v3[Track_no_v3.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v3 for particle: " << hit->GetParticleName() << ", Edep_v3.back() = " << Edep_v3.back() << G4endl; } else { ParticleName_v3.push_back(hit->GetParticleName()); Track_no_v3.push_back(hit->GetTrackNumber()); IsPrimaryParticle_v3.push_back(hit->GetIsPrimary()); X_pos_v3.push_back(x1); Y_pos_v3.push_back(y1); Z_pos_v3.push_back(z1); Theta_v3.push_back(th_v3); Edep_v3.push_back(edep); ParticleZ_v3.push_back(hit->GetParticleZ()); IsUp_v3.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v3 for particle: " << hit->GetParticleName() << ", Edep_v3[(hit->GetTrackNumber())-1] = " << Edep_v3.back() << G4endl; } Hit_mult_v3 = X_pos_v3.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det v3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det v3 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v3.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } //********************************** FOURTH TRACKER MODULE **********************************// nHits = hits_x4->entries(); for ( G4int h = 0 ; (h( hits_x4->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; th_x4 = std::atan( (x1-x0) / (z1-z0) ); th_x4 /= mrad; if( Track_no_x4.size()>0 && (hit->GetTrackNumber()) == (Track_no_x4.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_x4[Track_no_x4.size()-1] = hit->GetParticleName(); Track_no_x4[Track_no_x4.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_x4[Track_no_x4.size()-1] = hit->GetIsPrimary(); X_pos_x4[Track_no_x4.size()-1] = x1; Y_pos_x4[Track_no_x4.size()-1] = y1; Z_pos_x4[Track_no_x4.size()-1] = z1; Theta_x4[Track_no_x4.size()-1] = th_x4; Edep_x4[Track_no_x4.size()-1] += edep; ParticleZ_x4[Track_no_x4.size()-1] = hit->GetParticleZ(); IsUp_x4[Track_no_x4.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_x4.push_back(hit->GetParticleName()); Track_no_x4.push_back(hit->GetTrackNumber()); IsPrimaryParticle_x4.push_back(hit->GetIsPrimary()); X_pos_x4.push_back(x1); Y_pos_x4.push_back(y1); Z_pos_x4.push_back(z1); Theta_x4.push_back(th_x4); Edep_x4.push_back(edep); ParticleZ_x4.push_back(hit->GetParticleZ()); IsUp_x4.push_back((y1 > 0) ? true : false); } Hit_mult_x4 = X_pos_x4.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det x4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det x4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_x4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_u4->entries(); Edep_u4.clear(); for ( G4int h = 0 ; (h( hits_u4->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; //Edep_u4 += edep; th_u4 = std::atan( (x1-x0) / (z1-z0) ); th_u4 /= mrad; //if( (hit->GetTrackNumber()) <= (ParticleName_u4.size()) ) //ensures each element represents a track, element [0] = track 1 if( Track_no_u4.size()>0 && (hit->GetTrackNumber()) == (Track_no_u4.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_u4[Track_no_u4.size()-1] = hit->GetParticleName(); Track_no_u4[Track_no_u4.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_u4[Track_no_u4.size()-1] = hit->GetIsPrimary(); X_pos_u4[Track_no_u4.size()-1] = x1; Y_pos_u4[Track_no_u4.size()-1] = y1; Z_pos_u4[Track_no_u4.size()-1] = z1; Theta_u4[Track_no_u4.size()-1] = th_u4; Edep_u4[Track_no_u4.size()-1] += edep; ParticleZ_u4[Track_no_u4.size()-1] = hit->GetParticleZ(); IsUp_u4[Track_no_u4.size()-1] = (y1 > 0) ? true : false; } else { ParticleName_u4.push_back(hit->GetParticleName()); Track_no_u4.push_back(hit->GetTrackNumber()); IsPrimaryParticle_u4.push_back(hit->GetIsPrimary()); X_pos_u4.push_back(x1); Y_pos_u4.push_back(y1); Z_pos_u4.push_back(z1); Theta_u4.push_back(th_u4); Edep_u4.push_back(edep); ParticleZ_u4.push_back(hit->GetParticleZ()); IsUp_u4.push_back((y1 > 0) ? true : false); } Hit_mult_u4 = X_pos_u4.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det u4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det u4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_u4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } nHits = hits_v4->entries(); Edep_v4.clear(); for ( G4int h = 0 ; (h( hits_v4->GetHit( h ) ); //Uncomment this line if you want to record only //primary energy depositions //if ( hit->GetIsPrimary() == false ) continue; //Use to obtain exit KE of primary proton after passing through tracker and phantom if ( hit->GetIsPrimary() == true ){KE_out = static_cast(hit->GetKE());} G4ThreeVector particle_pos = hit->GetPosition(); //We save energy in MeV Float_t edep = static_cast(hit->GetEdep()); edep /= MeV; x0 = static_cast(primPos.x()); y0 = static_cast(primPos.y()); z0 = static_cast(primPos.z()); x1 = static_cast(particle_pos.x()); y1 = static_cast(particle_pos.y()); z1 = static_cast(particle_pos.z()); //We save positions in mm (world coordinates) x0 /= mm; y0 /= mm; z0 /= mm; x1 /= mm; y1 /= mm; z1 /= mm; ///Edep_v4 += edep; th_v4 = std::atan( (x1-x0) / (z1-z0) ); th_v4 /= mrad; if( Track_no_v4.size()>0 && (hit->GetTrackNumber()) == (Track_no_v4.back()) ) //ensures each element represents a track, element [0] = track 1 { ParticleName_v4[Track_no_v4.size()-1] = hit->GetParticleName(); Track_no_v4[Track_no_v4.size()-1] = hit->GetTrackNumber(); IsPrimaryParticle_v4[Track_no_v4.size()-1] = hit->GetIsPrimary(); X_pos_v4[Track_no_v4.size()-1] = x1; Y_pos_v4[Track_no_v4.size()-1] = y1; Z_pos_v4[Track_no_v4.size()-1] = z1; Theta_v4[Track_no_v4.size()-1] = th_v4; Edep_v4[Track_no_v4.size()-1] += edep; ParticleZ_v4[Track_no_v4.size()-1] = hit->GetParticleZ(); IsUp_v4[Track_no_v4.size()-1] = (y1 > 0) ? true : false; //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v4 for particle: " << hit->GetParticleName() << ", Edep_v4.back() = " << Edep_v4.back() << G4endl; } else { ParticleName_v4.push_back(hit->GetParticleName()); Track_no_v4.push_back(hit->GetTrackNumber()); IsPrimaryParticle_v4.push_back(hit->GetIsPrimary()); X_pos_v4.push_back(x1); Y_pos_v4.push_back(y1); Z_pos_v4.push_back(z1); Theta_v4.push_back(th_v4); Edep_v4.push_back(edep); ParticleZ_v4.push_back(hit->GetParticleZ()); IsUp_v4.push_back((y1 > 0) ? true : false); //if(edep==0)G4cout << "ZERO ENERGY DEPOSIT in det. v4 for particle: " << hit->GetParticleName() << ", Edep_v4[(hit->GetTrackNumber())-1] = " << Edep_v4.back() << G4endl; } Hit_mult_v4 = X_pos_v4.size(); //G4int planeNum = hit->GetPlaneNumber(); /*if((hit->GetIsPrimary() == true)) {G4cout << "The position of the primary " << hit->GetParticleName() << " hit in det v4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;} else {G4cout << "The position of the secondary " << hit->GetParticleName() << " hit in det v4 = (" << x1 << "," << y1 << "," << z1 << ") mm, it's energy is: " << edep << " MeV, energy sum is: " << Edep_v4.back() << " MeV, it's track no. is: " << hit->GetTrackNumber() << G4endl;}*/ } } //******************************************************************************************// else { G4cerr << "Error: No hits collection passed to RootSaver for this event" << G4endl; } Truth_x_pos = static_cast( primPos.x() ); Truth_y_pos = static_cast( primPos.y() ); Truth_z_pos = static_cast( primPos.z() ); //Measure angle of the beam in xz plane measured from z+ direction // -pi= 0 ) ? +1 : -1; Float_t sign_x = ( primMom.x()>= 0 ) ? +1 : -1; TruthTheta = ( primMom.z() != 0 ) ? TMath::PiOver2()*sign_x*(1-sign_z)+std::atan( primMom.x()/primMom.z() ) : sign_x*TMath::PiOver2(); //beam perpendicular to z */ TruthTheta_x = std::atan( primMom.x()/primMom.z() ); TruthTheta_x /= mrad; TruthTheta_y = std::atan( primMom.y()/primMom.z() ); TruthTheta_y /= mrad; trackerTree->Fill(); } void HistoManager::CreateStripRTTree(const std::string& treeName) { stripRTTree = new TTree(treeName.data(),treeName.data()); stripRTTree->Branch( "EventNumber" , &stripRT_evtNumber ); stripRTTree->Branch( "LayerNumber", &stripRT_plane); stripRTTree->Branch( "Strip", &stripRT_strip); stripRTTree->Branch( "Signal", &stripRT_signal); } void HistoManager::AddStripRTEvent(const G4int event, const StripHitCollection* const hits) { G4cout << "HistoManager::AddStripRTEvent" << G4endl; stripRT_evtNumber = event; // start by resetting everything for(int l(0); l<24; l++){ for(int s(0); s<1024; s++){ stripRT_Signal[l][s] = NULL; } } stripRT_plane.clear(); stripRT_strip.clear(); stripRT_signal.clear(); // now we need to loop over the hits collection and sum up the energy deposited per strip per channels G4int nHits = hits->entries(); for ( G4int h = 0 ; (h( hits->GetHit( h ) ); stripRT_Signal[hit->GetPlaneNumber()][hit->GetStripNumber()] += (hit->GetEdep()/eV)/3.6; } // now loop and print out all strips with signal and amount for(int l(0); l<24; l++){ for(int s(0); s<1024; s++){ if(stripRT_Signal[l][s]>0) { G4cout << l << " " << s << " " << stripRT_Signal[l][s] << G4endl; stripRT_plane.push_back(l); stripRT_strip.push_back(s); stripRT_signal.push_back(stripRT_Signal[l][s]); } } } stripRTTree->Fill(); }