Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:18 

0001 // system include files
0002 #include <cmath>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008 
0009 //
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 
0017 #include "SimG4Core/Notification/interface/TrackInformation.h"
0018 #include "SimG4Core/Notification/interface/Observer.h"
0019 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0020 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0021 #include "SimG4Core/Notification/interface/EndOfRun.h"
0022 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0023 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0024 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0025 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0026 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0027 
0028 // to retreive hits
0029 #include "SimG4CMS/Forward/interface/BscG4Hit.h"
0030 #include "SimG4CMS/Forward/interface/BscNumberingScheme.h"
0031 #include "SimG4CMS/Forward/interface/BscG4HitCollection.h"
0032 
0033 // G4 stuff
0034 #include "G4SDManager.hh"
0035 #include "G4Step.hh"
0036 #include "G4Track.hh"
0037 #include "G4VProcess.hh"
0038 #include "G4HCofThisEvent.hh"
0039 #include "G4UserEventAction.hh"
0040 #include "G4TransportationManager.hh"
0041 #include "G4ProcessManager.hh"
0042 #include "G4VTouchable.hh"
0043 
0044 #include <CLHEP/Vector/ThreeVector.h>
0045 #include <CLHEP/Vector/LorentzVector.h>
0046 #include <CLHEP/Random/Randomize.h>
0047 #include <CLHEP/Units/SystemOfUnits.h>
0048 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0049 
0050 // ----------------------------------------------------------------
0051 // Includes needed for Root ntupling
0052 //
0053 #include "TROOT.h"
0054 #include "TFile.h"
0055 #include "TH1.h"
0056 #include "TH2.h"
0057 #include "TProfile.h"
0058 #include "TNtuple.h"
0059 #include "TRandom.h"
0060 #include "TLorentzVector.h"
0061 #include "TUnixSystem.h"
0062 #include "TSystem.h"
0063 #include "TMath.h"
0064 #include "TF1.h"
0065 #include "TObjArray.h"
0066 #include "TObjString.h"
0067 #include "TNamed.h"
0068 
0069 //#define EDM_ML_DEBUG
0070 
0071 class BscAnalysisHistManager : public TNamed {
0072 public:
0073   BscAnalysisHistManager(const TString& managername);
0074   ~BscAnalysisHistManager() override;
0075 
0076   TH1F* getHisto(Int_t Number);
0077   TH1F* getHisto(const TObjString& histname);
0078 
0079   TH2F* getHisto2(Int_t Number);
0080   TH2F* getHisto2(const TObjString& histname);
0081 
0082   void writeToFile(const TString& fOutputFile, const TString& fRecreateFile);
0083 
0084 private:
0085   void bookHistos();
0086   void storeWeights();
0087   void histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup);
0088   void histInit(
0089       const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup);
0090 
0091   const char* fTypeTitle;
0092   TObjArray* fHistArray;
0093   TObjArray* fHistNamesArray;
0094 };
0095 
0096 class BscTest : public SimWatcher,
0097                 public Observer<const BeginOfJob*>,
0098                 public Observer<const BeginOfRun*>,
0099                 public Observer<const EndOfRun*>,
0100                 public Observer<const BeginOfEvent*>,
0101                 public Observer<const BeginOfTrack*>,
0102                 public Observer<const G4Step*>,
0103                 public Observer<const EndOfTrack*>,
0104                 public Observer<const EndOfEvent*> {
0105 public:
0106   BscTest(const edm::ParameterSet& p);
0107   ~BscTest() override;
0108 
0109 private:
0110   // observer classes
0111   void update(const BeginOfJob* run) override;
0112   void update(const BeginOfRun* run) override;
0113   void update(const EndOfRun* run) override;
0114   void update(const BeginOfEvent* evt) override;
0115   void update(const BeginOfTrack* trk) override;
0116   void update(const G4Step* step) override;
0117   void update(const EndOfTrack* trk) override;
0118   void update(const EndOfEvent* evt) override;
0119 
0120   int iev;
0121   int itrk;
0122   G4double entot0, tracklength0;
0123 
0124   // Utilities to get detector levels during a step
0125 
0126   int detLevels(const G4VTouchable*) const;
0127   G4String detName(const G4VTouchable*, int, int) const;
0128   void detectorLevel(const G4VTouchable*, int&, int*, G4String*) const;
0129 
0130   double rinCalo, zinCalo;
0131   int lastTrackID;
0132   int verbosity;
0133 
0134   // sumEnerDeposit - all deposited energy on all steps ;  sumStepl - length in steel !!!
0135   G4double sumEnerDeposit, sumStepl, sumStepc;
0136   // numofpart - # particles produced along primary track
0137   int numofpart;
0138   // last point of the primary track
0139   G4ThreeVector lastpo;
0140 
0141   // z:
0142   double z1, z2, z3, z4;
0143 
0144 private:
0145   Float_t bsceventarray[1];
0146   TNtuple* bsceventntuple;
0147   TFile bscOutputFile;
0148   int whichevent;
0149 
0150   BscAnalysisHistManager* TheHistManager;  //Histogram Manager of the analysis
0151   std::string fDataLabel;                  // Data type label
0152   std::string fOutputFile;                 // The output file name
0153   std::string fRecreateFile;               // Recreate the file flag, default="RECREATE"
0154 };
0155 
0156 enum ntbsc_elements { ntbsc_evt };
0157 
0158 //================================================================
0159 BscTest::BscTest(const edm::ParameterSet& p) {
0160   //constructor
0161   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("BscTest");
0162   verbosity = m_Anal.getParameter<int>("Verbosity");
0163   //verbosity    = 1;
0164 
0165   fDataLabel = m_Anal.getParameter<std::string>("FDataLabel");
0166   fOutputFile = m_Anal.getParameter<std::string>("FOutputFile");
0167   fRecreateFile = m_Anal.getParameter<std::string>("FRecreateFile");
0168 
0169   if (verbosity > 0) {
0170     edm::LogVerbatim("BscTest") << "============================================================================";
0171     edm::LogVerbatim("BscTest") << "BscTestconstructor :: Initialized as observer";
0172   }
0173   // Initialization:
0174 
0175   bsceventntuple = new TNtuple("NTbscevent", "NTbscevent", "evt");
0176   whichevent = 0;
0177   TheHistManager = new BscAnalysisHistManager(fDataLabel);
0178 
0179   if (verbosity > 0) {
0180     edm::LogVerbatim("BscTest") << "BscTest constructor :: Initialized BscAnalysisHistManager";
0181   }
0182 }
0183 
0184 BscTest::~BscTest() {
0185   //  delete UserNtuples;
0186 
0187   TFile bscOutputFile("newntbsc.root", "RECREATE");
0188   edm::LogVerbatim("BscTest") << "Bsc output root file has been created";
0189   bsceventntuple->Write();
0190   edm::LogVerbatim("BscTest") << ", written";
0191   bscOutputFile.Close();
0192   edm::LogVerbatim("BscTest") << ", closed";
0193   delete bsceventntuple;
0194   edm::LogVerbatim("BscTest") << ", and deleted";
0195 
0196   //------->while end
0197 
0198   // Write histograms to file
0199   TheHistManager->writeToFile(fOutputFile, fRecreateFile);
0200   if (verbosity > 0) {
0201     edm::LogVerbatim("BscTest") << std::endl << "BscTest Destructor  -------->  End of BscTest : ";
0202   }
0203 
0204   edm::LogVerbatim("BscTest") << "BscTest: End of process";
0205 }
0206 
0207 //================================================================
0208 // Histoes:
0209 //-----------------------------------------------------------------------------
0210 
0211 BscAnalysisHistManager::BscAnalysisHistManager(const TString& managername) {
0212   // The Constructor
0213 
0214   fTypeTitle = managername;
0215   fHistArray = new TObjArray();       // Array to store histos
0216   fHistNamesArray = new TObjArray();  // Array to store histos's names
0217 
0218   bookHistos();
0219 
0220   fHistArray->Compress();  // Removes empty space
0221   fHistNamesArray->Compress();
0222 }
0223 //-----------------------------------------------------------------------------
0224 
0225 BscAnalysisHistManager::~BscAnalysisHistManager() {
0226   // The Destructor
0227 
0228   if (fHistArray) {
0229     fHistArray->Delete();
0230     delete fHistArray;
0231   }
0232 
0233   if (fHistNamesArray) {
0234     fHistNamesArray->Delete();
0235     delete fHistNamesArray;
0236   }
0237 }
0238 //-----------------------------------------------------------------------------
0239 
0240 void BscAnalysisHistManager::bookHistos() {
0241   // at Start: (mm)
0242   histInit("TrackPhi", "Primary Phi", 100, 0., 360.);
0243   histInit("TrackTheta", "Primary Theta", 100, 0., 180.);
0244   histInit("TrackP", "Track XY position Z+ ", 80, -80., 80., 80, -80., 80.);
0245   histInit("TrackM", "Track XY position Z-", 80, -80., 80., 80, -80., 80.);
0246   histInit("DetIDs", "Track DetId - vs +", 16, -0.5, 15.5, 16, 15.5, 31.5);
0247 }
0248 
0249 //-----------------------------------------------------------------------------
0250 
0251 void BscAnalysisHistManager::writeToFile(const TString& fOutputFile, const TString& fRecreateFile) {
0252   //Write to file = fOutputFile
0253 
0254   edm::LogVerbatim("BscTest") << "================================================================";
0255   edm::LogVerbatim("BscTest") << " Write this Analysis to File " << fOutputFile;
0256   edm::LogVerbatim("BscTest") << "================================================================";
0257 
0258   TFile* file = new TFile(fOutputFile, fRecreateFile);
0259 
0260   fHistArray->Write();
0261   file->Close();
0262 }
0263 //-----------------------------------------------------------------------------
0264 
0265 void BscAnalysisHistManager::histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup) {
0266   // Add histograms and histograms names to the array
0267 
0268   char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
0269   strcpy(newtitle, title);
0270   strcat(newtitle, " (");
0271   strcat(newtitle, fTypeTitle);
0272   strcat(newtitle, ") ");
0273   fHistArray->AddLast((new TH1F(name, newtitle, nbinsx, xlow, xup)));
0274   fHistNamesArray->AddLast(new TObjString(name));
0275 }
0276 //-----------------------------------------------------------------------------
0277 
0278 void BscAnalysisHistManager::histInit(
0279     const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup) {
0280   // Add histograms and histograms names to the array
0281 
0282   char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
0283   strcpy(newtitle, title);
0284   strcat(newtitle, " (");
0285   strcat(newtitle, fTypeTitle);
0286   strcat(newtitle, ") ");
0287   fHistArray->AddLast((new TH2F(name, newtitle, nbinsx, xlow, xup, nbinsy, ylow, yup)));
0288   fHistNamesArray->AddLast(new TObjString(name));
0289 }
0290 //-----------------------------------------------------------------------------
0291 
0292 TH1F* BscAnalysisHistManager::getHisto(Int_t Number) {
0293   // Get a histogram from the array with index = Number
0294 
0295   if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
0296     return (TH1F*)(fHistArray->At(Number));
0297 
0298   } else {
0299     edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!getHisto!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0300     edm::LogVerbatim("BscTest") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
0301     edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0302 
0303     return (TH1F*)(fHistArray->At(0));
0304   }
0305 }
0306 //-----------------------------------------------------------------------------
0307 
0308 TH2F* BscAnalysisHistManager::getHisto2(Int_t Number) {
0309   // Get a histogram from the array with index = Number
0310 
0311   if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
0312     return (TH2F*)(fHistArray->At(Number));
0313 
0314   } else {
0315     edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!getHisto2!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0316     edm::LogVerbatim("BscTest") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
0317     edm::LogVerbatim("BscTest") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0318 
0319     return (TH2F*)(fHistArray->At(0));
0320   }
0321 }
0322 //-----------------------------------------------------------------------------
0323 
0324 TH1F* BscAnalysisHistManager::getHisto(const TObjString& histname) {
0325   // Get a histogram from the array with name = histname
0326 
0327   Int_t index = fHistNamesArray->IndexOf(&histname);
0328   return getHisto(index);
0329 }
0330 //-----------------------------------------------------------------------------
0331 
0332 TH2F* BscAnalysisHistManager::getHisto2(const TObjString& histname) {
0333   // Get a histogram from the array with name = histname
0334 
0335   Int_t index = fHistNamesArray->IndexOf(&histname);
0336   return getHisto2(index);
0337 }
0338 //-----------------------------------------------------------------------------
0339 
0340 void BscAnalysisHistManager::storeWeights() {
0341   // Add structure to each histogram to store the weights
0342 
0343   for (int i = 0; i < fHistArray->GetEntries(); i++) {
0344     ((TH1F*)(fHistArray->At(i)))->Sumw2();
0345   }
0346 }
0347 
0348 //==================================================================== per JOB
0349 void BscTest::update(const BeginOfJob* job) {
0350   //job
0351   edm::LogVerbatim("BscTest") << "BscTest:beggining of job";
0352   ;
0353 }
0354 
0355 //==================================================================== per RUN
0356 void BscTest::update(const BeginOfRun* run) {
0357   //run
0358 
0359   edm::LogVerbatim("BscTest") << std::endl << "BscTest:: Begining of Run";
0360 }
0361 
0362 void BscTest::update(const EndOfRun* run) { ; }
0363 
0364 //=================================================================== per EVENT
0365 void BscTest::update(const BeginOfEvent* evt) {
0366   iev = (*evt)()->GetEventID();
0367   if (verbosity > 0) {
0368     edm::LogVerbatim("BscTest") << "BscTest:update Event number = " << iev;
0369   }
0370   whichevent++;
0371 }
0372 
0373 //=================================================================== per Track
0374 void BscTest::update(const BeginOfTrack* trk) {
0375   itrk = (*trk)()->GetTrackID();
0376   if (verbosity > 1) {
0377     edm::LogVerbatim("BscTest") << "BscTest:update BeginOfTrack number = " << itrk;
0378   }
0379   if (itrk == 1) {
0380     sumEnerDeposit = 0.;
0381     numofpart = 0;
0382     sumStepl = 0.;
0383     sumStepc = 0.;
0384     tracklength0 = 0.;
0385   }
0386 }
0387 
0388 //=================================================================== per EndOfTrack
0389 void BscTest::update(const EndOfTrack* trk) {
0390   itrk = (*trk)()->GetTrackID();
0391   if (verbosity > 1) {
0392     edm::LogVerbatim("BscTest") << "BscTest:update EndOfTrack number = " << itrk;
0393   }
0394   if (itrk == 1) {
0395     G4double tracklength = (*trk)()->GetTrackLength();  // Accumulated track length
0396 
0397     TheHistManager->getHisto("SumEDep")->Fill(sumEnerDeposit);
0398     TheHistManager->getHisto("TrackL")->Fill(tracklength);
0399 
0400     // direction !!!
0401     G4ThreeVector vert_mom = (*trk)()->GetVertexMomentumDirection();
0402     G4ThreeVector vert_pos = (*trk)()->GetVertexPosition();  // vertex ,where this track was created
0403 
0404     // last step information
0405     const G4Step* aStep = (*trk)()->GetStep();
0406     G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0407     lastpo = preStepPoint->GetPosition();
0408   }
0409 }
0410 
0411 // ====================================================
0412 
0413 //=================================================================== each STEP
0414 void BscTest::update(const G4Step* aStep) {
0415   // ==========================================================================
0416 
0417   if (verbosity > 2) {
0418     G4int stepnumber = aStep->GetTrack()->GetCurrentStepNumber();
0419     edm::LogVerbatim("BscTest") << "BscTest:update Step number = " << stepnumber;
0420   }
0421   // track on aStep:                                                                                         !
0422   G4Track* theTrack = aStep->GetTrack();
0423   TrackInformation* trkInfo = dynamic_cast<TrackInformation*>(theTrack->GetUserInformation());
0424   if (trkInfo == nullptr) {
0425     edm::LogVerbatim("BscTest") << "BscTest on aStep: No trk info !!!! abort ";
0426   }
0427   G4int id = theTrack->GetTrackID();
0428   G4String particleType = theTrack->GetDefinition()->GetParticleName();  //   !!!
0429   G4int parentID = theTrack->GetParentID();                              //   !!!
0430   G4TrackStatus trackstatus = theTrack->GetTrackStatus();                //   !!!
0431   G4double tracklength = theTrack->GetTrackLength();                     // Accumulated track length
0432   G4ThreeVector trackmom = theTrack->GetMomentum();
0433   G4double entot = theTrack->GetTotalEnergy();  //   !!! deposited on step
0434   G4int curstepnumber = theTrack->GetCurrentStepNumber();
0435   G4double stepl = aStep->GetStepLength();
0436   G4double EnerDeposit = aStep->GetTotalEnergyDeposit();
0437   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0438   const G4ThreeVector& preposition = preStepPoint->GetPosition();
0439   G4ThreeVector prelocalpoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(preposition);
0440   G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
0441   const G4String& prename = currentPV->GetName();
0442 
0443   const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
0444   int pre_levels = detLevels(pre_touch);
0445   G4String name1[20];
0446   int copyno1[20];
0447   for (int i = 0; i < 20; ++i) {
0448     name1[i] = "";
0449     copyno1[i] = 0;
0450   }
0451   if (pre_levels > 0) {
0452     detectorLevel(pre_touch, pre_levels, copyno1, name1);
0453   }
0454 
0455   if (id == 1) {
0456     // on 1-st step:
0457     if (curstepnumber == 1) {
0458       entot0 = entot;
0459       //UserNtuples->fillg519(entot0,1.);
0460     }
0461 
0462     // on every step:
0463 
0464     // for Copper:
0465     if (prename == "SBST") {
0466       sumStepc += stepl;
0467       // =========
0468     }
0469     // for ststeel:
0470     //   if(prename == "SBSTs") {
0471     if (prename == "SBSTs") {
0472       sumStepl += stepl;
0473       // =========
0474     }
0475     // =========
0476     // =========
0477 
0478     // exclude last track point if it is in SD (MI was started their)
0479     if (trackstatus != 2) {
0480       // for SD:   Si Det.:   SISTATION:SIPLANE:(SIDETL+BOUNDDET        +SIDETR + CERAMDET)
0481       if (prename == "SIDETL" || prename == "SIDETR") {
0482         if (prename == "SIDETL") {
0483           //UserNtuples->fillg569(EnerDeposit,1.);
0484         }
0485         if (prename == "SIDETR") {
0486           //UserNtuples->fillg570(EnerDeposit,1.);
0487         }
0488 
0489         G4String posname = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
0490         if ((prename == "SIDETL" && posname != "SIDETL") || (prename == "SIDETR" && posname != "SIDETR")) {
0491           if (name1[2] == "SISTATION") {
0492             //UserNtuples->fillg539(copyno1[2],1.);
0493           }
0494           if (name1[3] == "SIPLANE") {
0495             //UserNtuples->fillg540(copyno1[3],1.);
0496           }
0497 
0498           if (prename == "SIDETL") {
0499             //UserNtuples->fillg541(EnerDeposit,1.);
0500             //UserNtuples->fillg561(numbcont,1.);
0501             if (copyno1[2] < 2) {
0502               //UserNtuples->fillg571(dx,1.);
0503             } else if (copyno1[2] < 3) {
0504               //UserNtuples->fillg563(dx,1.);
0505               if (copyno1[3] < 2) {
0506               } else if (copyno1[3] < 3) {
0507                 //UserNtuples->fillg572(dx,1.);
0508               } else if (copyno1[3] < 4) {
0509                 //UserNtuples->fillg573(dx,1.);
0510               } else if (copyno1[3] < 5) {
0511                 //UserNtuples->fillg574(dx,1.);
0512               } else if (copyno1[3] < 6) {
0513                 //UserNtuples->fillg575(dx,1.);
0514               } else if (copyno1[3] < 7) {
0515                 //UserNtuples->fillg576(dx,1.);
0516               } else if (copyno1[3] < 8) {
0517                 //UserNtuples->fillg577(dx,1.);
0518               } else if (copyno1[3] < 9) {
0519                 //UserNtuples->fillg578(dx,1.);
0520               } else if (copyno1[3] < 10) {
0521                 //UserNtuples->fillg579(dx,1.);
0522               }
0523             } else if (copyno1[2] < 4) {
0524               //UserNtuples->fillg565(dx,1.);
0525             } else if (copyno1[2] < 5) {
0526               //UserNtuples->fillg567(dx,1.);
0527             }
0528           }
0529           if (prename == "SIDETR") {
0530             //UserNtuples->fillg542(EnerDeposit,1.);
0531             //UserNtuples->fillg562(numbcont,1.);
0532             if (copyno1[2] < 2) {
0533               //UserNtuples->fillg581(dy,1.);
0534             } else if (copyno1[2] < 3) {
0535               //UserNtuples->fillg564(dy,1.);
0536               if (copyno1[3] < 2) {
0537               } else if (copyno1[3] < 3) {
0538                 //UserNtuples->fillg582(dy,1.);
0539               } else if (copyno1[3] < 4) {
0540                 //UserNtuples->fillg583(dy,1.);
0541               } else if (copyno1[3] < 5) {
0542                 //UserNtuples->fillg584(dy,1.);
0543               } else if (copyno1[3] < 6) {
0544                 //UserNtuples->fillg585(dy,1.);
0545               } else if (copyno1[3] < 7) {
0546                 //UserNtuples->fillg586(dy,1.);
0547               } else if (copyno1[3] < 8) {
0548                 //UserNtuples->fillg587(dy,1.);
0549               } else if (copyno1[3] < 9) {
0550                 //UserNtuples->fillg588(dy,1.);
0551               } else if (copyno1[3] < 10) {
0552                 //UserNtuples->fillg589(dy,1.);
0553               }
0554             } else if (copyno1[2] < 4) {
0555               //UserNtuples->fillg566(dy,1.);
0556             } else if (copyno1[2] < 5) {
0557               //UserNtuples->fillg568(dy,1.);
0558             }
0559           }
0560         }
0561       }
0562       // end of prenames SIDETL // SIDETR
0563     }
0564     // end of trackstatus != 2
0565 
0566     sumEnerDeposit += EnerDeposit;
0567     if (trackstatus == 2) {
0568       // primary track length
0569       //      //UserNtuples->fillg508(tracklength,1.);
0570       tracklength0 = tracklength;
0571     }
0572   }
0573   // end of primary track !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0574 
0575   if (parentID == 1 && curstepnumber == 1) {
0576     // particles deposit their energy along primary track
0577     numofpart += 1;
0578     if (prename == "SBST") {
0579       //UserNtuples->filld225(vert_pos.y(),vert_pos.z(),1.);
0580     }
0581     if (prename == "SBSTs") {
0582       //UserNtuples->filld226(vert_pos.y(),vert_pos.z(),1.);
0583     }
0584   }
0585 }
0586 // ==========================================================================
0587 // ==========================================================================
0588 int BscTest::detLevels(const G4VTouchable* touch) const {
0589   //Return number of levels
0590   if (touch)
0591     return ((touch->GetHistoryDepth()) + 1);
0592   else
0593     return 0;
0594 }
0595 // ==========================================================================
0596 
0597 G4String BscTest::detName(const G4VTouchable* touch, int level, int currentlevel) const {
0598   //Go down to current level
0599   if (level > 0 && level >= currentlevel) {
0600     int ii = level - currentlevel;
0601     return touch->GetVolume(ii)->GetName();
0602   } else {
0603     return "NotFound";
0604   }
0605 }
0606 
0607 void BscTest::detectorLevel(const G4VTouchable* touch, int& level, int* copyno, G4String* name) const {
0608   //Get name and copy numbers
0609   if (level > 0) {
0610     for (int ii = 0; ii < level; ii++) {
0611       int i = level - ii - 1;
0612       G4VPhysicalVolume* pv = touch->GetVolume(i);
0613       if (pv != nullptr)
0614         name[ii] = pv->GetName();
0615       else
0616         name[ii] = "Unknown";
0617       copyno[ii] = touch->GetReplicaNumber(i);
0618     }
0619   }
0620 }
0621 // ==========================================================================
0622 
0623 //===================================================================   End Of Event
0624 void BscTest::update(const EndOfEvent* evt) {
0625   // ==========================================================================
0626 
0627   if (verbosity > 1) {
0628     iev = (*evt)()->GetEventID();
0629     edm::LogVerbatim("BscTest") << "BscTest:update EndOfEvent = " << iev;
0630   }
0631   // Fill-in ntuple
0632   bsceventarray[ntbsc_evt] = (float)whichevent;
0633 
0634   //
0635   int trackID = 0;
0636   G4PrimaryParticle* thePrim = nullptr;
0637 
0638   // prim.vertex:
0639   G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0640   if (nvertex != 1)
0641     edm::LogVerbatim("BscTest") << "BscTest: My warning: NumberOfPrimaryVertex != 1  -->  = " << nvertex;
0642 
0643   for (int i = 0; i < nvertex; i++) {
0644     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0645     if (avertex == nullptr) {
0646       edm::LogVerbatim("BscTest") << "BscTest  End Of Event ERR: pointer to vertex = 0";
0647       continue;
0648     }
0649     G4int npart = avertex->GetNumberOfParticle();
0650     if (npart != 1)
0651       edm::LogVerbatim("BscTest") << "BscTest: My warning: NumberOfPrimaryPart != 1  -->  = " << npart;
0652     if (npart == 0)
0653       edm::LogVerbatim("BscTest") << "BscTest End Of Event ERR: no NumberOfParticle";
0654 
0655     if (thePrim == nullptr)
0656       thePrim = avertex->GetPrimary(trackID);
0657 
0658     if (thePrim != nullptr) {
0659       // primary vertex:
0660       G4double vx = 0., vy = 0., vz = 0.;
0661       vx = avertex->GetX0();
0662       vy = avertex->GetY0();
0663       vz = avertex->GetZ0();
0664       //UserNtuples->fillh01(vx);
0665       //UserNtuples->fillh02(vy);
0666       //UserNtuples->fillh03(vz);
0667       TheHistManager->getHisto("VtxX")->Fill(vx);
0668       TheHistManager->getHisto("VtxY")->Fill(vy);
0669       TheHistManager->getHisto("VtxZ")->Fill(vz);
0670     }
0671   }
0672   // prim.vertex loop end
0673 
0674   //=========================== thePrim != 0 ================================================================================
0675   if (thePrim != nullptr) {
0676     //
0677     // number of secondary particles deposited their energy along primary track
0678     //UserNtuples->fillg518(numofpart,1.);
0679     if (lastpo.z() < z4 && lastpo.perp() < 100.) {
0680       //UserNtuples->fillg536(numofpart,1.);
0681     }
0682     //
0683 
0684     // direction !!!
0685     G4ThreeVector mom = thePrim->GetMomentum();
0686 
0687     double phi = atan2(mom.y(), mom.x());
0688     if (phi < 0.)
0689       phi += twopi;
0690     double phigrad = phi * 180. / pi;
0691 
0692     double th = mom.theta();
0693     double eta = -log(tan(th / 2));
0694     TheHistManager->getHisto("PrimaryEta")->Fill(eta);
0695     TheHistManager->getHisto("PrimaryPhigrad")->Fill(phigrad);
0696     TheHistManager->getHisto("PrimaryTh")->Fill(th * 180. / pi);
0697 
0698     TheHistManager->getHisto("PrimaryLastpoZ")->Fill(lastpo.z());
0699     if (lastpo.z() < z4) {
0700       TheHistManager->getHisto("PrimaryLastpoX")->Fill(lastpo.x());
0701       TheHistManager->getHisto("PrimaryLastpoY")->Fill(lastpo.y());
0702     }
0703     if (numofpart > 4) {
0704       TheHistManager->getHisto("XLastpoNumofpart")->Fill(lastpo.x());
0705       TheHistManager->getHisto("YLastpoNumofpart")->Fill(lastpo.y());
0706     }
0707 
0708     // ==========================================================================
0709 
0710     // hit map for Bsc
0711     // ==================================
0712 
0713     std::map<int, float, std::less<int> > themap;
0714     std::map<int, float, std::less<int> > themap1;
0715 
0716     std::map<int, float, std::less<int> > themapxy;
0717     std::map<int, float, std::less<int> > themapz;
0718     // access to the G4 hit collections:  -----> this work OK:
0719 #ifdef EDM_ML_DEBUG
0720     edm::LogVerbatim("BscTest") << "1";
0721 #endif
0722     G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0723 #ifdef EDM_ML_DEBUG
0724     edm::LogVerbatim("BscTest") << "2";
0725 #endif
0726     if (verbosity > 0) {
0727       edm::LogVerbatim("BscTest") << "BscTest:  accessed all HC";
0728       ;
0729     }
0730     int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("BSCHits");
0731 
0732     BscG4HitCollection* theCAFI = (BscG4HitCollection*)allHC->GetHC(CAFIid);
0733     if (verbosity > 0) {
0734       edm::LogVerbatim("BscTest") << "BscTest: theCAFI->entries = " << theCAFI->entries();
0735     }
0736     int varia;  // = 0 -all; =1 - MI; =2 - noMI
0737     //varia = 0;
0738     if (lastpo.z() < z4) {
0739       varia = 1;
0740     } else {
0741       varia = 2;
0742     }  // no MI end:
0743     int nhits = theCAFI->entries();
0744     for (int j = 0; j < nhits; j++) {
0745       BscG4Hit* aHit = (*theCAFI)[j];
0746       const CLHEP::Hep3Vector& hitPoint = aHit->getEntry();
0747       double zz = hitPoint.z();
0748       TheHistManager->getHisto("zHits")->Fill(zz);
0749       if (tracklength0 > 8300.)
0750         TheHistManager->getHisto("zHitsTrLoLe")->Fill(zz);
0751     }
0752 
0753     if (varia == 2) {
0754       //int nhit11 = 0, nhit12 = 0, nhit13 = 0;
0755       double totallosenergy = 0.;
0756       for (int j = 0; j < nhits; j++) {
0757         BscG4Hit* aHit = (*theCAFI)[j];
0758 
0759         const CLHEP::Hep3Vector& hitEntryLocalPoint = aHit->getEntryLocalP();
0760         const CLHEP::Hep3Vector& hitExitLocalPoint = aHit->getExitLocalP();
0761         const CLHEP::Hep3Vector& hitPoint = aHit->getEntry();
0762         //int trackIDhit = aHit->getTrackID();
0763         unsigned int unitID = aHit->getUnitID();
0764         double losenergy = aHit->getEnergyLoss();
0765 
0766         double zz = hitPoint.z();
0767 
0768         TheHistManager->getHisto("zHitsnoMI")->Fill(zz);
0769 
0770         if (verbosity > 2) {
0771           edm::LogVerbatim("BscTest") << "BscTest:zHits = " << zz;
0772         }
0773 
0774         themap[unitID] += losenergy;
0775         totallosenergy += losenergy;
0776 
0777         int zside;
0778         //int sector;
0779         BscNumberingScheme::unpackBscIndex(unitID);
0780         zside = (unitID & 32) >> 5;
0781         //sector = (unitID & 7);
0782 
0783         //
0784         //=======================================
0785         G4ThreeVector middle = (hitExitLocalPoint + hitEntryLocalPoint) / 2.;
0786         themapz[unitID] = hitPoint.z() + middle.z();
0787         //=======================================
0788         // Y
0789         if (zside == 1) {
0790           //UserNtuples->fillg24(losenergy,1.);
0791           if (losenergy > 0.00003) {
0792             themap1[unitID] += 1.;
0793           }
0794         }
0795         //X
0796         else if (zside == 2) {
0797           //UserNtuples->fillg25(losenergy,1.);
0798           if (losenergy > 0.00005) {
0799             themap1[unitID] += 1.;
0800           }
0801         }
0802         //     }
0803         //
0804         /*
0805         if (sector == 1) {
0806           nhit11 += 1;
0807           //UserNtuples->fillg33(rr,1.);
0808           //UserNtuples->fillg11(yy,1.);
0809         }
0810         if (sector == 2) {
0811           nhit12 += 1;
0812           //UserNtuples->fillg34(rr,1.);
0813           //UserNtuples->fillg86(yy,1.);
0814         }
0815         if (sector == 3) {
0816           nhit13 += 1;
0817           //UserNtuples->fillg35(rr,1.);
0818           //UserNtuples->fillg87(yy,1.);
0819         }
0820 
0821         if (lastpo.z() < z4 && lastpo.perp() < 120.) {
0822           // MIonly:
0823           //UserNtuples->fillg16(lastpo.z(),1.);
0824           //UserNtuples->fillg18(zz,1.);
0825           //                                                                     Station I
0826           if (zz < z2) {
0827             //UserNtuples->fillg54(dx,1.);
0828             //UserNtuples->fillg55(dy,1.);
0829           }
0830           //                                                                     Station II
0831           if (zz < z3 && zz > z2) {
0832             //UserNtuples->fillg50(dx,1.);
0833             //UserNtuples->fillg51(dy,1.);
0834           }
0835           //                                                                     Station III
0836           if (zz < z4 && zz > z3) {
0837             //UserNtuples->fillg64(dx,1.);
0838             //UserNtuples->fillg65(dy,1.);
0839             //UserNtuples->filld209(xx,yy,1.);
0840           }
0841         } else {
0842           // no MIonly:
0843           //UserNtuples->fillg17(lastpo.z(),1.);
0844           //UserNtuples->fillg19(zz,1.);
0845           //UserNtuples->fillg74(incidentEnergyHit,1.);
0846           //UserNtuples->fillg75(float(trackIDhit),1.);
0847           //                                                                     Station I
0848           if (zz < z2) {
0849             //UserNtuples->fillg56(dx,1.);
0850             //UserNtuples->fillg57(dy,1.);
0851             //UserNtuples->fillg20(numofpart,1.);
0852             //UserNtuples->fillg21(sumEnerDeposit,1.);
0853             if (zside == 1) {
0854               //UserNtuples->fillg26(losenergy,1.);
0855             }
0856             if (zside == 2) {
0857               //UserNtuples->fillg76(losenergy,1.);
0858             }
0859             if (trackIDhit == 1) {
0860               //UserNtuples->fillg70(dx,1.);
0861               //UserNtuples->fillg71(incidentEnergyHit,1.);
0862               //UserNtuples->fillg79(losenergy,1.);
0863             } else {
0864               //UserNtuples->fillg82(dx,1.);
0865             }
0866           }
0867           //                                                                     Station II
0868           if (zz < z3 && zz > z2) {
0869             //UserNtuples->fillg52(dx,1.);
0870             //UserNtuples->fillg53(dy,1.);
0871             //UserNtuples->fillg22(numofpart,1.);
0872             //UserNtuples->fillg23(sumEnerDeposit,1.);
0873             //UserNtuples->fillg80(incidentEnergyHit,1.);
0874             //UserNtuples->fillg81(float(trackIDhit),1.);
0875             if (zside == 1) {
0876               //UserNtuples->fillg27(losenergy,1.);
0877             }
0878             if (zside == 2) {
0879               //UserNtuples->fillg77(losenergy,1.);
0880             }
0881             if (trackIDhit == 1) {
0882               //UserNtuples->fillg72(dx,1.);
0883               //UserNtuples->fillg73(incidentEnergyHit,1.);
0884             } else {
0885               //UserNtuples->fillg83(dx,1.);
0886             }
0887           }
0888           //                                                                     Station III
0889           if (zz < z4 && zz > z3) {
0890             if (zside == 1) {
0891               //UserNtuples->fillg28(losenergy,1.);
0892             }
0893             if (zside == 2) {
0894               //UserNtuples->fillg78(losenergy,1.);
0895             }
0896           }
0897         }
0898         */
0899       }  // MIonly or noMIonly ENDED
0900       if (totallosenergy == 0.0) {
0901         edm::LogVerbatim("BscTest") << "BscTest:     number of hits = " << theCAFI->entries();
0902         for (int j = 0; j < nhits; j++) {
0903           BscG4Hit* aHit = (*theCAFI)[j];
0904           double losenergy = aHit->getEnergyLoss();
0905           edm::LogVerbatim("BscTest") << " j hits = " << j << "losenergy = " << losenergy;
0906         }
0907       }
0908       //   FIBRE Hit collected analysis
0909       /*
0910       double totalEnergy = 0.;
0911       int nhitsX = 0, nhitsY = 0, nsumhit = 0;
0912       for (int sector = 1; sector < 4; sector++) {
0913         int nhitsecX = 0, nhitsecY = 0;
0914         for (int zmodule = 1; zmodule < 11; zmodule++) {
0915           for (int zside = 1; zside < 3; zside++) {
0916             int det = 1;  // nhit = 0;
0917             //  int sScale = 20;
0918             int index = BscNumberingScheme::packBscIndex(det, zside, sector);
0919             double theTotalEnergy = themap[index];
0920             //   X planes
0921             if (zside < 2) {
0922               //UserNtuples->fillg47(theTotalEnergy,1.);
0923               if (theTotalEnergy > 0.00003) {
0924                 nhitsX += 1;
0925                 //      nhitsecX += themap1[index];
0926                 //      nhit=1;
0927               }
0928             }
0929             //   Y planes
0930             else {
0931               //UserNtuples->fillg49(theTotalEnergy,1.);
0932               if (theTotalEnergy > 0.00005) {
0933                 nhitsY += 1;
0934                 //      nhitsecY += themap1[index];
0935                 //      nhit=1;
0936               }
0937             }
0938 
0939             totalEnergy += themap[index];
0940           }  // for
0941         }    // for
0942              //UserNtuples->fillg39(nhitsecY,1.);
0943         if (nhitsecX > 10 || nhitsecY > 10) {
0944           nsumhit += 1;
0945           //UserNtuples->fillp213(float(sector),float(1.),1.);
0946         } else {  //UserNtuples->fillp213(float(sector),float(0.),1.);
0947         }
0948       }  // for
0949 
0950       if (nsumhit >= 2) {  //UserNtuples->fillp212(vy,float(1.),1.);
0951       } else {             //UserNtuples->fillp212(vy,float(0.),1.);
0952       }
0953       */
0954     }  // MI or no MI or all  - end
0955   }    // primary end
0956 
0957   if (verbosity > 0) {
0958     edm::LogVerbatim("BscTest") << "BscTest:  END OF Event " << (*evt)()->GetEventID();
0959   }
0960 }
0961 
0962 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0963 #include "FWCore/PluginManager/interface/ModuleDef.h"
0964 
0965 DEFINE_SIMWATCHER(BscTest);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!