Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-13 23:27:39 

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