File indexing completed on 2025-01-30 09:19:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 #include "TFile.h"
0052 #include "TSystem.h"
0053 #include "TTree.h"
0054 #include "TH1.h"
0055
0056 #include "Event.hh"
0057 #include "CalorimeterHit.hh"
0058
0059 int main(int argc, char** argv)
0060 {
0061
0062 TSystem ts;
0063 gSystem->Load("libCaTSClassesDict");
0064 if(argc < 3)
0065 {
0066 G4cout
0067 << "Program requires 2 arguments: name of input file, name of output file"
0068 << G4endl;
0069 exit(1);
0070 }
0071 TFile* outfile = new TFile(argv[2], "RECREATE");
0072 TFile fo(argv[1]);
0073 fo.GetListOfKeys()->Print();
0074 Event* event = new Event();
0075 TTree* Tevt = (TTree*) fo.Get("Events");
0076 Tevt->SetBranchAddress("event.", &event);
0077 TBranch* fevtbranch = Tevt->GetBranch("event.");
0078 Int_t nevent = fevtbranch->GetEntries();
0079 std::cout << " Nr. of Events: " << nevent << std::endl;
0080 double max = 0.;
0081 double min = 1000000;
0082 std::string CollectionName = argv[3];
0083 CollectionName = CollectionName + "_Calorimeter_HC";
0084 for(Int_t i = 0; i < nevent; i++)
0085 {
0086 fevtbranch->GetEntry(i);
0087 auto* hcmap = event->GetHCMap();
0088 for(const auto& ele : *hcmap)
0089 {
0090 if(ele.first.compare(CollectionName) == 0)
0091 {
0092 auto hits = ele.second;
0093 G4int NbHits = hits.size();
0094 for(G4int ii = 0; ii < NbHits; ii++)
0095 {
0096 CalorimeterHit* drcaloHit =
0097 dynamic_cast<CalorimeterHit*>(hits.at(ii));
0098 const double ed = drcaloHit->GetEdep();
0099 if(ed > max)
0100 max = ed;
0101 if(ed < min)
0102 min = ed;
0103 }
0104 }
0105 }
0106 }
0107 outfile->cd();
0108 TH1F* hedep =
0109 new TH1F("energy", "edep", 100, min - 0.1 * min, max + 0.1 * max);
0110 for(Int_t i = 0; i < nevent; i++)
0111 {
0112 fevtbranch->GetEntry(i);
0113 auto* hcmap = event->GetHCMap();
0114 for(const auto& ele : *hcmap)
0115 {
0116 if(ele.first.compare(CollectionName) == 0)
0117 {
0118 auto hits = ele.second;
0119 G4int NbHits = hits.size();
0120 for(G4int ii = 0; ii < NbHits; ii++)
0121 {
0122 CalorimeterHit* drcaloHit =
0123 dynamic_cast<CalorimeterHit*>(hits.at(ii));
0124 hedep->Fill(drcaloHit->GetEdep());
0125 }
0126 }
0127 }
0128 }
0129
0130 outfile->Write();
0131 }