// // ******************************************************************** // * 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.hh,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...... #ifndef HistoManager_h #define HistoManager_h 1 #include "globals.hh" #include "vector" #include "cassert" #include "TFile.h" #include "TTree.h" #include "TBranch.h" #include "TH1F.h" #include "TH2F.h" #include "TProfile.h" #include "StripHit.hh" #include "StripDigi.hh" #include "G4Types.hh" #include //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... class HistoMessenger; const G4int MaxHisto = 9999; const G4int MaxHisto2D = 9999; // Structure to allow ease of use for histograms struct HistoInfo { TH1F* histo; G4bool exist; G4String Label; G4String Title; G4int Nbins; G4double Vmin ; G4double Vmax ; G4double Unit ; G4double Width; G4bool ascii; }; struct Histo2DInfo { TH2F* histo; G4bool exist; G4String Label; G4String Title; G4int NbinsX; G4int NbinsY; G4double Xmin ; G4double Xmax ; G4double Ymin ; G4double Ymax ; G4double UnitX ; G4double UnitY ; }; struct ProfileInfo { TProfile* profile; G4bool exist; G4String Label; G4String Title; G4int Nbins; G4double Vmin ; G4double Vmax ; G4double Unit ; G4double Width; G4bool ascii; }; struct TreeInfo{ TTree* tree; G4String Label; G4String Title; G4bool exist; double rndmSeed; double evtNumber; double x; double y; double z; double px; double py; double pz; double trackID; double parentID; bool isIncident; bool isExit; double kineticEnergy; double edep; double niel; double pdgid; double charge; double processType; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... class HistoManager { public: HistoManager(); ~HistoManager(); static HistoManager* GetHistoManager(); void SetFileName (const G4String& name) { fileName[0] = name;}; void SetFileType (const G4String& name) { fileType = name;}; void SetFileOption (const G4String& name) { fileOption = name;}; void book(); void save(); // functions to use TTrees void SetTree(const G4String& label, const G4String& title); void 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=-1); void AddBranch( TTree* input_tree, G4String title, std::vector data ); void AddBranch( TTree* input_tree, G4String title, std::vector data ); void AddBranch( TTree* input_tree, G4String title, std::vector data ); G4int GetTreeAddress(const G4String& label); //functions to use TH1F void SetHisto (const G4String& label, const G4String& title, G4int,G4double,G4double,const G4String& unit="none"); void FillHisto(const G4String& label, G4double e, G4double weight = 1.0); G4int GetHistoAddress(const G4String& label); // functions to use TH2F void SetHisto2D (const G4String& label, const G4String& title, G4int nbinsx, G4double vminx, G4double vmaxx, G4int nbinsy, G4double vminy, G4double vmaxy, const G4String& unitx="none", const G4String& unity="none"); void FillHisto2D(const G4String& label, G4double x, G4double y, G4double weight = 1.0); G4int GetHisto2DAddress(const G4String& label); // functions to use TProfile void SetProfile(const G4String& label, const G4String& title, G4int,G4double,G4double,const G4String& unit="none"); void FillProfile(const G4String& label, G4double e, G4double weight = 1.0); G4int GetProfileAddress(const G4String& label); // functions for HistoMessenger (not yet implimented) void Scale (const G4String&, G4double); void RemoveHisto (const G4String&); void PrintHisto (const G4String&); G4bool HistoExist (G4int id) {return histoInfo[id].exist;} G4double GetHistoUnit(G4int id) {return histoInfo[id].Unit;} G4double GetBinWidth (G4int id) {return histoInfo[id].Width;} G4bool Histo2DExist(G4int id) {return histo2DInfo[id].exist;} TFile* GetFile(){ return file;}; private: //static G4ThreadLocal static HistoManager* fHistoManager; G4String fileName[2]; G4String fileType; G4String fileOption; TFile* file; //static G4ThreadLocal G4int threadID; std::vector histoInfo; std::vector histo2DInfo; std::vector treeInfo; std::vector profileInfo; G4bool factoryOn; HistoMessenger* histoMessenger; private: void saveAscii(); //----------------------------- // Code to be able to save Michelas RT in my root files //---------------------------- public: virtual void CreateRTTree(const std::string& treeName = "rangetel"); virtual void AddRTEvent(const G4int event, const G4int pdgid, const G4int layerNumber, const G4int x, const G4int y, const G4int nElectrons); virtual void StoreRTInfo(); private: TTree* rtTree; G4int rtEvtNumber; G4int rtPDGID; G4int rtLayerNumber; G4int rtX; G4int rtY; G4int rtElectrons; std::vector evtsnumbers,pdgids,layernumbers,xs,ys,electrons; std::fstream binout; // ---------------------------------- // Code to impliment Jon's strips with Michelas RT // ---------------------------------- public: virtual void CreateStripRTTree(const std::string& treeName="StripRT"); virtual void AddStripRTEvent(const G4int event,const StripHitCollection* const hits); private: TTree* stripRTTree; G4int stripRT_evtNumber; G4int stripRT_Signal[24][1024]; std::vector stripRT_signal; std::vector stripRT_plane; std::vector stripRT_strip; //---------------------------- // Code from Jon to try and get his trees working in my root files //---------------------------- public: virtual void CreateTrackerTree(const std::string& treeName = "tracker", const int planes = 0, const int strips = 0 ); virtual void 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& primaryPos, const G4ThreeVector& primaryMom, const G4float KE_in, const G4float KE_out); private: TTree* trackerTree; // Pointer to the ROOT TTree Int_t nStrips; // Number of strips in each det. plane //*** TTree variables ***// Int_t Event_no; //Hit multiplicity for det. (no. of tracks (particles) / det. / event ) Int_t Hit_mult_x1; Int_t Hit_mult_u1; Int_t Hit_mult_v1; Int_t Hit_mult_x2; Int_t Hit_mult_u2; Int_t Hit_mult_v2; Int_t Hit_mult_TP1; Int_t Hit_mult_TP2; Int_t Hit_mult_TP3; Int_t Hit_mult_TP4; Int_t Hit_mult_x3; Int_t Hit_mult_u3; Int_t Hit_mult_v3; Int_t Hit_mult_x4; Int_t Hit_mult_u4; Int_t Hit_mult_v4; //Signal in electrons for each strip on each det. Float_t* Signal_x1; Float_t* Signal_u1; Float_t* Signal_v1; Float_t* Signal_x2; Float_t* Signal_u2; Float_t* Signal_v2; Float_t* Signal_TP1; Float_t* Signal_TP2; Float_t* Signal_TP3; Float_t* Signal_TP4; Float_t* Signal_x3; Float_t* Signal_u3; Float_t* Signal_v3; Float_t* Signal_x4; Float_t* Signal_u4; Float_t* Signal_v4; //strips fired / det. / event std::vector Strip_no_x1; std::vector Strip_no_u1; std::vector Strip_no_v1; std::vector Strip_no_x2; std::vector Strip_no_u2; std::vector Strip_no_v2; std::vector Strip_no_TP1; std::vector Strip_no_TP2; std::vector Strip_no_TP3; std::vector Strip_no_TP4; std::vector Strip_no_x3; std::vector Strip_no_u3; std::vector Strip_no_v3; std::vector Strip_no_x4; std::vector Strip_no_u4; std::vector Strip_no_v4; //total no. of strips fired / det. / event Int_t ClusterSize_x1; Int_t ClusterSize_u1; Int_t ClusterSize_v1; Int_t ClusterSize_x2; Int_t ClusterSize_u2; Int_t ClusterSize_v2; Int_t ClusterSize_TP1; Int_t ClusterSize_TP2; Int_t ClusterSize_TP3; Int_t ClusterSize_TP4; Int_t ClusterSize_x3; Int_t ClusterSize_u3; Int_t ClusterSize_v3; Int_t ClusterSize_x4; Int_t ClusterSize_u4; Int_t ClusterSize_v4; //total no. of det. fired / event Int_t Det_mult; //Position and angle of primary particle at origin Float_t Truth_x_pos; Float_t Truth_y_pos; Float_t Truth_z_pos; Float_t TruthTheta_x; Float_t TruthTheta_y; //Kinetic Energy of primary at origin Float_t KE_in; //Kinetic Energy of primary at exit //(only works partially, gets last track stored in G4TrackingManager usually electron, sometimes proton, // have ensured it only record information for particles exceeding 100 MeV for now.) Float_t KE_out; //For prelim pCT reconstruction (with perfect energy info) KE imediately before/after phantom is //required so that no energy loss in detector materials/compensator has to be calcualted Float_t KE_TP1; Float_t KE_TP2; Float_t KE_TP3; Float_t KE_TP4; //Sum of edep (MeV) in det. std::vector Edep_x1; std::vector Edep_u1; std::vector Edep_v1; std::vector Edep_x2; std::vector Edep_u2; std::vector Edep_v2; std::vector Edep_TP1; std::vector Edep_TP2; std::vector Edep_TP3; std::vector Edep_TP4; std::vector Edep_x3; std::vector Edep_u3; std::vector Edep_v3; std::vector Edep_x4; std::vector Edep_u4; std::vector Edep_v4; //Polar angle in the xz plane (measured from z-axis) std::vector Theta_x1; std::vector Theta_u1; std::vector Theta_v1; std::vector Theta_x2; std::vector Theta_u2; std::vector Theta_v2; std::vector Theta_TP1; std::vector Theta_TP2; std::vector Theta_TP3; std::vector Theta_TP4; std::vector Theta_x3; std::vector Theta_u3; std::vector Theta_v3; std::vector Theta_x4; std::vector Theta_u4; std::vector Theta_v4; //x position of particle in det. (must be double as float/Float_t dont work! (bug)) std::vector X_pos_x1; std::vector X_pos_u1; std::vector X_pos_v1; std::vector X_pos_x2; std::vector X_pos_u2; std::vector X_pos_v2; std::vector X_pos_TP1; std::vector X_pos_TP2; std::vector X_pos_TP3; std::vector X_pos_TP4; std::vector X_pos_x3; std::vector X_pos_u3; std::vector X_pos_v3; std::vector X_pos_x4; std::vector X_pos_u4; std::vector X_pos_v4; //y position of particle in det. (must be double as float/Float_t dont work! (bug)) std::vector Y_pos_x1; std::vector Y_pos_u1; std::vector Y_pos_v1; std::vector Y_pos_x2; std::vector Y_pos_u2; std::vector Y_pos_v2; std::vector Y_pos_TP1; std::vector Y_pos_TP2; std::vector Y_pos_TP3; std::vector Y_pos_TP4; std::vector Y_pos_x3; std::vector Y_pos_u3; std::vector Y_pos_v3; std::vector Y_pos_x4; std::vector Y_pos_u4; std::vector Y_pos_v4; //z position of particle in det. (must be double as float/Float_t dont work! (bug)) std::vector Z_pos_x1; std::vector Z_pos_u1; std::vector Z_pos_v1; std::vector Z_pos_x2; std::vector Z_pos_u2; std::vector Z_pos_v2; std::vector Z_pos_TP1; std::vector Z_pos_TP2; std::vector Z_pos_TP3; std::vector Z_pos_TP4; std::vector Z_pos_x3; std::vector Z_pos_u3; std::vector Z_pos_v3; std::vector Z_pos_x4; std::vector Z_pos_u4; std::vector Z_pos_v4; //Does hit in det. occur on left side? std::vector IsUp_x1; std::vector IsUp_u1; std::vector IsUp_v1; std::vector IsUp_x2; std::vector IsUp_u2; std::vector IsUp_v2; std::vector IsUp_TP1; std::vector IsUp_TP2; std::vector IsUp_TP3; std::vector IsUp_TP4; std::vector IsUp_x3; std::vector IsUp_u3; std::vector IsUp_v3; std::vector IsUp_x4; std::vector IsUp_u4; std::vector IsUp_v4; //Does hit in det. contain primary event? std::vector IsPrimaryParticle_x1; std::vector IsPrimaryParticle_u1; std::vector IsPrimaryParticle_v1; std::vector IsPrimaryParticle_x2; std::vector IsPrimaryParticle_u2; std::vector IsPrimaryParticle_v2; std::vector IsPrimaryParticle_TP1; std::vector IsPrimaryParticle_TP2; std::vector IsPrimaryParticle_TP3; std::vector IsPrimaryParticle_TP4; std::vector IsPrimaryParticle_x3; std::vector IsPrimaryParticle_u3; std::vector IsPrimaryParticle_v3; std::vector IsPrimaryParticle_x4; std::vector IsPrimaryParticle_u4; std::vector IsPrimaryParticle_v4; //Particle track no. (for counting number of particles / det. / event) std::vector Track_no_x1; std::vector Track_no_u1; std::vector Track_no_v1; std::vector Track_no_x2; std::vector Track_no_u2; std::vector Track_no_v2; std::vector Track_no_TP1; std::vector Track_no_TP2; std::vector Track_no_TP3; std::vector Track_no_TP4; std::vector Track_no_x3; std::vector Track_no_u3; std::vector Track_no_v3; std::vector Track_no_x4; std::vector Track_no_u4; std::vector Track_no_v4; //Particle name / det. / event (must be TString otherwise won't work), //At command line use e.g. tv__tree->Draw("hit_mult_x1","ParticleName_x1==\"e-\""); to select as condition std::vector ParticleName_x1; std::vector ParticleName_u1; std::vector ParticleName_v1; std::vector ParticleName_x2; std::vector ParticleName_u2; std::vector ParticleName_v2; std::vector ParticleName_TP1; std::vector ParticleName_TP2; std::vector ParticleName_TP3; std::vector ParticleName_TP4; std::vector ParticleName_x3; std::vector ParticleName_u3; std::vector ParticleName_v3; std::vector ParticleName_x4; std::vector ParticleName_u4; std::vector ParticleName_v4; //Particle Charge (Z) / det. / event (dE proportional to Z^2) std::vector ParticleZ_x1; std::vector ParticleZ_u1; std::vector ParticleZ_v1; std::vector ParticleZ_x2; std::vector ParticleZ_u2; std::vector ParticleZ_v2; std::vector ParticleZ_TP1; std::vector ParticleZ_TP2; std::vector ParticleZ_TP3; std::vector ParticleZ_TP4; std::vector ParticleZ_x3; std::vector ParticleZ_u3; std::vector ParticleZ_v3; std::vector ParticleZ_x4; std::vector ParticleZ_u4; std::vector ParticleZ_v4; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... #endif