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