Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-12 23:11:55 

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