composite_calorimeter/src/CCalEventAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCalEndOfEventAction.cc
0028 // Description: CCalEndOfEventAction provides User actions at end of event
0029 ///////////////////////////////////////////////////////////////////////////////
0030 #include "CCalEventAction.hh"
0031 #include "CCaloSD.hh"
0032 #include "CCalPrimaryGeneratorAction.hh"
0033 #include "CCalSteppingAction.hh"
0034 #include "CCalG4HitCollection.hh"
0035 #include "CCalG4Hit.hh"
0036 #include "CCaloOrganization.hh"
0037 #include "CCalSDList.hh"
0038
0039 #include "G4AnalysisManager.hh"
0040 #include "G4ios.hh"
0041 #include "G4SystemOfUnits.hh"
0042 #include "G4Event.hh"
0043 #include "G4SDManager.hh"
0044 #include "G4HCofThisEvent.hh"
0045 #include "G4RunManager.hh"
0046 #include "G4UserSteppingAction.hh"
0047 #include "G4UImanager.hh"
0048
0049 #include <iostream>
0050 #include <vector>
0051 #include <map>
0052
0053 //#define debug
0054 //#define ddebug
0055
0056
0057 CCalEventAction::CCalEventAction(CCalPrimaryGeneratorAction* pg,
0058                                  CCalSteppingAction* sa):
0059   isInitialized(false),fPrimaryGenerator(pg),
0060   fSteppingAction(sa),SDnames(nullptr),numberOfSD(0)
0061 {
0062 #ifdef debug
0063   G4cout << "Instantiate CCalEndOfEventAction" << G4endl;
0064 #endif
0065
0066   G4cout << "Get Calorimter organisation" << G4endl;
0067   theOrg = new CCaloOrganization;
0068 }
0069
0070 CCalEventAction::~CCalEventAction() {
0071   delete theOrg;
0072   delete[] SDnames;
0073 }
0074
0075 void CCalEventAction::initialize() {
0076
0077   isInitialized = true;
0078   numberOfSD = CCalSDList::getInstance()->getNumberOfCaloSD();
0079 #ifdef debug
0080   G4cout << "CCalEndOfEventAction look for " << numberOfSD
0081        << " calorimeter-like SD" << G4endl;
0082 #endif
0083   if (numberOfSD > 0) {
0084     G4int n = numberOfSD;
0085     n = std::min(n, 2);
0086     SDnames = new nameType[n];
0087   }
0088   for (G4int i=0; i<numberOfSD; ++i) {
0089     SDnames[i] = G4String(CCalSDList::getInstance()->getCaloSDName(i));
0090 #ifdef debug
0091     G4cout << "CCalEndOfEventAction: found SD " << i << " name "
0092          << SDnames[i] << G4endl;
0093 #endif
0094   }
0095 }
0096
0097 void CCalEventAction::BeginOfEventAction(const G4Event* evt) {
0098
0099   if (!isInitialized) { initialize(); }
0100   G4cout << " --- Begin of event: " << evt->GetEventID() << G4endl;
0101   /*
0102   if(15 == evt->GetEventID()) {
0103     G4UImanager * UImanager = G4UImanager::GetUIpointer();
0104     UImanager->ApplyCommand("/tracking/verbose 2");
0105   }
0106   */
0107 }
0108
0109 void CCalEventAction::EndOfEventAction(const G4Event* evt){
0110
0111   fSteppingAction->endOfEvent();
0112   //
0113   // Look for the Hit Collection
0114   //
0115   G4HCofThisEvent* allHC = evt->GetHCofThisEvent();
0116   if (allHC == 0) {
0117 #ifdef debug
0118     G4cout << "CCalEndOfEventAction: No Hit Collection in this event"
0119          << G4endl;
0120 #endif
0121     return;
0122   }
0123
0124   //
0125   // hits info
0126   //
0127
0128   //Now make summary
0129   G4float hcalE[28], ecalE[49], fullE=0., edec=0, edhc=0;
0130   G4int i = 0;
0131   for (i = 0; i < 28; i++) {hcalE[i]=0.;}
0132   for (i = 0; i < 49; i++) {ecalE[i]=0.;}
0133
0134   G4float* edep = new G4float[numberOfSD];
0135   for (i = 0; i < numberOfSD; ++i){
0136
0137     //
0138     // Look for the Hit Collection
0139     //
0140     edep[i] = 0;
0141     G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
0142
0143     if (caloHCid >= 0) {
0144       CCalG4HitCollection* theHC =
0145         (CCalG4HitCollection*) allHC->GetHC(caloHCid);
0146
0147       if (theHC != 0) {
0148
0149         G4int nentries = theHC->entries();
0150 #ifdef debug
0151         G4cout << " There are " << nentries << " hits in " << SDnames[i]
0152                << " :" << G4endl;
0153 #endif
0154
0155         if (nentries > 0) {
0156
0157           G4int j;
0158           for (j=0; j<nentries; j++){
0159 #ifdef ddebug
0160             G4cout << "Hit " << j;
0161 #endif
0162             CCalG4Hit* aHit =  (*theHC)[j];
0163             G4float En = aHit->getEnergyDeposit();
0164             G4int unitID = aHit->getUnitID();
0165             G4int id=-1;
0166             if (unitID > 0 && unitID < 29) {
0167               id = unitID - 1; // HCal
0168               hcalE[id] += En/GeV;
0169             } else {
0170               G4int i0 = unitID/4096;
0171               G4int i1 = (unitID/64)%64;
0172               G4int i2 = unitID%64;
0173               if (i0 == 1 && i1 < 8 && i2 < 8) {
0174                 id = i1*7 + i2; // ECal
0175                 ecalE[id] += En/GeV;
0176               }
0177             }
0178 #ifdef ddebug
0179             G4cout << " with Energy = " << En/MeV << " MeV in Unit " << unitID
0180                    << " " << id << G4endl;
0181 #endif
0182             fullE   += En/GeV;
0183             edep[i] += En/GeV;
0184           }
0185 #ifdef ddebug
0186           G4cout << " ===> Total Energy Deposit in this Calorimeter = "
0187                  << edep[i]*1000.0 << "  MeV " << G4endl;
0188 #endif
0189         }
0190       }
0191     }
0192     if (SDnames[i] == "HadronCalorimeter") {
0193       edhc = edep[i];
0194     } else if (SDnames[i] == "CrystalMatrix") {
0195       edec = edep[i];
0196     }
0197   }
0198
0199   delete[] edep;
0200
0201   G4ThreeVector pos = fPrimaryGenerator->GetParticlePosition();
0202   G4float ener = fPrimaryGenerator->GetParticleEnergy()/GeV;
0203   G4float x    = pos.x()/mm;
0204   G4float y    = pos.y()/mm;
0205   G4float z    = pos.z()/mm;
0206
0207   //Save results
0208   G4AnalysisManager* man = G4AnalysisManager::Instance();
0209   //1)
0210   static G4int IDenergy = -1;
0211   if (IDenergy < 0)
0212     IDenergy = man->GetH1Id("h4000");
0213   man->FillH1(IDenergy,fullE);
0214   //2)
0215 #ifdef debug
0216   G4double totalFilledEnergyHcal = 0.0;
0217 #endif
0218   static G4int IDhcalE = -1;
0219   if (IDhcalE < 0)
0220     IDhcalE = man->GetH1Id("h100");
0221   for (G4int j=0; j<28; j++) {
0222     man->FillH1(IDhcalE+j,hcalE[j]);
0223 #ifdef debug
0224     G4cout << "Fill Hcal histo " << j << " with " << hcalE[j] << G4endl;
0225     totalFilledEnergyHcal += hcalE[j];
0226 #endif
0227   }
0228 #ifdef debug
0229   G4cout << "CCalAnalysis::InsertEnergyHcal: Total filled Energy Hcal histo "
0230          << totalFilledEnergyHcal << G4endl;
0231   totalFilledEnergyEcal = 0.0;
0232 #endif
0233
0234   //3)
0235   static G4int IDecalE = -1;
0236   if (IDecalE < 0)
0237     IDecalE = man->GetH1Id("h200");
0238   for (G4int j=0; j<49; j++) {
0239     man->FillH1(IDecalE+j,ecalE[j]);
0240 #ifdef debug
0241     G4cout << "Fill Ecal histo " << j << " with " << ecalE[j] << G4endl;
0242     totalFilledEnergyEcal += ecalE[j];
0243 #endif
0244   }
0245 #ifdef debug
0246   G4cout << "CCalAnalysis::InsertEnergyEal: Total filled Energy Ecal histo "
0247          << totalFilledEnergyEcal << G4endl;
0248 #endif
0249   // 4)
0250   G4int counter=0;
0251   for (G4int j=0; j<28; j++)
0252     {
0253       man->FillNtupleFColumn(counter,hcalE[j]);
0254       counter++;
0255     }
0256   for (G4int j=0; j<49; j++)
0257     {
0258       man->FillNtupleFColumn(counter,ecalE[j]);
0259       counter++;
0260     }
0261   man->FillNtupleFColumn(counter,ener);
0262   man->FillNtupleFColumn(counter+1,x);
0263   man->FillNtupleFColumn(counter+2,y);
0264   man->FillNtupleFColumn(counter+3,z);
0265   man->FillNtupleFColumn(counter+4,fullE);
0266   man->FillNtupleFColumn(counter+5,edec);
0267   man->FillNtupleFColumn(counter+6,edhc);
0268   man->AddNtupleRow();
0269 #ifdef debug
0270   G4cout << "CCalAnalysis:: Fill Ntuple " << G4endl;
0271 #endif
0272
0273   // 5)
0274   static G4int IDtimeProfile = -1;
0275   if (IDtimeProfile < 0)
0276     IDtimeProfile = man->GetH1Id("h901");
0277   for (i = 0; i < numberOfSD; i++){
0278     G4int caloHCid = G4SDManager::GetSDMpointer()->GetCollectionID(SDnames[i]);
0279     if (caloHCid >= 0) {
0280       CCalG4HitCollection* theHC =
0281         (CCalG4HitCollection*) allHC->GetHC(caloHCid);
0282       if (theHC != 0) {
0283         G4int nentries = theHC->entries();
0284         if (nentries > 0) {
0285           for (G4int k=0; k<nentries; k++) {
0286             CCalG4Hit* aHit =  (*theHC)[k];
0287             man->FillH1(IDtimeProfile,aHit->getTimeSlice(),aHit->getEnergyDeposit()/GeV);
0288
0289 #ifdef debug
0290             G4cout << "CCalAnalysis:: Fill Time Profile with Hit " << k
0291                    << " Edeposit " << edep << " Gev" << G4endl;
0292 #endif
0293           }
0294         }
0295       }
0296     }
0297   }
0298 }
0299
0300