Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-01 02:03:30 

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