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 <iostream>
0003 #include <map>
0004 #include <string>
0005 //
0006 // user include files
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/PluginManager/interface/ModuleDef.h"
0014 //
0015 #include "SimG4Core/Notification/interface/TrackInformation.h"
0016 #include "SimG4Core/Notification/interface/Observer.h"
0017 #include "SimG4Core/Notification/interface/BeginOfJob.h"
0018 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0019 #include "SimG4Core/Notification/interface/EndOfRun.h"
0020 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0021 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0022 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0023 #include "SimG4Core/Notification/interface/EndOfTrack.h"
0024 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0025 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0026 // to retreive hits
0027 #include "SimG4CMS/FP420/interface/FP420NumberingScheme.h"
0028 #include "SimG4CMS/FP420/interface/FP420G4HitCollection.h"
0029 // G4 stuff
0030 #include "G4HCofThisEvent.hh"
0031 #include "G4ProcessManager.hh"
0032 #include "G4SDManager.hh"
0033 #include "G4Step.hh"
0034 #include "G4Track.hh"
0035 #include "G4TransportationManager.hh"
0036 #include "G4UserEventAction.hh"
0037 #include "G4VProcess.hh"
0038 #include "G4VTouchable.hh"
0039 
0040 #include <CLHEP/Random/Randomize.h>
0041 #include <CLHEP/Units/GlobalSystemOfUnits.h>
0042 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0043 
0044 //================================================================
0045 // Root stuff
0046 #include "TROOT.h"
0047 #include "TFile.h"
0048 #include "TH1.h"
0049 #include "TH2.h"
0050 #include "TProfile.h"
0051 #include "TNtuple.h"
0052 #include "TRandom.h"
0053 #include "TLorentzVector.h"
0054 #include "TUnixSystem.h"
0055 #include "TSystem.h"
0056 #include "TMath.h"
0057 #include "TF1.h"
0058 #include "TObjArray.h"
0059 #include "TObjString.h"
0060 #include "TNamed.h"
0061 
0062 class Fp420AnalysisHistManager : public TNamed {
0063 public:
0064   Fp420AnalysisHistManager(const TString& managername);
0065   ~Fp420AnalysisHistManager() override;
0066 
0067   TH1F* getHisto(Int_t Number);
0068   TH1F* getHisto(const TObjString& histname);
0069 
0070   TH2F* getHisto2(Int_t Number);
0071   TH2F* getHisto2(const TObjString& histname);
0072 
0073   void writeToFile(const TString& fOutputFile, const TString& fRecreateFile);
0074 
0075 private:
0076   void bookHistos();
0077   void storeWeights();
0078   void histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup);
0079   void histInit(
0080       const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup);
0081 
0082   const char* fTypeTitle;
0083   TObjArray* fHistArray;
0084   TObjArray* fHistNamesArray;
0085 };
0086 
0087 class FP420Test : public SimWatcher,
0088                   public Observer<const BeginOfJob*>,
0089                   public Observer<const BeginOfRun*>,
0090                   public Observer<const EndOfRun*>,
0091                   public Observer<const BeginOfEvent*>,
0092                   public Observer<const BeginOfTrack*>,
0093                   public Observer<const G4Step*>,
0094                   public Observer<const EndOfTrack*>,
0095                   public Observer<const EndOfEvent*> {
0096 public:
0097   FP420Test(const edm::ParameterSet& p);
0098   ~FP420Test() override;
0099 
0100 private:
0101   // observer classes
0102   void update(const BeginOfJob* run) override;
0103   void update(const BeginOfRun* run) override;
0104   void update(const EndOfRun* run) override;
0105   void update(const BeginOfEvent* evt) override;
0106   void update(const BeginOfTrack* trk) override;
0107   void update(const G4Step* step) override;
0108   void update(const EndOfTrack* trk) override;
0109   void update(const EndOfEvent* evt) override;
0110 
0111 private:
0112   // Utilities to get detector levels during a step
0113 
0114   int detLevels(const G4VTouchable*) const;
0115   G4String detName(const G4VTouchable*, int, int) const;
0116   void detectorLevel(const G4VTouchable*, int&, int*, G4String*) const;
0117 
0118   //UHB_Analysis* UserNtuples;
0119 
0120   int iev;
0121   int itrk;
0122   G4double entot0, tracklength0;
0123 
0124   double rinCalo, zinCalo;
0125   int lastTrackID;
0126   int verbosity;
0127 
0128   // sumEnerDeposit - all deposited energy on all steps ;  sumStepl - length in steel !!!
0129   G4double sumEnerDeposit, sumStepl, sumStepc;
0130   // numofpart - # particles produced along primary track
0131   int numofpart;
0132   // last point of the primary track
0133   G4ThreeVector lastpo;
0134 
0135   // z:
0136   double z1, z2, z3, z4, zD2, zD3;
0137   int sn0, dn0, pn0, rn0;
0138   int rn00;
0139   double zSiDet, z420;
0140   double zGapLDet, zBoundDet, zSiStep, zSiPlane, zinibeg;
0141   double zBlade, gapBlade;
0142 
0143   Float_t fp420eventarray[1];
0144   TNtuple* fp420eventntuple;
0145   TFile fp420OutputFile;
0146   int whichevent;
0147 
0148   Fp420AnalysisHistManager* theHistManager;  //Histogram Manager of the analysis
0149   std::string fDataLabel;                    // Data type label
0150   std::string fOutputFile;                   // The output file name
0151   std::string fRecreateFile;                 // Recreate the file flag, default="RECREATE"
0152 };
0153 
0154 //================================================================
0155 enum ntfp420_elements { ntfp420_evt };
0156 //================================================================
0157 FP420Test::FP420Test(const edm::ParameterSet& p) {
0158   //constructor
0159   edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("FP420Test");
0160   verbosity = m_Anal.getParameter<int>("Verbosity");
0161   //verbosity    = 1;
0162 
0163   fDataLabel = m_Anal.getParameter<std::string>("FDataLabel");
0164   fOutputFile = m_Anal.getParameter<std::string>("FOutputFile");
0165   fRecreateFile = m_Anal.getParameter<std::string>("FRecreateFile");
0166 
0167   if (verbosity > 0) {
0168     edm::LogVerbatim("FP420Test") << "============================================================================";
0169     edm::LogVerbatim("FP420Test") << "FP420Testconstructor :: Initialized as observer";
0170   }
0171   // Initialization:
0172 
0173   pn0 = 6;
0174   sn0 = 3;
0175   rn00 = 7;
0176 
0177   z420 = 420000.0;  // mm
0178   zD2 = 4000.0;     // mm
0179   zD3 = 8000.0;     // mm
0180   //
0181   zBlade = 5.00;
0182   gapBlade = 1.6;
0183   double gapSupplane = 1.6;
0184   zSiPlane = 2 * zBlade + gapBlade + gapSupplane;
0185 
0186   double zKapton = 0.1;
0187   zSiStep = zSiPlane + zKapton;
0188 
0189   double zBoundDet = 0.020;
0190   double zSiElectr = 0.250;
0191   double zCeramDet = 0.500;
0192   //
0193   zSiDet = 0.250;
0194   zGapLDet = zBlade / 2 - (zSiDet + zSiElectr + zBoundDet + zCeramDet / 2);
0195   //
0196   //  zSiStation = 5*(2*(5.+1.6)+0.1)+2*6.+1.0 =  79.5
0197   double zSiStation = (pn0 - 1) * (2 * (zBlade + gapBlade) + zKapton) + 2 * 6. + 0.0;  // =  78.5
0198   // 11.=e1, 12.=e2 in zzzrectangle.xml
0199   double eee1 = 11.;
0200   double eee2 = 12.;
0201 
0202   zinibeg = (eee1 - eee2) / 2.;
0203 
0204   z1 = zinibeg + (zSiStation + 10.) / 2 + z420;  // z1 - right after 1st Station
0205   z2 = z1 + zD2;                                 //z2 - right after middle Station
0206   z3 = z1 + zD3;                                 //z3 - right after last   Station
0207   z4 = z1 + 2 * zD3;
0208   //==================================
0209 
0210   fp420eventntuple = new TNtuple("NTfp420event", "NTfp420event", "evt");
0211 
0212   whichevent = 0;
0213 
0214   //   fDataLabel      = "defaultData";
0215   //       fOutputFile     = "TheAnlysis.root";
0216   //       fRecreateFile   = "RECREATE";
0217 
0218   theHistManager = new Fp420AnalysisHistManager(fDataLabel);
0219 
0220   if (verbosity > 0) {
0221     edm::LogVerbatim("FP420Test") << "FP420Test constructor :: Initialized Fp420AnalysisHistManager";
0222   }
0223 }
0224 
0225 FP420Test::~FP420Test() {
0226   //  delete UserNtuples;
0227 
0228   TFile fp420OutputFile("newntfp420.root", "RECREATE");
0229   edm::LogVerbatim("FP420Test") << "FP420 output root file has been created";
0230   fp420eventntuple->Write();
0231   edm::LogVerbatim("FP420Test") << "FP420 output root file has been  written";
0232   fp420OutputFile.Close();
0233   edm::LogVerbatim("FP420Test") << "FP420 output root file has been  closed";
0234   delete fp420eventntuple;
0235   edm::LogVerbatim("FP420Test") << "FP420 output root file pointer has been  deleted";
0236 
0237   //------->while end
0238 
0239   // Write histograms to file
0240   theHistManager->writeToFile(fOutputFile, fRecreateFile);
0241   if (verbosity > 0) {
0242     edm::LogVerbatim("FP420Test") << std::endl << "FP420Test Destructor  -------->  End of FP420Test : ";
0243   }
0244 }
0245 
0246 //================================================================
0247 // Histoes:
0248 //-----------------------------------------------------------------------------
0249 
0250 Fp420AnalysisHistManager::Fp420AnalysisHistManager(const TString& managername) {
0251   // The Constructor
0252 
0253   fTypeTitle = managername;
0254   fHistArray = new TObjArray();       // Array to store histos
0255   fHistNamesArray = new TObjArray();  // Array to store histos's names
0256 
0257   TH1::AddDirectory(kFALSE);
0258   bookHistos();
0259 
0260   fHistArray->Compress();  // Removes empty space
0261   fHistNamesArray->Compress();
0262 
0263   //      storeWeights();                    // Store the weights
0264 }
0265 //-----------------------------------------------------------------------------
0266 
0267 Fp420AnalysisHistManager::~Fp420AnalysisHistManager() {
0268   // The Destructor
0269 
0270   if (fHistArray) {
0271     fHistArray->Delete();
0272     delete fHistArray;
0273   }
0274 
0275   if (fHistNamesArray) {
0276     fHistNamesArray->Delete();
0277     delete fHistNamesArray;
0278   }
0279 }
0280 //-----------------------------------------------------------------------------
0281 
0282 void Fp420AnalysisHistManager::bookHistos() {
0283   // at Start: (mm)
0284   histInit("PrimaryEta", "Primary Eta", 100, 9., 12.);
0285   histInit("PrimaryPhigrad", "Primary Phigrad", 100, 0., 360.);
0286   histInit("PrimaryTh", "Primary Th", 100, 0., 180.);
0287   histInit("PrimaryLastpoZ", "Primary Lastpo Z", 100, -200., 430000.);
0288   histInit("PrimaryLastpoX", "Primary Lastpo X Z<z4", 100, -30., 30.);
0289   histInit("PrimaryLastpoY", "Primary Lastpo Y Z<z4", 100, -30., 30.);
0290   histInit("XLastpoNumofpart", "Primary Lastpo X n>10", 100, -30., 30.);
0291   histInit("YLastpoNumofpart", "Primary Lastpo Y n>10", 100, -30., 30.);
0292   histInit("VtxX", "Vtx X", 100, -50., 50.);
0293   histInit("VtxY", "Vtx Y", 100, -50., 50.);
0294   histInit("VtxZ", "Vtx Z", 100, -200., 430000.);
0295   // Book the histograms and add them to the array
0296   histInit("SumEDep", "This is sum Energy deposited", 100, -1, 199.);
0297   histInit("TrackL", "This is TrackL", 100, 0., 12000.);
0298   histInit("zHits", "z Hits all events", 100, 400000., 430000.);
0299   histInit("zHitsnoMI", "z Hits no MI", 100, 400000., 430000.);
0300   histInit("zHitsTrLoLe", "z Hits TrLength bigger 8300", 100, 400000., 430000.);
0301   histInit("NumberOfHits", "NumberOfHits", 100, 0., 300.);
0302 }
0303 
0304 //-----------------------------------------------------------------------------
0305 
0306 void Fp420AnalysisHistManager::writeToFile(const TString& fOutputFile, const TString& fRecreateFile) {
0307   //Write to file = fOutputFile
0308 
0309   edm::LogVerbatim("FP420Test") << "================================================================";
0310   edm::LogVerbatim("FP420Test") << " Write this Analysis to File " << fOutputFile;
0311   edm::LogVerbatim("FP420Test") << "================================================================";
0312 
0313   TFile* file = new TFile(fOutputFile, fRecreateFile);
0314 
0315   fHistArray->Write();
0316   file->Close();
0317 }
0318 //-----------------------------------------------------------------------------
0319 
0320 void Fp420AnalysisHistManager::histInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup) {
0321   // Add histograms and histograms names to the array
0322 
0323   char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
0324   strcpy(newtitle, title);
0325   strcat(newtitle, " (");
0326   strcat(newtitle, fTypeTitle);
0327   strcat(newtitle, ") ");
0328   fHistArray->AddLast((new TH1F(name, newtitle, nbinsx, xlow, xup)));
0329   fHistNamesArray->AddLast(new TObjString(name));
0330 }
0331 //-----------------------------------------------------------------------------
0332 
0333 void Fp420AnalysisHistManager::histInit(
0334     const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup) {
0335   // Add histograms and histograms names to the array
0336 
0337   char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
0338   strcpy(newtitle, title);
0339   strcat(newtitle, " (");
0340   strcat(newtitle, fTypeTitle);
0341   strcat(newtitle, ") ");
0342   fHistArray->AddLast((new TH2F(name, newtitle, nbinsx, xlow, xup, nbinsy, ylow, yup)));
0343   fHistNamesArray->AddLast(new TObjString(name));
0344 }
0345 //-----------------------------------------------------------------------------
0346 
0347 TH1F* Fp420AnalysisHistManager::getHisto(Int_t Number) {
0348   // Get a histogram from the array with index = Number
0349 
0350   if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
0351     return (TH1F*)(fHistArray->At(Number));
0352 
0353   } else {
0354     edm::LogVerbatim("FP420Test") << "!!!!!!!!!!!!!!!!!!getHisto!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0355     edm::LogVerbatim("FP420Test") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
0356     edm::LogVerbatim("FP420Test") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0357 
0358     return (TH1F*)(fHistArray->At(0));
0359   }
0360 }
0361 //-----------------------------------------------------------------------------
0362 
0363 TH2F* Fp420AnalysisHistManager::getHisto2(Int_t Number) {
0364   // Get a histogram from the array with index = Number
0365 
0366   if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
0367     return (TH2F*)(fHistArray->At(Number));
0368 
0369   } else {
0370     edm::LogVerbatim("FP420Test") << "!!!!!!!!!!!!!!!!getHisto2!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0371     edm::LogVerbatim("FP420Test") << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number;
0372     edm::LogVerbatim("FP420Test") << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
0373 
0374     return (TH2F*)(fHistArray->At(0));
0375   }
0376 }
0377 //-----------------------------------------------------------------------------
0378 
0379 TH1F* Fp420AnalysisHistManager::getHisto(const TObjString& histname) {
0380   // Get a histogram from the array with name = histname
0381 
0382   Int_t index = fHistNamesArray->IndexOf(&histname);
0383   return getHisto(index);
0384 }
0385 //-----------------------------------------------------------------------------
0386 
0387 TH2F* Fp420AnalysisHistManager::getHisto2(const TObjString& histname) {
0388   // Get a histogram from the array with name = histname
0389 
0390   Int_t index = fHistNamesArray->IndexOf(&histname);
0391   return getHisto2(index);
0392 }
0393 //-----------------------------------------------------------------------------
0394 
0395 void Fp420AnalysisHistManager::storeWeights() {
0396   // Add structure to each histogram to store the weights
0397 
0398   for (int i = 0; i < fHistArray->GetEntries(); i++) {
0399     ((TH1F*)(fHistArray->At(i)))->Sumw2();
0400   }
0401 }
0402 
0403 // Histoes end :
0404 
0405 //================================================================
0406 
0407 // using several observers
0408 
0409 //==================================================================== per JOB
0410 void FP420Test::update(const BeginOfJob* job) {
0411   //job
0412   edm::LogVerbatim("FP420Test") << "FP420Test:beggining of job";
0413   ;
0414 }
0415 
0416 //==================================================================== per RUN
0417 void FP420Test::update(const BeginOfRun* run) {
0418   //run
0419 
0420   edm::LogVerbatim("FP420Test") << std::endl << "FP420Test:: Begining of Run";
0421 }
0422 
0423 void FP420Test::update(const EndOfRun* run) { ; }
0424 
0425 //=================================================================== per EVENT
0426 void FP420Test::update(const BeginOfEvent* evt) {
0427   iev = (*evt)()->GetEventID();
0428   if (verbosity > 0) {
0429     edm::LogVerbatim("FP420Test") << "FP420Test:update Event number = " << iev;
0430   }
0431   whichevent++;
0432 }
0433 
0434 //=================================================================== per Track
0435 void FP420Test::update(const BeginOfTrack* trk) {
0436   itrk = (*trk)()->GetTrackID();
0437   if (verbosity > 1) {
0438     edm::LogVerbatim("FP420Test") << "FP420Test:update BeginOfTrack number = " << itrk;
0439   }
0440   if (itrk == 1) {
0441     sumEnerDeposit = 0.;
0442     numofpart = 0;
0443     sumStepl = 0.;
0444     sumStepc = 0.;
0445     tracklength0 = 0.;
0446   }
0447 }
0448 
0449 //=================================================================== per EndOfTrack
0450 void FP420Test::update(const EndOfTrack* trk) {
0451   itrk = (*trk)()->GetTrackID();
0452   if (verbosity > 1) {
0453     edm::LogVerbatim("FP420Test") << "FP420Test:update EndOfTrack number = " << itrk;
0454   }
0455   if (itrk == 1) {
0456     G4double tracklength = (*trk)()->GetTrackLength();  // Accumulated track length
0457 
0458     theHistManager->getHisto("SumEDep")->Fill(sumEnerDeposit);
0459     theHistManager->getHisto("TrackL")->Fill(tracklength);
0460 
0461     // direction !!!
0462     G4ThreeVector vert_mom = (*trk)()->GetVertexMomentumDirection();
0463     G4ThreeVector vert_pos = (*trk)()->GetVertexPosition();  // vertex ,where this track was created
0464 
0465     //    float eta = 0.5 * log( (1.+vert_mom.z()) / (1.-vert_mom.z()) );
0466     //    float phi = atan2(vert_mom.y(),vert_mom.x());
0467     //    if (phi < 0.) phi += twopi;
0468     //    float phigrad = phi*180./pi;
0469 
0470     //      float XV = vert_pos.x(); // mm
0471     //      float YV = vert_pos.y(); // mm
0472     //UserNtuples->fillg543(phigrad,1.);
0473     //UserNtuples->fillp203(phigrad,sumStepl,1.);
0474     //UserNtuples->fillg544(XV,1.);
0475     //UserNtuples->fillp201(XV,sumStepl,1.);
0476     //UserNtuples->fillg545(YV,1.);
0477     //UserNtuples->fillp202(YV,sumStepl,1.);
0478 
0479     //UserNtuples->fillg524(eta,1.);
0480 
0481     //UserNtuples->fillg534(sumStepl,1.);
0482     //UserNtuples->fillg535(eta,sumStepl);
0483     //UserNtuples->fillp207(eta,sumStepl,1.);
0484     //UserNtuples->filld217(eta,sumStepl,1.);
0485     //UserNtuples->filld220(phigrad,sumStepl,1.);
0486     //UserNtuples->filld221(XV,sumStepl,1.);
0487     //UserNtuples->filld222(YV,sumStepl,1.);
0488     //UserNtuples->fillg537(sumStepc,1.);
0489     //UserNtuples->fillg84(sumStepl,1.);
0490 
0491     // MI = (multiple interactions):
0492     if (tracklength < z4) {
0493       //        //UserNtuples->fillp208(eta,tracklength,1.);
0494       //UserNtuples->filld218(eta,tracklength,1.);
0495       //UserNtuples->fillg538(sumStepc,1.);
0496       //UserNtuples->fillg85(sumStepl,1.);
0497     }
0498 
0499     // last step information
0500     const G4Step* aStep = (*trk)()->GetStep();
0501     //   G4int csn = (*trk)()->GetCurrentStepNumber();
0502     //   G4double sl = (*trk)()->GetStepLength();
0503     // preStep
0504     G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0505     lastpo = preStepPoint->GetPosition();
0506 
0507     // Analysis:
0508     if (lastpo.z() < z1 && lastpo.perp() < 100.) {
0509       //UserNtuples->fillg525(eta,1.);
0510       //UserNtuples->fillg546(XV,1.);
0511       //UserNtuples->fillg551(YV,1.);
0512       //UserNtuples->fillg556(phigrad,1.);
0513     }
0514     if ((lastpo.z() > z1 && lastpo.z() < z2) && lastpo.perp() < 100.) {
0515       //UserNtuples->fillg526(eta,1.);
0516       //UserNtuples->fillg547(XV,1.);
0517       //UserNtuples->fillg552(YV,1.);
0518       //UserNtuples->fillg557(phigrad,1.);
0519     }
0520     if (lastpo.z() < z2 && lastpo.perp() < 100.) {
0521       //UserNtuples->fillg527(eta,1.);
0522       //UserNtuples->fillg548(XV,1.);
0523       //UserNtuples->fillg553(YV,1.);
0524       //UserNtuples->fillg558(phigrad,1.);
0525       //UserNtuples->fillg521(lastpo.x(),1.);
0526       //UserNtuples->fillg522(lastpo.y(),1.);
0527       //UserNtuples->fillg523(lastpo.z(),1.);
0528     }
0529     if (lastpo.z() < z3 && lastpo.perp() < 100.) {
0530       //UserNtuples->fillg528(eta,1.);
0531       //UserNtuples->fillg549(XV,1.);
0532       //UserNtuples->fillg554(YV,1.);
0533       //UserNtuples->fillg559(phigrad,1.);
0534     }
0535     if (lastpo.z() < z4 && lastpo.perp() < 100.) {
0536       //UserNtuples->fillg529(eta,1.);
0537       //UserNtuples->fillg550(XV,1.);
0538       //UserNtuples->fillg555(YV,1.);
0539       //UserNtuples->fillg560(phigrad,1.);
0540       //UserNtuples->fillg531(lastpo.x(),1.);
0541       //UserNtuples->fillg532(lastpo.y(),1.);
0542       //UserNtuples->fillg533(lastpo.z(),1.);
0543     }
0544   }
0545 }
0546 
0547 // ====================================================
0548 
0549 //=================================================================== each STEP
0550 void FP420Test::update(const G4Step* aStep) {
0551   // ==========================================================================
0552 
0553   if (verbosity > 2) {
0554     G4int stepnumber = aStep->GetTrack()->GetCurrentStepNumber();
0555     edm::LogVerbatim("FP420Test") << "FP420Test:update Step number = " << stepnumber;
0556   }
0557   // track on aStep:                                                                                         !
0558   G4Track* theTrack = aStep->GetTrack();
0559   TrackInformation* trkInfo = dynamic_cast<TrackInformation*>(theTrack->GetUserInformation());
0560   if (trkInfo == nullptr) {
0561     edm::LogVerbatim("FP420Test") << "FP420Test on aStep: No trk info !!!! abort ";
0562   }
0563   G4int id = theTrack->GetTrackID();
0564   G4String particleType = theTrack->GetDefinition()->GetParticleName();  //   !!!
0565   G4int parentID = theTrack->GetParentID();                              //   !!!
0566   G4TrackStatus trackstatus = theTrack->GetTrackStatus();                //   !!!
0567   G4double tracklength = theTrack->GetTrackLength();                     // Accumulated track length
0568   G4ThreeVector trackmom = theTrack->GetMomentum();
0569   G4double entot = theTrack->GetTotalEnergy();  //   !!! deposited on step
0570   G4int curstepnumber = theTrack->GetCurrentStepNumber();
0571   // const G4ThreeVector&   vert_pos       = theTrack->GetVertexPosition(); // vertex ,where this track was created
0572   // const G4ThreeVector&   vert_mom       = theTrack->GetVertexMomentumDirection();
0573 
0574   //   double costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+vert_mom.y()*vert_mom.y()+vert_mom.z()*vert_mom.z());
0575   //   double theta = acos(min(max(costheta,double(-1.)),double(1.)));
0576   //  float eta = -log(tan(theta/2));
0577   //   double phi = -1000.;
0578   //   if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x());
0579   //   if (phi < 0.) phi += twopi;
0580   //   double phigrad = phi*360./twopi;
0581 
0582   //G4double       trackvel       = theTrack->GetVelocity();
0583 
0584   //edm::LogVerbatim("FP420Test") << " trackstatus= " << trackstatus << " entot= " << entot ;
0585 
0586   // step points:                                                                                         !
0587   G4double stepl = aStep->GetStepLength();
0588   G4double EnerDeposit = aStep->GetTotalEnergyDeposit();
0589   // pointers:                                                                                         !
0590   //G4VPhysicalVolume*  physvol       = theTrack->GetVolume();
0591   //G4VPhysicalVolume*  nextphysvol   = theTrack->GetNextVolume();
0592   //G4Material*       materialtr     = theTrack->GetMaterial();
0593   //G4Material*       nextmaterialtr = theTrack->GetNextMaterial();
0594 
0595   // preStep
0596   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0597   const G4ThreeVector& preposition = preStepPoint->GetPosition();
0598   G4ThreeVector prelocalpoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(preposition);
0599   G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
0600   const G4String& prename = currentPV->GetName();
0601 
0602   const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
0603   int pre_levels = detLevels(pre_touch);
0604 
0605   G4String name1[20];
0606   int copyno1[20];
0607   for (int i = 0; i < 20; ++i) {
0608     name1[i] = "";
0609     copyno1[i] = 0;
0610   }
0611   if (pre_levels > 0) {
0612     detectorLevel(pre_touch, pre_levels, copyno1, name1);
0613   }
0614 
0615   // gen_track:
0616   //  id=1    parentID=1   trackstatus=0,2   tracklength(accumulated)  curstepnumber   entot  vert_pos
0617   if (id == 1) {
0618     // on 1-st step:
0619     if (curstepnumber == 1) {
0620       entot0 = entot;
0621       //UserNtuples->fillg519(entot0,1.);
0622     }
0623 
0624     // on every step:
0625 
0626     // for Copper:
0627     if (prename == "SBST") {
0628       sumStepc += stepl;
0629       // =========
0630     }
0631     // for ststeel:
0632     //   if(prename == "SBSTs") {
0633     if (prename == "SBSTs") {
0634       sumStepl += stepl;
0635       // =========
0636     }
0637     // =========
0638     // =========
0639 
0640     // exclude last track point if it is in SD (MI was started their)
0641     if (trackstatus != 2) {
0642       // for SD:   Si Det.:   SISTATION:SIPLANE:(SIDETL+BOUNDDET        +SIDETR + CERAMDET)
0643       if (prename == "SIDETL" || prename == "SIDETR") {
0644         if (prename == "SIDETL") {
0645           //UserNtuples->fillg569(EnerDeposit,1.);
0646         }
0647         if (prename == "SIDETR") {
0648           //UserNtuples->fillg570(EnerDeposit,1.);
0649         }
0650 
0651         //     double numbcont = 10.*(copyno1[2]-1)+copyno1[3];
0652 
0653         // =========
0654         //     double   xx    = preposition.x();
0655         //     double   yy    = preposition.y();
0656         //     double   zz    = preposition.z();
0657         // =========
0658         //UserNtuples->fillg580(theta,1.);
0659         //UserNtuples->fillg07(phigrad,1.);
0660         //       double xPrimAtZhit = vert_pos.x() + (zz-vert_pos.z())*tan(theta)*cos(phi);
0661         //       double yPrimAtZhit = vert_pos.y() + (zz-vert_pos.z())*tan(theta)*sin(phi);
0662         // =========
0663         //     double  dx = xPrimAtZhit - xx;
0664         //     double  dy = yPrimAtZhit - yy;
0665         // =========
0666 
0667         //     //UserNtuples->fillp212(numbcont,dx,1.);
0668         //     //UserNtuples->fillp213(numbcont,dy,1.);
0669         // =========
0670 
0671         // last step of current SD Volume:
0672         G4String posname = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
0673         if ((prename == "SIDETL" && posname != "SIDETL") || (prename == "SIDETR" && posname != "SIDETR")) {
0674           if (name1[2] == "SISTATION") {
0675             //UserNtuples->fillg539(copyno1[2],1.);
0676           }
0677           if (name1[3] == "SIPLANE") {
0678             //UserNtuples->fillg540(copyno1[3],1.);
0679           }
0680 
0681           if (prename == "SIDETL") {
0682             //UserNtuples->fillg541(EnerDeposit,1.);
0683             //UserNtuples->fillg561(numbcont,1.);
0684             if (copyno1[2] < 2) {
0685               //UserNtuples->fillg571(dx,1.);
0686             } else if (copyno1[2] < 3) {
0687               //UserNtuples->fillg563(dx,1.);
0688               if (copyno1[3] < 2) {
0689               } else if (copyno1[3] < 3) {
0690                 //UserNtuples->fillg572(dx,1.);
0691               } else if (copyno1[3] < 4) {
0692                 //UserNtuples->fillg573(dx,1.);
0693               } else if (copyno1[3] < 5) {
0694                 //UserNtuples->fillg574(dx,1.);
0695               } else if (copyno1[3] < 6) {
0696                 //UserNtuples->fillg575(dx,1.);
0697               } else if (copyno1[3] < 7) {
0698                 //UserNtuples->fillg576(dx,1.);
0699               } else if (copyno1[3] < 8) {
0700                 //UserNtuples->fillg577(dx,1.);
0701               } else if (copyno1[3] < 9) {
0702                 //UserNtuples->fillg578(dx,1.);
0703               } else if (copyno1[3] < 10) {
0704                 //UserNtuples->fillg579(dx,1.);
0705               }
0706             } else if (copyno1[2] < 4) {
0707               //UserNtuples->fillg565(dx,1.);
0708             } else if (copyno1[2] < 5) {
0709               //UserNtuples->fillg567(dx,1.);
0710             }
0711           }
0712           if (prename == "SIDETR") {
0713             //UserNtuples->fillg542(EnerDeposit,1.);
0714             //UserNtuples->fillg562(numbcont,1.);
0715             if (copyno1[2] < 2) {
0716               //UserNtuples->fillg581(dy,1.);
0717             } else if (copyno1[2] < 3) {
0718               //UserNtuples->fillg564(dy,1.);
0719               if (copyno1[3] < 2) {
0720               } else if (copyno1[3] < 3) {
0721                 //UserNtuples->fillg582(dy,1.);
0722               } else if (copyno1[3] < 4) {
0723                 //UserNtuples->fillg583(dy,1.);
0724               } else if (copyno1[3] < 5) {
0725                 //UserNtuples->fillg584(dy,1.);
0726               } else if (copyno1[3] < 6) {
0727                 //UserNtuples->fillg585(dy,1.);
0728               } else if (copyno1[3] < 7) {
0729                 //UserNtuples->fillg586(dy,1.);
0730               } else if (copyno1[3] < 8) {
0731                 //UserNtuples->fillg587(dy,1.);
0732               } else if (copyno1[3] < 9) {
0733                 //UserNtuples->fillg588(dy,1.);
0734               } else if (copyno1[3] < 10) {
0735                 //UserNtuples->fillg589(dy,1.);
0736               }
0737             } else if (copyno1[2] < 4) {
0738               //UserNtuples->fillg566(dy,1.);
0739             } else if (copyno1[2] < 5) {
0740               //UserNtuples->fillg568(dy,1.);
0741             }
0742           }
0743         }
0744       }
0745       // end of prenames SIDETL // SIDETR
0746     }
0747     // end of trackstatus != 2
0748 
0749     // deposition of energy on steps along primary track
0750     //UserNtuples->fillg500(EnerDeposit,1.);
0751     // collect sum deposited energy on all steps along primary track
0752     sumEnerDeposit += EnerDeposit;
0753     // position of step for primary track:
0754     //UserNtuples->fillg501(preposition.x(),1.);
0755     //UserNtuples->fillg502(preposition.y(),1.);
0756     //UserNtuples->fillg503(preposition.z(),1.);
0757     //UserNtuples->fillg504(preposition.x(),EnerDeposit);
0758     //UserNtuples->fillg505(preposition.y(),EnerDeposit);
0759     //UserNtuples->fillg506(preposition.z(),EnerDeposit);
0760     // 2D step coordinates weighted by energy deposited on step
0761     //      //UserNtuples->fillp201(preposition.x(),preposition.y(),EnerDeposit);
0762     //      //UserNtuples->fillp202(preposition.x(),preposition.z(),EnerDeposit);
0763     //      //UserNtuples->fillp203(preposition.y(),preposition.z(),EnerDeposit);
0764     //UserNtuples->filld204(preposition.x(),preposition.y(),EnerDeposit);
0765     //UserNtuples->filld205(preposition.x(),preposition.z(),EnerDeposit);
0766     //UserNtuples->filld206(preposition.y(),preposition.z(),EnerDeposit);
0767     //UserNtuples->filld223(preposition.y(),preposition.z(),EnerDeposit);
0768     // last step of primary track
0769     if (trackstatus == 2) {
0770       // primary track length
0771       //      //UserNtuples->fillg508(tracklength,1.);
0772       tracklength0 = tracklength;
0773       // how many steps primary track consist
0774       //UserNtuples->fillg509(curstepnumber,1.);
0775       // tot. energy of primary track at the end of trajectory(before disappeare)
0776       //UserNtuples->fillg510((entot0-entot),1.);
0777       //UserNtuples->fillg520((entot0-entot),1.);
0778     }
0779   }
0780   // end of primary track !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0781 
0782   if (parentID == 1 && curstepnumber == 1) {
0783     // particles deposit their energy along primary track
0784     numofpart += 1;
0785     // energy of radiated particle
0786     //UserNtuples->fillg511(entot,1.);
0787     // coordinates  of radiated particle
0788     //UserNtuples->fillg512(vert_pos.x(),1.);
0789     //UserNtuples->fillg513(vert_pos.y(),1.);
0790     //UserNtuples->fillg514(vert_pos.z(),1.);
0791     //UserNtuples->fillg515(vert_pos.x(),entot);
0792     //UserNtuples->fillg516(vert_pos.y(),entot);
0793     //UserNtuples->fillg517(vert_pos.z(),entot);
0794 
0795     //UserNtuples->filld214(vert_pos.x(),vert_pos.y(),1.);
0796     //UserNtuples->filld215(vert_pos.x(),vert_pos.z(),1.);
0797     //UserNtuples->filld216(vert_pos.y(),vert_pos.z(),1.);
0798     //UserNtuples->filld219(vert_pos.y(),vert_pos.z(),1.);
0799     //UserNtuples->filld224(vert_pos.y(),vert_pos.z(),1.);
0800     if (prename == "SBST") {
0801       //UserNtuples->filld225(vert_pos.y(),vert_pos.z(),1.);
0802     }
0803     if (prename == "SBSTs") {
0804       //UserNtuples->filld226(vert_pos.y(),vert_pos.z(),1.);
0805     }
0806   }
0807 
0808   // ==========================================================================
0809 }
0810 // ==========================================================================
0811 // ==========================================================================
0812 int FP420Test::detLevels(const G4VTouchable* touch) const {
0813   //Return number of levels
0814   if (touch)
0815     return ((touch->GetHistoryDepth()) + 1);
0816   else
0817     return 0;
0818 }
0819 // ==========================================================================
0820 
0821 G4String FP420Test::detName(const G4VTouchable* touch, int level, int currentlevel) const {
0822   //Go down to current level
0823   if (level > 0 && level >= currentlevel) {
0824     int ii = level - currentlevel;
0825     return touch->GetVolume(ii)->GetName();
0826   } else {
0827     return "NotFound";
0828   }
0829 }
0830 
0831 void FP420Test::detectorLevel(const G4VTouchable* touch, int& level, int* copyno, G4String* name) const {
0832   //Get name and copy numbers
0833   if (level > 0) {
0834     for (int ii = 0; ii < level; ii++) {
0835       int i = level - ii - 1;
0836       G4VPhysicalVolume* pv = touch->GetVolume(i);
0837       if (pv != nullptr)
0838         name[ii] = pv->GetName();
0839       else
0840         name[ii] = "Unknown";
0841       copyno[ii] = touch->GetReplicaNumber(i);
0842     }
0843   }
0844 }
0845 // ==========================================================================
0846 
0847 //===================================================================   End Of Event
0848 void FP420Test::update(const EndOfEvent* evt) {
0849   // ==========================================================================
0850 
0851   if (verbosity > 1) {
0852     iev = (*evt)()->GetEventID();
0853     edm::LogVerbatim("FP420Test") << "FP420Test:update EndOfEvent = " << iev;
0854   }
0855   // Fill-in ntuple
0856   fp420eventarray[ntfp420_evt] = (float)whichevent;
0857 
0858   //
0859   int trackID = 0;
0860   G4PrimaryParticle* thePrim = nullptr;
0861 
0862   // prim.vertex:
0863   G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0864   if (nvertex != 1)
0865     edm::LogVerbatim("FP420Test") << "FP420Test: My warning: NumberOfPrimaryVertex != 1  -->  = " << nvertex;
0866 
0867   for (int i = 0; i < nvertex; i++) {
0868     G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
0869     if (avertex == nullptr) {
0870       edm::LogVerbatim("FP420Test") << "FP420Test  End Of Event ERR: pointer to vertex = 0";
0871       continue;
0872     }
0873     G4int npart = avertex->GetNumberOfParticle();
0874     if (npart != 1)
0875       edm::LogVerbatim("FP420Test") << "FP420Test: My warning: NumberOfPrimaryPart != 1  -->  = " << npart;
0876     if (npart == 0)
0877       edm::LogVerbatim("FP420Test") << "FP420Test End Of Event ERR: no NumberOfParticle";
0878 
0879     // find just primary track:                                                             track pointer: thePrim
0880     if (thePrim == nullptr)
0881       thePrim = avertex->GetPrimary(trackID);
0882 
0883     if (thePrim != nullptr) {
0884       // primary vertex:
0885       G4double vx = 0., vy = 0., vz = 0.;
0886       vx = avertex->GetX0();
0887       vy = avertex->GetY0();
0888       vz = avertex->GetZ0();
0889       //UserNtuples->fillh01(vx);
0890       //UserNtuples->fillh02(vy);
0891       //UserNtuples->fillh03(vz);
0892       theHistManager->getHisto("VtxX")->Fill(vx);
0893       theHistManager->getHisto("VtxY")->Fill(vy);
0894       theHistManager->getHisto("VtxZ")->Fill(vz);
0895     }
0896   }
0897   // prim.vertex loop end
0898 
0899   //=========================== thePrim != 0 ================================================================================
0900   if (thePrim != nullptr) {
0901     //      inline G4ParticleDefinition * GetG4code() const
0902     //      inline G4PrimaryParticle * GetNext() const
0903     //      inline G4PrimaryParticle * GetDaughter() const
0904 
0905     // prim.vertex
0906     //int ivert = 0 ;
0907     //G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(ivert);
0908     //G4double vx=0.,vy=0.,vz=0.;
0909     //vx = avertex->GetX0();
0910     //vy = avertex->GetY0();
0911     //vz = avertex->GetZ0();
0912 
0913     //
0914     // number of secondary particles deposited their energy along primary track
0915     //UserNtuples->fillg518(numofpart,1.);
0916     if (lastpo.z() < z4 && lastpo.perp() < 100.) {
0917       //UserNtuples->fillg536(numofpart,1.);
0918     }
0919     //
0920 
0921     // direction !!!
0922     G4ThreeVector mom = thePrim->GetMomentum();
0923 
0924     double phi = atan2(mom.y(), mom.x());
0925     if (phi < 0.)
0926       phi += twopi;
0927     double phigrad = phi * 180. / pi;
0928 
0929     double th = mom.theta();
0930     double eta = -log(tan(th / 2));
0931     // works OK:
0932     //      double  costheta =mom.z()/sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z());
0933     //      double theta = acos(min(max(costheta,double(-1.)),double(1.)));
0934 
0935     //UserNtuples->fillh04(eta);
0936     //UserNtuples->fillh05(phigrad);
0937     //UserNtuples->fillh06(th);
0938 
0939     //UserNtuples->fillg13(lastpo.x(),1.);
0940     //UserNtuples->fillg14(lastpo.y(),1.);
0941     //UserNtuples->fillg15(lastpo.z(),1.);
0942 
0943     theHistManager->getHisto("PrimaryEta")->Fill(eta);
0944     theHistManager->getHisto("PrimaryPhigrad")->Fill(phigrad);
0945     theHistManager->getHisto("PrimaryTh")->Fill(th * 180. / pi);
0946 
0947     theHistManager->getHisto("PrimaryLastpoZ")->Fill(lastpo.z());
0948     if (lastpo.z() < z4) {
0949       theHistManager->getHisto("PrimaryLastpoX")->Fill(lastpo.x());
0950       theHistManager->getHisto("PrimaryLastpoY")->Fill(lastpo.y());
0951     }
0952     if (numofpart > 4) {
0953       theHistManager->getHisto("XLastpoNumofpart")->Fill(lastpo.x());
0954       theHistManager->getHisto("YLastpoNumofpart")->Fill(lastpo.y());
0955     }
0956 
0957     // ==========================================================================
0958 
0959     // hit map for FP420
0960     // ==================================
0961 
0962     std::map<int, float, std::less<int> > themap;
0963     std::map<int, float, std::less<int> > themap1;
0964 
0965     std::map<int, float, std::less<int> > themapxy;
0966     std::map<int, float, std::less<int> > themapz;
0967     // access to the G4 hit collections:  -----> this work OK:
0968 
0969     //  edm::LogInfo("FP420Test") << "1";
0970     G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0971     //  edm::LogInfo("FP420Test") << "2";
0972     if (verbosity > 0) {
0973       edm::LogVerbatim("FP420Test") << "FP420Test:  accessed all HC";
0974       ;
0975     }
0976     int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("FP420SI");
0977     // edm::LogInfo("FP420Test") << "3";
0978     // edm::LogVerbatim("FP420Test") << " CAFIid = " << CAFIid;;
0979 
0980     FP420G4HitCollection* theCAFI = (FP420G4HitCollection*)allHC->GetHC(CAFIid);
0981     //  CaloG4HitCollection* theCAFI = (CaloG4HitCollection*) allHC->GetHC(CAFIid);
0982     if (verbosity > 0) {
0983       //edm::LogVerbatim("FP420Test") << "FP420Test: theCAFI = " << theCAFI;
0984       edm::LogVerbatim("FP420Test") << "FP420Test: theCAFI->entries = " << theCAFI->entries();
0985     }
0986     // edm::LogInfo("FP420Test") << "theCAFI->entries="<< theCAFI->entries();
0987     theHistManager->getHisto("NumberOfHits")->Fill(theCAFI->entries());
0988 
0989     // access to the G4 hit collections ----> variant 2: give 0 hits
0990     //  FP420G4HitCollection *   theCAFI;
0991     //  theCAFI = new FP420G4HitCollection();
0992     // ==========================================================================
0993     //   Silicon Hit collection start
0994     //0) if particle goes into flat beam pipe below detector:
0995     int varia;  // = 0 -all; =1 - MI; =2 - noMI
0996     //                      Select MI or noMI over all 3 stations
0997     // 1)MI:
0998     //     if particle goes through window into detector:
0999     // lastpoint of track in lateral dir. outside the detector and inside in z
1000     // lastpoint of track in lateral dir. outside the detector and inside in z
1001     // for all except zzzmarta.xml
1002     //    if(  lastpo.z()<z4  ||  abs(lastpo.x())> 5. || lastpo.y()< 10.2 || lastpo.y()>30.2   ) {
1003     // for zzzmarta.xml
1004     //   if(  lastpo.z()<z4  ||  abs(lastpo.x())> 10. || lastpo.y()< 3.2 || lastpo.y()>43.2   ) {
1005     if (lastpo.z() < z4) {
1006       //  if(  lastpo.z()<z4 && lastpo.perp()< 100. ) {
1007       //  if(lastpo.z()<z4  || lastpo.perp()> 45.) {
1008       //UserNtuples->fillg66(theCAFI->entries(),1.);
1009       varia = 1;
1010     } else {
1011       // 2)   no MI start in detector:
1012       //UserNtuples->fillg31(numofpart,1.);
1013       //UserNtuples->fillg67(theCAFI->entries(),1.);
1014       varia = 2;
1015     }  // no MI end:
1016     int nhits = theCAFI->entries();
1017     for (int j = 0; j < nhits; j++) {
1018       FP420G4Hit* aHit = (*theCAFI)[j];
1019       G4ThreeVector hitPoint = aHit->getEntry();
1020       double zz = hitPoint.z();
1021       theHistManager->getHisto("zHits")->Fill(zz);
1022       if (tracklength0 > 8300.)
1023         theHistManager->getHisto("zHitsTrLoLe")->Fill(zz);
1024     }
1025     // varia = 0;
1026     //     if( varia == 0) {
1027     if (varia == 2) {
1028       // .............
1029       // number of hits < 50
1030       //    if(theCAFI->entries() <50) {
1031       //    if(theCAFI->entries() > 0) {
1032       //    if(theCAFI->entries() > -1) {
1033       // .............
1034       //int nhit11 = 0, nhit12 = 0, nhit13 = 0;
1035       double totallosenergy = 0.;
1036       for (int j = 0; j < nhits; j++) {
1037         FP420G4Hit* aHit = (*theCAFI)[j];
1038 
1039         G4ThreeVector hitEntryLocalPoint = aHit->getEntryLocalP();
1040         G4ThreeVector hitExitLocalPoint = aHit->getExitLocalP();
1041         G4ThreeVector hitPoint = aHit->getEntry();
1042         //    double  elmenergy =  aHit->getEM();
1043         //    double  hadrenergy =  aHit->getHadr();
1044         //    double incidentEnergyHit  = aHit->getIncidentEnergy();
1045         int trackIDhit = aHit->getTrackID();
1046         unsigned int unitID = aHit->getUnitID();
1047         //    double   timeslice = aHit->getTimeSlice();
1048         //    int     timesliceID = aHit->getTimeSliceID();
1049         //    double  depenergy = aHit->getEnergyDeposit();
1050         //    float   pabs = aHit->getPabs();
1051         //    float   tof = aHit->getTof();
1052         double losenergy = aHit->getEnergyLoss();
1053         //    int   particletype = aHit->getParticleType();
1054         //    float thetaEntry = aHit->getThetaAtEntry();
1055         //    float phiEntry = aHit->getPhiAtEntry();
1056         //    float xEntry = aHit->getX();
1057         //    float yEntry = aHit->getY();
1058         //    float zEntry = aHit->getZ();
1059         //    int  parentID = aHit->getParentId();
1060         //    float vxget = aHit->getVx();
1061         //    float vyget = aHit->getVy();
1062         //    float vzget = aHit->getVz();
1063 
1064         //    double th_hit    = hitPoint.theta();
1065         //    double eta_hit = -log(tan(th_hit/2));
1066         //    double phi_hit   = hitPoint.phi();
1067         //    if (phi_hit < 0.) phi_hit += twopi;
1068         //    double phigrad_hit = phi_hit*180./pi;
1069         //UserNtuples->fillg60(eta_hit,losenergy);
1070         //UserNtuples->fillg61(eta_hit,1.);
1071         //UserNtuples->fillg62(phigrad_hit,losenergy);
1072         //UserNtuples->fillg63(phigrad_hit,1.);
1073 
1074         //    double   xx    = hitPoint.x();
1075         //    double   yy    = hitPoint.y();
1076         double zz = hitPoint.z();
1077 
1078         theHistManager->getHisto("zHitsnoMI")->Fill(zz);
1079 
1080         if (verbosity > 2) {
1081           edm::LogVerbatim("FP420Test") << "FP420Test:zHits = " << zz;
1082         }
1083         //   double   rr    = hitPoint.perp();
1084         /*
1085       if(aHit->getTrackID() == 1) {
1086       emu += aHit->getEnergyDeposit();} else {
1087       erest += aHit->getEnergyDeposit();}
1088     */
1089 
1090         //collect lost in Si en.of hits in every plane and put it into themap[]
1091         //UserNtuples->fillg30(losenergy,1.);
1092         themap[unitID] += losenergy;
1093         totallosenergy += losenergy;
1094 
1095         int det, zside, sector, zmodule;
1096         //    CaloNumberingPacker::unpackCastorIndex(unitID, det, zside, sector, zmodule);
1097         FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
1098         int justlayer = FP420NumberingScheme::unpackLayerIndex(rn00, zside);  // 1,2
1099         if (justlayer < 1 || justlayer > 2) {
1100           edm::LogVerbatim("FP420Test") << "FP420Test:WRONG  justlayer= " << justlayer;
1101         }
1102         // zside=1,2 ; zmodule=1,10 ; sector=1,3
1103         //UserNtuples->fillg44(float(sector),1.);
1104         //UserNtuples->fillg45(float(zmodule),1.);
1105         //UserNtuples->fillg46(float(zside),1.);
1106         //      int sScale = 20;
1107         // intindex is a continues numbering of FP420
1108         //int zScale = 2; unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside; //intindex=1-30:X,Y,X,Y,X,Y...
1109         // int zScale = 10;   unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule; //intindex=1-30:XXXXXXXXXX,YYYYYYYYYY,...
1110         //UserNtuples->fillg40(float(intindex),1.);
1111         //UserNtuples->fillg48(float(intindex),losenergy);
1112         //
1113         //=======================================
1114         //   G4ThreeVector middle = (hitExitLocalPoint+hitEntryLocalPoint)/2.;
1115         G4ThreeVector middle = (hitExitLocalPoint - hitEntryLocalPoint) / 2.;
1116         themapz[unitID] = hitPoint.z() + fabs(middle.z());
1117         if (verbosity > 2) {
1118           edm::LogVerbatim("FP420Test") << "1111111111111111111111111111111111111111111111111111111111111111111111111 ";
1119           edm::LogVerbatim("FP420Test") << "FP420Test: det, zside, sector, zmodule = " << det << zside << sector
1120                                         << zmodule;
1121           edm::LogVerbatim("FP420Test") << "FP420Test: justlayer = " << justlayer;
1122           edm::LogVerbatim("FP420Test") << "FP420Test: hitExitLocalPoint = " << hitExitLocalPoint;
1123           edm::LogVerbatim("FP420Test") << "FP420Test: hitEntryLocalPoint = " << hitEntryLocalPoint;
1124           edm::LogVerbatim("FP420Test") << "FP420Test:  middle= " << middle;
1125           edm::LogVerbatim("FP420Test") << "FP420Test:  hitPoint.z()-419000.= " << hitPoint.z() - 419000.;
1126 
1127           edm::LogVerbatim("FP420Test") << "FP420Test:zHits-419000. = " << themapz[unitID] - 419000.;
1128         }
1129         //=======================================
1130         // Y
1131         if (zside == 1) {
1132           //UserNtuples->fillg24(losenergy,1.);
1133           if (losenergy > 0.00003) {
1134             themap1[unitID] += 1.;
1135           }
1136         }
1137         //X
1138         if (zside == 2) {
1139           //UserNtuples->fillg25(losenergy,1.);
1140           if (losenergy > 0.00005) {
1141             themap1[unitID] += 1.;
1142           }
1143         }
1144         //     }
1145         //
1146         /*
1147         if (sector == 1) {
1148           nhit11 += 1;
1149           //UserNtuples->fillg33(rr,1.);
1150           //UserNtuples->fillg11(yy,1.);
1151         }
1152         if (sector == 2) {
1153           nhit12 += 1;
1154           //UserNtuples->fillg34(rr,1.);
1155           //UserNtuples->fillg86(yy,1.);
1156         }
1157         if (sector == 3) {
1158           nhit13 += 1;
1159           //UserNtuples->fillg35(rr,1.);
1160           //UserNtuples->fillg87(yy,1.);
1161         }
1162         */
1163         //UserNtuples->fillg10(xx,1.);
1164         //UserNtuples->fillg12(zz,1.);
1165         //UserNtuples->fillg32(rr,1.);
1166 
1167         // =========
1168         //    double xPrimAtZhit = vx + (zz-vz)*tan(th)*cos(phi);
1169         //    double yPrimAtZhit = vy + (zz-vz)*tan(th)*sin(phi);
1170 
1171         //       double  dx = xPrimAtZhit - xx;
1172         //       double  dy = yPrimAtZhit - yy;
1173 
1174         //                      Select SD hits
1175         //    if(rr<120.) {
1176         //                      Select MI or noMI over all 3 stations
1177         if (lastpo.z() < z4 && lastpo.perp() < 120.) {
1178           // MIonly:
1179           //UserNtuples->fillg16(lastpo.z(),1.);
1180           //UserNtuples->fillg18(zz,1.);
1181           //                                                                     Station I
1182           if (zz < z2) {
1183             //UserNtuples->fillg54(dx,1.);
1184             //UserNtuples->fillg55(dy,1.);
1185           }
1186           //                                                                     Station II
1187           if (zz < z3 && zz > z2) {
1188             //UserNtuples->fillg50(dx,1.);
1189             //UserNtuples->fillg51(dy,1.);
1190           }
1191           //                                                                     Station III
1192           if (zz < z4 && zz > z3) {
1193             //UserNtuples->fillg64(dx,1.);
1194             //UserNtuples->fillg65(dy,1.);
1195             //UserNtuples->filld209(xx,yy,1.);
1196           }
1197         } else {
1198           // no MIonly:
1199           //UserNtuples->fillg17(lastpo.z(),1.);
1200           //UserNtuples->fillg19(zz,1.);
1201           //UserNtuples->fillg74(incidentEnergyHit,1.);
1202           //UserNtuples->fillg75(float(trackIDhit),1.);
1203           //                                                                     Station I
1204           if (zz < z2) {
1205             //UserNtuples->fillg56(dx,1.);
1206             //UserNtuples->fillg57(dy,1.);
1207             //UserNtuples->fillg20(numofpart,1.);
1208             //UserNtuples->fillg21(sumEnerDeposit,1.);
1209             if (zside == 1) {
1210               //UserNtuples->fillg26(losenergy,1.);
1211             }
1212             if (zside == 2) {
1213               //UserNtuples->fillg76(losenergy,1.);
1214             }
1215             if (trackIDhit == 1) {
1216               //UserNtuples->fillg70(dx,1.);
1217               //UserNtuples->fillg71(incidentEnergyHit,1.);
1218               //UserNtuples->fillg79(losenergy,1.);
1219             } else {
1220               //UserNtuples->fillg82(dx,1.);
1221             }
1222           }
1223           //                                                                     Station II
1224           if (zz < z3 && zz > z2) {
1225             //UserNtuples->fillg52(dx,1.);
1226             //UserNtuples->fillg53(dy,1.);
1227             //UserNtuples->fillg22(numofpart,1.);
1228             //UserNtuples->fillg23(sumEnerDeposit,1.);
1229             //UserNtuples->fillg80(incidentEnergyHit,1.);
1230             //UserNtuples->fillg81(float(trackIDhit),1.);
1231             if (zside == 1) {
1232               //UserNtuples->fillg27(losenergy,1.);
1233             }
1234             if (zside == 2) {
1235               //UserNtuples->fillg77(losenergy,1.);
1236             }
1237             if (trackIDhit == 1) {
1238               //UserNtuples->fillg72(dx,1.);
1239               //UserNtuples->fillg73(incidentEnergyHit,1.);
1240             } else {
1241               //UserNtuples->fillg83(dx,1.);
1242             }
1243           }
1244           //                                                                     Station III
1245           if (zz < z4 && zz > z3) {
1246             //UserNtuples->fillg68(dx,1.);
1247             //UserNtuples->fillg69(dy,1.);
1248             //UserNtuples->filld210(xx,yy,1.);
1249             //UserNtuples->fillg22(numofpart,1.);
1250             //UserNtuples->fillg23(sumEnerDeposit,1.);
1251             if (zside == 1) {
1252               //UserNtuples->fillg28(losenergy,1.);
1253             }
1254             if (zside == 2) {
1255               //UserNtuples->fillg78(losenergy,1.);
1256             }
1257           }
1258         }  // MIonly or noMIonly ENDED
1259            //    }
1260 
1261         //     !!!!!!!!!!!!!
1262 
1263       }  // for loop on all hits ENDED  ENDED  ENDED  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1264 
1265       //     !!!!!!!!!!!!!
1266 
1267       //======================================================================================================SUMHIT
1268       //UserNtuples->fillg29(totallosenergy,1.);
1269       //UserNtuples->fillg36(nhit11,1.);
1270       //UserNtuples->fillg37(nhit12,1.);
1271       //UserNtuples->fillg38(nhit13,1.);
1272       //======================================================================================================SUMHIT
1273       //   int rn00=3;//test only with 2 sensors in superlayer, not 4
1274       //      int rn00=rn0;//always
1275       if (verbosity > 2) {
1276         edm::LogVerbatim("FP420Test") << "22222222222222222222222222222222222222222222222222222222222222222222222222 ";
1277       }
1278       int det = 1;
1279       int allplacesforsensors = 7;
1280       for (int sector = 1; sector < sn0; sector++) {
1281         for (int zmodule = 1; zmodule < pn0; zmodule++) {
1282           for (int zsideinorder = 1; zsideinorder < allplacesforsensors; zsideinorder++) {
1283             int zside = FP420NumberingScheme::realzside(rn00, zsideinorder);  //1,3,5,2,4,6
1284             if (verbosity > 2) {
1285               edm::LogVerbatim("FP420Test") << "FP420Test:  sector= " << sector << " zmodule= " << zmodule
1286                                             << " zsideinorder= " << zsideinorder << " zside= " << zside;
1287             }
1288             if (zside != 0) {
1289               int justlayer = FP420NumberingScheme::unpackLayerIndex(rn00, zside);  // 1,2
1290               if (justlayer < 1 || justlayer > 2) {
1291                 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG  justlayer= " << justlayer;
1292               }
1293               int copyinlayer = FP420NumberingScheme::unpackCopyIndex(rn00, zside);  // 1,2,3
1294               if (copyinlayer < 1 || copyinlayer > 3) {
1295                 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG  copyinlayer= " << copyinlayer;
1296               }
1297               int orientation = FP420NumberingScheme::unpackOrientation(rn00, zside);  // Front: = 1; Back: = 2
1298               if (orientation < 1 || orientation > 2) {
1299                 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG  orientation= " << orientation;
1300               }
1301 
1302               // iu is a continues numbering of planes(!)  over two arm FP420 set up
1303               int detfixed =
1304                   1;  // use this treatment for each set up arm, hence no sense to do it defferently for +FP420 and -FP420;
1305               //                                                                    and  ...[ii] massives have prepared in such a way
1306               unsigned int ii =
1307                   FP420NumberingScheme::packMYIndex(rn00, pn0, sn0, detfixed, justlayer, sector, zmodule) - 1;
1308               // ii = 0-19   --> 20 items
1309               if (verbosity > 2) {
1310                 edm::LogVerbatim("FP420Test")
1311                     << "FP420Test:  justlayer = " << justlayer << " copyinlayer = " << copyinlayer
1312                     << " orientation = " << orientation << " ii= " << ii;
1313               }
1314               double zdiststat = 0.;
1315               if (sn0 < 4) {
1316                 if (sector == 2)
1317                   zdiststat = zD3;
1318               } else {
1319                 if (sector == 2)
1320                   zdiststat = zD2;
1321                 if (sector == 3)
1322                   zdiststat = zD3;
1323               }
1324               double kplane = -(pn0 - 1) / 2 - 0.5 + (zmodule - 1);  //-3.5 +0...5 = -3.5,-2.5,-1.5,+2.5,+1.5
1325               double zcurrent = zinibeg + z420 + (zSiStep - zSiPlane) / 2 + kplane * zSiStep + zdiststat;
1326               //double zcurrent = zinibeg +(zSiStep-zSiPlane)/2  + kplane*zSiStep + (sector-1)*zUnit;
1327               if (verbosity > 2) {
1328                 edm::LogVerbatim("FP420Test") << "FP420Test:  Leftzcurrent-419000. = " << zcurrent - 419000.;
1329                 edm::LogVerbatim("FP420Test") << "FP420Test:  zGapLDet = " << zGapLDet;
1330               }
1331               if (justlayer == 1) {
1332                 if (orientation == 1)
1333                   zcurrent += (zGapLDet + zSiDet / 2);
1334                 if (orientation == 2)
1335                   zcurrent += zBlade - (zGapLDet + zSiDet / 2);
1336               }
1337               if (justlayer == 2) {
1338                 if (orientation == 1)
1339                   zcurrent += (zGapLDet + zSiDet / 2) + zBlade + gapBlade;
1340                 if (orientation == 2)
1341                   zcurrent += 2 * zBlade + gapBlade - (zGapLDet + zSiDet / 2);
1342               }
1343               //   .
1344               //
1345               if (det == 2)
1346                 zcurrent = -zcurrent;
1347               //
1348               if (verbosity > 2) {
1349                 edm::LogVerbatim("FP420Test") << "FP420Test:  zcurrent-419000. = " << zcurrent - 419000.;
1350               }
1351               //================================== end of for loops in continuius number iu:
1352             }  //if(zside!=0
1353           }    // for superlayer
1354         }      // for zmodule
1355       }        // for sector
1356 
1357       if (verbosity > 2) {
1358         edm::LogVerbatim("FP420Test")
1359             << "----------------------------------------------------------------------------- ";
1360       }
1361 
1362       //======================================================================================================CHECK
1363       if (totallosenergy == 0.0) {
1364         edm::LogVerbatim("FP420Test") << "FP420Test:     number of hits = " << theCAFI->entries();
1365         for (int j = 0; j < nhits; j++) {
1366           FP420G4Hit* aHit = (*theCAFI)[j];
1367           double losenergy = aHit->getEnergyLoss();
1368           edm::LogVerbatim("FP420Test") << " j hits = " << j << "losenergy = " << losenergy;
1369         }
1370       }
1371       //======================================================================================================CHECK
1372 
1373       //====================================================================================================== HIT  START
1374 
1375       //   FIBRE Hit collected analysis
1376       /*
1377       double totalEnergy = 0.;
1378       int  nhitsX = 0, nhitsY = 0, nsumhit = 0;
1379       for (int sector = 1; sector < 4; sector++) {
1380         int nhitsecX = 0, nhitsecY = 0;
1381         for (int zmodule = 1; zmodule < 11; zmodule++) {
1382           for (int zside = 1; zside < 3; zside++) {
1383             int det = 1;
1384             //      int nhit = 0;
1385             //  int sScale = 20;
1386             int index = FP420NumberingScheme::packFP420Index(det, zside, sector, zmodule);
1387             double theTotalEnergy = themap[index];
1388             //   X planes
1389             if (zside < 2) {
1390               //UserNtuples->fillg47(theTotalEnergy,1.);
1391               if (theTotalEnergy > 0.00003) {
1392                 nhitsX += 1;
1393                 //      nhitsecX += themap1[index];
1394                 //      nhit=1;
1395               }
1396             }
1397             //   Y planes
1398             else {
1399               //UserNtuples->fillg49(theTotalEnergy,1.);
1400               if (theTotalEnergy > 0.00005) {
1401                 nhitsY += 1;
1402                 //      nhitsecY += themap1[index];
1403                 //      nhit=1;
1404               }
1405             }
1406             // intindex is a continues numbering of FP420
1407             //        int zScale=2;  unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
1408             // int zScale=10;       unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule;
1409             //UserNtuples->fillg41(float(intindex),theTotalEnergy);
1410             //UserNtuples->fillg42(float(intindex),1.);
1411             //UserNtuples->fillp208(float(intindex),float(nhit),1.);
1412             //UserNtuples->fillp211(float(intindex),float(themap1[index]),1.);
1413             totalEnergy += themap[index];
1414           }  // for
1415         }    // for
1416         //UserNtuples->fillg39(nhitsecY,1.);
1417         if (nhitsecX > 10 || nhitsecY > 10) {
1418           nsumhit += 1;
1419           //UserNtuples->fillp213(float(sector),float(1.),1.);
1420         } else {  //UserNtuples->fillp213(float(sector),float(0.),1.);
1421         }
1422       }  // for
1423       //====================================================================================================== HIT  END
1424 
1425       //====================================================================================================== HIT  ALL
1426       //UserNtuples->fillg43(totalEnergy,1.);
1427       //UserNtuples->fillg58(nhitsX,1.);
1428       //UserNtuples->fillg59(nhitsY,1.);
1429       //  if( nsumhit !=0 ) { //UserNtuples->fillp212(vy,float(1.),1.);
1430       if (nsumhit >= 2) {  //UserNtuples->fillp212(vy,float(1.),1.);
1431       } else {             //UserNtuples->fillp212(vy,float(0.),1.);
1432       }
1433       */
1434       //====================================================================================================== HIT  ALL
1435 
1436       //====================================================================================================== number of hits
1437       // .............
1438       //    } // number of hits < 50
1439       // .............
1440     }  // MI or no MI or all  - end
1441 
1442   }  // primary end
1443      //=========================== thePrim != 0  end   ===
1444   // ==========================================================================
1445   if (verbosity > 0) {
1446     edm::LogVerbatim("FP420Test") << "FP420Test:  END OF Event " << (*evt)()->GetEventID();
1447   }
1448 }
1449 
1450 DEFINE_SIMWATCHER(FP420Test);  //=
[ 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!