File indexing completed on 2021-08-12 23:11:55
0001
0002
0003
0004
0005 #include <iostream>
0006 #include <map>
0007 #include <string>
0008
0009
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
0031 #include "SimG4CMS/FP420/interface/FP420NumberingScheme.h"
0032 #include "SimG4CMS/FP420/interface/FP420G4HitCollection.h"
0033
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
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
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
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
0123
0124 int iev;
0125 int itrk;
0126 G4double entot0, tracklength0;
0127
0128 double rinCalo, zinCalo;
0129 int lastTrackID;
0130 int verbosity;
0131
0132
0133 G4double sumEnerDeposit, sumStepl, sumStepc;
0134
0135 int numofpart;
0136
0137 G4ThreeVector lastpo;
0138
0139
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;
0153 std::string fDataLabel;
0154 std::string fOutputFile;
0155 std::string fRecreateFile;
0156 };
0157
0158
0159 enum ntfp420_elements { ntfp420_evt };
0160
0161 FP420Test::FP420Test(const edm::ParameterSet& p) {
0162
0163 edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("FP420Test");
0164 verbosity = m_Anal.getParameter<int>("Verbosity");
0165
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
0176
0177 pn0 = 6;
0178 sn0 = 3;
0179 rn00 = 7;
0180
0181 z420 = 420000.0;
0182 zD2 = 4000.0;
0183 zD3 = 8000.0;
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
0201 double zSiStation = (pn0 - 1) * (2 * (zBlade + gapBlade) + zKapton) + 2 * 6. + 0.0;
0202
0203 double eee1 = 11.;
0204 double eee2 = 12.;
0205
0206 zinibeg = (eee1 - eee2) / 2.;
0207
0208 z1 = zinibeg + (zSiStation + 10.) / 2 + z420;
0209 z2 = z1 + zD2;
0210 z3 = z1 + zD3;
0211 z4 = z1 + 2 * zD3;
0212
0213
0214 fp420eventntuple = new TNtuple("NTfp420event", "NTfp420event", "evt");
0215
0216 whichevent = 0;
0217
0218
0219
0220
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
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
0242
0243
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
0252
0253
0254 Fp420AnalysisHistManager::Fp420AnalysisHistManager(const TString& managername) {
0255
0256
0257 fTypeTitle = managername;
0258 fHistArray = new TObjArray();
0259 fHistNamesArray = new TObjArray();
0260
0261 TH1::AddDirectory(kFALSE);
0262 bookHistos();
0263
0264 fHistArray->Compress();
0265 fHistNamesArray->Compress();
0266
0267
0268 }
0269
0270
0271 Fp420AnalysisHistManager::~Fp420AnalysisHistManager() {
0272
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
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
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
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
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
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
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
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
0385
0386 Int_t index = fHistNamesArray->IndexOf(&histname);
0387 return getHisto(index);
0388 }
0389
0390
0391 TH2F* Fp420AnalysisHistManager::getHisto2(const TObjString& histname) {
0392
0393
0394 Int_t index = fHistNamesArray->IndexOf(&histname);
0395 return getHisto2(index);
0396 }
0397
0398
0399 void Fp420AnalysisHistManager::storeWeights() {
0400
0401
0402 for (int i = 0; i < fHistArray->GetEntries(); i++) {
0403 ((TH1F*)(fHistArray->At(i)))->Sumw2();
0404 }
0405 }
0406
0407
0408
0409
0410
0411
0412
0413
0414 void FP420Test::update(const BeginOfJob* job) {
0415
0416 edm::LogVerbatim("FP420Test") << "FP420Test:beggining of job";
0417 ;
0418 }
0419
0420
0421 void FP420Test::update(const BeginOfRun* run) {
0422
0423
0424 edm::LogVerbatim("FP420Test") << std::endl << "FP420Test:: Begining of Run";
0425 }
0426
0427 void FP420Test::update(const EndOfRun* run) { ; }
0428
0429
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
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
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();
0461
0462 theHistManager->getHisto("SumEDep")->Fill(sumEnerDeposit);
0463 theHistManager->getHisto("TrackL")->Fill(tracklength);
0464
0465
0466 G4ThreeVector vert_mom = (*trk)()->GetVertexMomentumDirection();
0467 G4ThreeVector vert_pos = (*trk)()->GetVertexPosition();
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496 if (tracklength < z4) {
0497
0498
0499
0500
0501 }
0502
0503
0504 const G4Step* aStep = (*trk)()->GetStep();
0505
0506
0507
0508 G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0509 lastpo = preStepPoint->GetPosition();
0510
0511
0512 if (lastpo.z() < z1 && lastpo.perp() < 100.) {
0513
0514
0515
0516
0517 }
0518 if ((lastpo.z() > z1 && lastpo.z() < z2) && lastpo.perp() < 100.) {
0519
0520
0521
0522
0523 }
0524 if (lastpo.z() < z2 && lastpo.perp() < 100.) {
0525
0526
0527
0528
0529
0530
0531
0532 }
0533 if (lastpo.z() < z3 && lastpo.perp() < 100.) {
0534
0535
0536
0537
0538 }
0539 if (lastpo.z() < z4 && lastpo.perp() < 100.) {
0540
0541
0542
0543
0544
0545
0546
0547 }
0548 }
0549 }
0550
0551
0552
0553
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
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();
0572 G4ThreeVector trackmom = theTrack->GetMomentum();
0573 G4double entot = theTrack->GetTotalEnergy();
0574 G4int curstepnumber = theTrack->GetCurrentStepNumber();
0575
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589
0590
0591 G4double stepl = aStep->GetStepLength();
0592 G4double EnerDeposit = aStep->GetTotalEnergyDeposit();
0593
0594
0595
0596
0597
0598
0599
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
0620
0621 if (id == 1) {
0622
0623 if (curstepnumber == 1) {
0624 entot0 = entot;
0625
0626 }
0627
0628
0629
0630
0631 if (prename == "SBST") {
0632 sumStepc += stepl;
0633
0634 }
0635
0636
0637 if (prename == "SBSTs") {
0638 sumStepl += stepl;
0639
0640 }
0641
0642
0643
0644
0645 if (trackstatus != 2) {
0646
0647 if (prename == "SIDETL" || prename == "SIDETR") {
0648 if (prename == "SIDETL") {
0649
0650 }
0651 if (prename == "SIDETR") {
0652
0653 }
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676 G4String posname = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
0677 if ((prename == "SIDETL" && posname != "SIDETL") || (prename == "SIDETR" && posname != "SIDETR")) {
0678 if (name1[2] == "SISTATION") {
0679
0680 }
0681 if (name1[3] == "SIPLANE") {
0682
0683 }
0684
0685 if (prename == "SIDETL") {
0686
0687
0688 if (copyno1[2] < 2) {
0689
0690 } else if (copyno1[2] < 3) {
0691
0692 if (copyno1[3] < 2) {
0693 } else if (copyno1[3] < 3) {
0694
0695 } else if (copyno1[3] < 4) {
0696
0697 } else if (copyno1[3] < 5) {
0698
0699 } else if (copyno1[3] < 6) {
0700
0701 } else if (copyno1[3] < 7) {
0702
0703 } else if (copyno1[3] < 8) {
0704
0705 } else if (copyno1[3] < 9) {
0706
0707 } else if (copyno1[3] < 10) {
0708
0709 }
0710 } else if (copyno1[2] < 4) {
0711
0712 } else if (copyno1[2] < 5) {
0713
0714 }
0715 }
0716 if (prename == "SIDETR") {
0717
0718
0719 if (copyno1[2] < 2) {
0720
0721 } else if (copyno1[2] < 3) {
0722
0723 if (copyno1[3] < 2) {
0724 } else if (copyno1[3] < 3) {
0725
0726 } else if (copyno1[3] < 4) {
0727
0728 } else if (copyno1[3] < 5) {
0729
0730 } else if (copyno1[3] < 6) {
0731
0732 } else if (copyno1[3] < 7) {
0733
0734 } else if (copyno1[3] < 8) {
0735
0736 } else if (copyno1[3] < 9) {
0737
0738 } else if (copyno1[3] < 10) {
0739
0740 }
0741 } else if (copyno1[2] < 4) {
0742
0743 } else if (copyno1[2] < 5) {
0744
0745 }
0746 }
0747 }
0748 }
0749
0750 }
0751
0752
0753
0754
0755
0756 sumEnerDeposit += EnerDeposit;
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773 if (trackstatus == 2) {
0774
0775
0776 tracklength0 = tracklength;
0777
0778
0779
0780
0781
0782 }
0783 }
0784
0785
0786 if (parentID == 1 && curstepnumber == 1) {
0787
0788 numofpart += 1;
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804 if (prename == "SBST") {
0805
0806 }
0807 if (prename == "SBSTs") {
0808
0809 }
0810 }
0811
0812
0813 }
0814
0815
0816 int FP420Test::detLevels(const G4VTouchable* touch) const {
0817
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
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
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
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
0860 fp420eventarray[ntfp420_evt] = (float)whichevent;
0861
0862
0863 int trackID = 0;
0864 G4PrimaryParticle* thePrim = nullptr;
0865
0866
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
0882 if (thePrim == nullptr)
0883 thePrim = avertex->GetPrimary(trackID);
0884
0885 if (thePrim != nullptr) {
0886
0887 G4double vx = 0., vy = 0., vz = 0.;
0888 vx = avertex->GetX0();
0889 vy = avertex->GetY0();
0890 vz = avertex->GetZ0();
0891
0892
0893
0894 theHistManager->getHisto("VtxX")->Fill(vx);
0895 theHistManager->getHisto("VtxY")->Fill(vy);
0896 theHistManager->getHisto("VtxZ")->Fill(vz);
0897 }
0898 }
0899
0900
0901
0902 if (thePrim != nullptr) {
0903
0904
0905
0906
0907
0908
0909
0910
0911
0912
0913
0914
0915
0916
0917
0918 if (lastpo.z() < z4 && lastpo.perp() < 100.) {
0919
0920 }
0921
0922
0923
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
0934
0935
0936
0937
0938
0939
0940
0941
0942
0943
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
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
0970
0971
0972 G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
0973
0974 if (verbosity > 0) {
0975 edm::LogVerbatim("FP420Test") << "FP420Test: accessed all HC";
0976 ;
0977 }
0978 int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("FP420SI");
0979
0980
0981
0982 FP420G4HitCollection* theCAFI = (FP420G4HitCollection*)allHC->GetHC(CAFIid);
0983
0984 if (verbosity > 0) {
0985
0986 edm::LogVerbatim("FP420Test") << "FP420Test: theCAFI->entries = " << theCAFI->entries();
0987 }
0988
0989 theHistManager->getHisto("NumberOfHits")->Fill(theCAFI->entries());
0990
0991
0992
0993
0994
0995
0996
0997 int varia;
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007 if (lastpo.z() < z4) {
1008
1009
1010
1011 varia = 1;
1012 } else {
1013
1014
1015
1016 varia = 2;
1017 }
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
1028
1029 if (varia == 2) {
1030
1031
1032
1033
1034
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
1045
1046
1047 int trackIDhit = aHit->getTrackID();
1048 unsigned int unitID = aHit->getUnitID();
1049
1050
1051
1052
1053
1054 double losenergy = aHit->getEnergyLoss();
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
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
1086
1087
1088
1089
1090
1091
1092
1093
1094 themap[unitID] += losenergy;
1095 totallosenergy += losenergy;
1096
1097 int det, zside, sector, zmodule;
1098
1099 FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
1100 int justlayer = FP420NumberingScheme::unpackLayerIndex(rn00, zside);
1101 if (justlayer < 1 || justlayer > 2) {
1102 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG justlayer= " << justlayer;
1103 }
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
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
1133 if (zside == 1) {
1134
1135 if (losenergy > 0.00003) {
1136 themap1[unitID] += 1.;
1137 }
1138 }
1139
1140 if (zside == 2) {
1141
1142 if (losenergy > 0.00005) {
1143 themap1[unitID] += 1.;
1144 }
1145 }
1146
1147
1148 if (sector == 1) {
1149 nhit11 += 1;
1150
1151
1152 }
1153 if (sector == 2) {
1154 nhit12 += 1;
1155
1156
1157 }
1158 if (sector == 3) {
1159 nhit13 += 1;
1160
1161
1162 }
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177 if (lastpo.z() < z4 && lastpo.perp() < 120.) {
1178
1179
1180
1181
1182 if (zz < z2) {
1183
1184
1185 }
1186
1187 if (zz < z3 && zz > z2) {
1188
1189
1190 }
1191
1192 if (zz < z4 && zz > z3) {
1193
1194
1195
1196 }
1197 } else {
1198
1199
1200
1201
1202
1203
1204 if (zz < z2) {
1205
1206
1207
1208
1209 if (zside == 1) {
1210
1211 }
1212 if (zside == 2) {
1213
1214 }
1215 if (trackIDhit == 1) {
1216
1217
1218
1219 } else {
1220
1221 }
1222 }
1223
1224 if (zz < z3 && zz > z2) {
1225
1226
1227
1228
1229
1230
1231 if (zside == 1) {
1232
1233 }
1234 if (zside == 2) {
1235
1236 }
1237 if (trackIDhit == 1) {
1238
1239
1240 } else {
1241
1242 }
1243 }
1244
1245 if (zz < z4 && zz > z3) {
1246
1247
1248
1249
1250
1251 if (zside == 1) {
1252
1253 }
1254 if (zside == 2) {
1255
1256 }
1257 }
1258 }
1259
1260
1261
1262
1263 }
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
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);
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);
1290 if (justlayer < 1 || justlayer > 2) {
1291 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG justlayer= " << justlayer;
1292 }
1293 int copyinlayer = FP420NumberingScheme::unpackCopyIndex(rn00, zside);
1294 if (copyinlayer < 1 || copyinlayer > 3) {
1295 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG copyinlayer= " << copyinlayer;
1296 }
1297 int orientation = FP420NumberingScheme::unpackOrientation(rn00, zside);
1298 if (orientation < 1 || orientation > 2) {
1299 edm::LogVerbatim("FP420Test") << "FP420Test:WRONG orientation= " << orientation;
1300 }
1301
1302
1303 int detfixed =
1304 1;
1305
1306 unsigned int ii =
1307 FP420NumberingScheme::packMYIndex(rn00, pn0, sn0, detfixed, justlayer, sector, zmodule) - 1;
1308
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);
1325 double zcurrent = zinibeg + z420 + (zSiStep - zSiPlane) / 2 + kplane * zSiStep + zdiststat;
1326
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
1352 }
1353 }
1354 }
1355 }
1356
1357 if (verbosity > 2) {
1358 edm::LogVerbatim("FP420Test")
1359 << "----------------------------------------------------------------------------- ";
1360 }
1361
1362
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
1372
1373
1374
1375
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
1384
1385 int index = FP420NumberingScheme::packFP420Index(det, zside, sector, zmodule);
1386 double theTotalEnergy = themap[index];
1387
1388 if (zside < 2) {
1389
1390 if (theTotalEnergy > 0.00003) {
1391 nhitsX += 1;
1392
1393
1394 }
1395 }
1396
1397 else {
1398
1399 if (theTotalEnergy > 0.00005) {
1400 nhitsY += 1;
1401
1402
1403 }
1404 }
1405
1406
1407
1408
1409
1410
1411
1412 totalEnergy += themap[index];
1413 }
1414 }
1415
1416 if (nhitsecX > 10 || nhitsecY > 10) {
1417 nsumhit += 1;
1418
1419 } else {
1420 }
1421 }
1422
1423
1424
1425
1426
1427
1428
1429 if (nsumhit >= 2) {
1430 } else {
1431 }
1432
1433
1434
1435
1436
1437
1438
1439 }
1440
1441 }
1442
1443
1444 if (verbosity > 0) {
1445 edm::LogVerbatim("FP420Test") << "FP420Test: END OF Event " << (*evt)()->GetEventID();
1446 }
1447 }
1448
1449 DEFINE_SIMWATCHER(FP420Test);