CMS 3D CMS Logo
FP420Test.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 
4 // system include files
5 #include <iostream>
6 #include <iomanip>
7 #include <cmath>
8 #include <vector>
9 //
10 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
11 #include "SimG4Core/Notification/interface/EndOfEvent.h"
12 #include "SimG4Core/Notification/interface/TrackWithHistory.h"
13 #include "SimG4Core/Notification/interface/TrackInformation.h"
14 
15 #include "FWCore/Framework/interface/Event.h"
16 #include "FWCore/Framework/interface/EventSetup.h"
17 #include "FWCore/MessageLogger/interface/MessageLogger.h"
18 // to retreive hits
19 #include "SimG4CMS/FP420/interface/FP420NumberingScheme.h"
20 #include "SimG4CMS/FP420/interface/FP420G4HitCollection.h"
21 #include "SimG4CMS/FP420/interface/FP420Test.h"
22 
23 // G4 stuff
24 #include "G4SDManager.hh"
25 #include "G4Step.hh"
26 #include "G4Track.hh"
27 #include "G4VProcess.hh"
28 #include "G4HCofThisEvent.hh"
29 #include "G4UserEventAction.hh"
30 #include "G4TransportationManager.hh"
31 #include "G4ProcessManager.hh"
32 
33 #include "CLHEP/Units/GlobalSystemOfUnits.h"
34 #include "CLHEP/Units/GlobalPhysicalConstants.h"
35 
36 //================================================================
37 // Root stuff
38 
39 // Include the standard header <cassert> to effectively include
40 // the standard header <assert.h> within the std namespace.
41 #include <cassert>
42 
43 using namespace edm;
44 using namespace std;
45 //================================================================
46 
47 //UserVerbosity FP420Test::std::cout("FP420Test","info","FP420Test");
48 
49 enum ntfp420_elements { ntfp420_evt };
50 //================================================================
51 FP420Test::FP420Test(const edm::ParameterSet& p) {
52  //constructor
53  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("FP420Test");
54  verbosity = m_Anal.getParameter<int>("Verbosity");
55  //verbosity = 1;
56 
57  fDataLabel = m_Anal.getParameter<std::string>("FDataLabel");
58  fOutputFile = m_Anal.getParameter<std::string>("FOutputFile");
59  fRecreateFile = m_Anal.getParameter<std::string>("FRecreateFile");
60 
61  if (verbosity > 0) {
62  std::cout << "============================================================================" << std::endl;
63  std::cout << "FP420Testconstructor :: Initialized as observer" << std::endl;
64  }
65  // Initialization:
66 
67  pn0 = 6;
68  sn0 = 3;
69  rn00 = 7;
70 
71  z420 = 420000.0; // mm
72  zD2 = 4000.0; // mm
73  zD3 = 8000.0; // mm
74  //
75  zBlade = 5.00;
76  gapBlade = 1.6;
77  double gapSupplane = 1.6;
78  ZSiPlane = 2 * zBlade + gapBlade + gapSupplane;
79 
80  double ZKapton = 0.1;
81  ZSiStep = ZSiPlane + ZKapton;
82 
83  double ZBoundDet = 0.020;
84  double ZSiElectr = 0.250;
85  double ZCeramDet = 0.500;
86  //
87  ZSiDet = 0.250;
88  ZGapLDet = zBlade / 2 - (ZSiDet + ZSiElectr + ZBoundDet + ZCeramDet / 2);
89  //
90  // ZSiStation = 5*(2*(5.+1.6)+0.1)+2*6.+1.0 = 79.5
91  double ZSiStation = (pn0 - 1) * (2 * (zBlade + gapBlade) + ZKapton) + 2 * 6. + 0.0; // = 78.5
92  // 11.=e1, 12.=e2 in zzzrectangle.xml
93  double eee1 = 11.;
94  double eee2 = 12.;
95 
96  zinibeg = (eee1 - eee2) / 2.;
97 
98  z1 = zinibeg + (ZSiStation + 10.) / 2 + z420; // z1 - right after 1st Station
99  z2 = z1 + zD2; //z2 - right after middle Station
100  z3 = z1 + zD3; //z3 - right after last Station
101  z4 = z1 + 2 * zD3;
102  //==================================
103 
104  fp420eventntuple = new TNtuple("NTfp420event", "NTfp420event", "evt");
105 
106  whichevent = 0;
107 
108  // fDataLabel = "defaultData";
109  // fOutputFile = "TheAnlysis.root";
110  // fRecreateFile = "RECREATE";
111 
112  TheHistManager = new Fp420AnalysisHistManager(fDataLabel);
113 
114  if (verbosity > 0) {
115  std::cout << "FP420Test constructor :: Initialized Fp420AnalysisHistManager" << std::endl;
116  }
117 }
118 
119 FP420Test::~FP420Test() {
120  // delete UserNtuples;
121 
122  TFile fp420OutputFile("newntfp420.root", "RECREATE");
123  std::cout << "FP420 output root file has been created";
124  fp420eventntuple->Write();
125  std::cout << ", written";
126  fp420OutputFile.Close();
127  std::cout << ", closed";
128  delete fp420eventntuple;
129  std::cout << ", and deleted" << std::endl;
130 
131  //------->while end
132 
133  // Write histograms to file
134  TheHistManager->WriteToFile(fOutputFile, fRecreateFile);
135  if (verbosity > 0) {
136  std::cout << std::endl << "FP420Test Destructor --------> End of FP420Test : " << std::endl;
137  }
138 }
139 
140 //================================================================
141 // Histoes:
142 //-----------------------------------------------------------------------------
143 
144 Fp420AnalysisHistManager::Fp420AnalysisHistManager(const TString& managername) {
145  // The Constructor
146 
147  fTypeTitle = managername;
148  fHistArray = new TObjArray(); // Array to store histos
149  fHistNamesArray = new TObjArray(); // Array to store histos's names
150 
151  TH1::AddDirectory(kFALSE);
152  BookHistos();
153 
154  fHistArray->Compress(); // Removes empty space
155  fHistNamesArray->Compress();
156 
157  // StoreWeights(); // Store the weights
158 }
159 //-----------------------------------------------------------------------------
160 
161 Fp420AnalysisHistManager::~Fp420AnalysisHistManager() {
162  // The Destructor
163 
164  if (fHistArray) {
165  fHistArray->Delete();
166  delete fHistArray;
167  }
168 
169  if (fHistNamesArray) {
170  fHistNamesArray->Delete();
171  delete fHistNamesArray;
172  }
173 }
174 //-----------------------------------------------------------------------------
175 
176 void Fp420AnalysisHistManager::BookHistos() {
177  // at Start: (mm)
178  HistInit("PrimaryEta", "Primary Eta", 100, 9., 12.);
179  HistInit("PrimaryPhigrad", "Primary Phigrad", 100, 0., 360.);
180  HistInit("PrimaryTh", "Primary Th", 100, 0., 180.);
181  HistInit("PrimaryLastpoZ", "Primary Lastpo Z", 100, -200., 430000.);
182  HistInit("PrimaryLastpoX", "Primary Lastpo X Z<z4", 100, -30., 30.);
183  HistInit("PrimaryLastpoY", "Primary Lastpo Y Z<z4", 100, -30., 30.);
184  HistInit("XLastpoNumofpart", "Primary Lastpo X n>10", 100, -30., 30.);
185  HistInit("YLastpoNumofpart", "Primary Lastpo Y n>10", 100, -30., 30.);
186  HistInit("VtxX", "Vtx X", 100, -50., 50.);
187  HistInit("VtxY", "Vtx Y", 100, -50., 50.);
188  HistInit("VtxZ", "Vtx Z", 100, -200., 430000.);
189  // Book the histograms and add them to the array
190  HistInit("SumEDep", "This is sum Energy deposited", 100, -1, 199.);
191  HistInit("TrackL", "This is TrackL", 100, 0., 12000.);
192  HistInit("zHits", "z Hits all events", 100, 400000., 430000.);
193  HistInit("zHitsnoMI", "z Hits no MI", 100, 400000., 430000.);
194  HistInit("zHitsTrLoLe", "z Hits TrLength bigger 8300", 100, 400000., 430000.);
195  HistInit("NumberOfHits", "NumberOfHits", 100, 0., 300.);
196 }
197 
198 //-----------------------------------------------------------------------------
199 
200 void Fp420AnalysisHistManager::WriteToFile(const TString& fOutputFile, const TString& fRecreateFile) {
201  //Write to file = fOutputFile
202 
203  std::cout << "================================================================" << std::endl;
204  std::cout << " Write this Analysis to File " << fOutputFile << std::endl;
205  std::cout << "================================================================" << std::endl;
206 
207  TFile* file = new TFile(fOutputFile, fRecreateFile);
208 
209  fHistArray->Write();
210  file->Close();
211 }
212 //-----------------------------------------------------------------------------
213 
214 void Fp420AnalysisHistManager::HistInit(const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup) {
215  // Add histograms and histograms names to the array
216 
217  char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
218  strcpy(newtitle, title);
219  strcat(newtitle, " (");
220  strcat(newtitle, fTypeTitle);
221  strcat(newtitle, ") ");
222  fHistArray->AddLast((new TH1F(name, newtitle, nbinsx, xlow, xup)));
223  fHistNamesArray->AddLast(new TObjString(name));
224 }
225 //-----------------------------------------------------------------------------
226 
227 void Fp420AnalysisHistManager::HistInit(
228  const char* name, const char* title, Int_t nbinsx, Axis_t xlow, Axis_t xup, Int_t nbinsy, Axis_t ylow, Axis_t yup) {
229  // Add histograms and histograms names to the array
230 
231  char* newtitle = new char[strlen(title) + strlen(fTypeTitle) + 5];
232  strcpy(newtitle, title);
233  strcat(newtitle, " (");
234  strcat(newtitle, fTypeTitle);
235  strcat(newtitle, ") ");
236  fHistArray->AddLast((new TH2F(name, newtitle, nbinsx, xlow, xup, nbinsy, ylow, yup)));
237  fHistNamesArray->AddLast(new TObjString(name));
238 }
239 //-----------------------------------------------------------------------------
240 
241 TH1F* Fp420AnalysisHistManager::GetHisto(Int_t Number) {
242  // Get a histogram from the array with index = Number
243 
244  if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
245  return (TH1F*)(fHistArray->At(Number));
246 
247  } else {
248  std::cout << "!!!!!!!!!!!!!!!!!!GetHisto!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
249  std::cout << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number << std::endl;
250  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
251 
252  return (TH1F*)(fHistArray->At(0));
253  }
254 }
255 //-----------------------------------------------------------------------------
256 
257 TH2F* Fp420AnalysisHistManager::GetHisto2(Int_t Number) {
258  // Get a histogram from the array with index = Number
259 
260  if (Number <= fHistArray->GetLast() && fHistArray->At(Number) != (TObject*)nullptr) {
261  return (TH2F*)(fHistArray->At(Number));
262 
263  } else {
264  std::cout << "!!!!!!!!!!!!!!!!GetHisto2!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
265  std::cout << " WARNING ERROR - HIST ID INCORRECT (TOO HIGH) - " << Number << std::endl;
266  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
267 
268  return (TH2F*)(fHistArray->At(0));
269  }
270 }
271 //-----------------------------------------------------------------------------
272 
273 TH1F* Fp420AnalysisHistManager::GetHisto(const TObjString& histname) {
274  // Get a histogram from the array with name = histname
275 
276  Int_t index = fHistNamesArray->IndexOf(&histname);
277  return GetHisto(index);
278 }
279 //-----------------------------------------------------------------------------
280 
281 TH2F* Fp420AnalysisHistManager::GetHisto2(const TObjString& histname) {
282  // Get a histogram from the array with name = histname
283 
284  Int_t index = fHistNamesArray->IndexOf(&histname);
285  return GetHisto2(index);
286 }
287 //-----------------------------------------------------------------------------
288 
289 void Fp420AnalysisHistManager::StoreWeights() {
290  // Add structure to each histogram to store the weights
291 
292  for (int i = 0; i < fHistArray->GetEntries(); i++) {
293  ((TH1F*)(fHistArray->At(i)))->Sumw2();
294  }
295 }
296 
297 // Histoes end :
298 
299 //================================================================
300 
301 // using several observers
302 
303 //==================================================================== per JOB
304 void FP420Test::update(const BeginOfJob* job) {
305  //job
306  std::cout << "FP420Test:beggining of job" << std::endl;
307  ;
308 }
309 
310 //==================================================================== per RUN
311 void FP420Test::update(const BeginOfRun* run) {
312  //run
313 
314  std::cout << std::endl << "FP420Test:: Begining of Run" << std::endl;
315 }
316 
317 void FP420Test::update(const EndOfRun* run) { ; }
318 
319 //=================================================================== per EVENT
320 void FP420Test::update(const BeginOfEvent* evt) {
321  iev = (*evt)()->GetEventID();
322  if (verbosity > 0) {
323  std::cout << "FP420Test:update Event number = " << iev << std::endl;
324  }
325  whichevent++;
326 }
327 
328 //=================================================================== per Track
329 void FP420Test::update(const BeginOfTrack* trk) {
330  itrk = (*trk)()->GetTrackID();
331  if (verbosity > 1) {
332  std::cout << "FP420Test:update BeginOfTrack number = " << itrk << std::endl;
333  }
334  if (itrk == 1) {
335  SumEnerDeposit = 0.;
336  numofpart = 0;
337  SumStepl = 0.;
338  SumStepc = 0.;
339  tracklength0 = 0.;
340  }
341 }
342 
343 //=================================================================== per EndOfTrack
344 void FP420Test::update(const EndOfTrack* trk) {
345  itrk = (*trk)()->GetTrackID();
346  if (verbosity > 1) {
347  std::cout << "FP420Test:update EndOfTrack number = " << itrk << std::endl;
348  }
349  if (itrk == 1) {
350  G4double tracklength = (*trk)()->GetTrackLength(); // Accumulated track length
351 
352  TheHistManager->GetHisto("SumEDep")->Fill(SumEnerDeposit);
353  TheHistManager->GetHisto("TrackL")->Fill(tracklength);
354 
355  // direction !!!
356  G4ThreeVector vert_mom = (*trk)()->GetVertexMomentumDirection();
357  G4ThreeVector vert_pos = (*trk)()->GetVertexPosition(); // vertex ,where this track was created
358 
359  // float eta = 0.5 * log( (1.+vert_mom.z()) / (1.-vert_mom.z()) );
360  // float phi = atan2(vert_mom.y(),vert_mom.x());
361  // if (phi < 0.) phi += twopi;
362  // float phigrad = phi*180./pi;
363 
364  // float XV = vert_pos.x(); // mm
365  // float YV = vert_pos.y(); // mm
366  //UserNtuples->fillg543(phigrad,1.);
367  //UserNtuples->fillp203(phigrad,SumStepl,1.);
368  //UserNtuples->fillg544(XV,1.);
369  //UserNtuples->fillp201(XV,SumStepl,1.);
370  //UserNtuples->fillg545(YV,1.);
371  //UserNtuples->fillp202(YV,SumStepl,1.);
372 
373  //UserNtuples->fillg524(eta,1.);
374 
375  //UserNtuples->fillg534(SumStepl,1.);
376  //UserNtuples->fillg535(eta,SumStepl);
377  //UserNtuples->fillp207(eta,SumStepl,1.);
378  //UserNtuples->filld217(eta,SumStepl,1.);
379  //UserNtuples->filld220(phigrad,SumStepl,1.);
380  //UserNtuples->filld221(XV,SumStepl,1.);
381  //UserNtuples->filld222(YV,SumStepl,1.);
382  //UserNtuples->fillg537(SumStepc,1.);
383  //UserNtuples->fillg84(SumStepl,1.);
384 
385  // MI = (multiple interactions):
386  if (tracklength < z4) {
387  // //UserNtuples->fillp208(eta,tracklength,1.);
388  //UserNtuples->filld218(eta,tracklength,1.);
389  //UserNtuples->fillg538(SumStepc,1.);
390  //UserNtuples->fillg85(SumStepl,1.);
391  }
392 
393  // last step information
394  const G4Step* aStep = (*trk)()->GetStep();
395  // G4int csn = (*trk)()->GetCurrentStepNumber();
396  // G4double sl = (*trk)()->GetStepLength();
397  // preStep
398  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
399  lastpo = preStepPoint->GetPosition();
400 
401  // Analysis:
402  if (lastpo.z() < z1 && lastpo.perp() < 100.) {
403  //UserNtuples->fillg525(eta,1.);
404  //UserNtuples->fillg546(XV,1.);
405  //UserNtuples->fillg551(YV,1.);
406  //UserNtuples->fillg556(phigrad,1.);
407  }
408  if ((lastpo.z() > z1 && lastpo.z() < z2) && lastpo.perp() < 100.) {
409  //UserNtuples->fillg526(eta,1.);
410  //UserNtuples->fillg547(XV,1.);
411  //UserNtuples->fillg552(YV,1.);
412  //UserNtuples->fillg557(phigrad,1.);
413  }
414  if (lastpo.z() < z2 && lastpo.perp() < 100.) {
415  //UserNtuples->fillg527(eta,1.);
416  //UserNtuples->fillg548(XV,1.);
417  //UserNtuples->fillg553(YV,1.);
418  //UserNtuples->fillg558(phigrad,1.);
419  //UserNtuples->fillg521(lastpo.x(),1.);
420  //UserNtuples->fillg522(lastpo.y(),1.);
421  //UserNtuples->fillg523(lastpo.z(),1.);
422  }
423  if (lastpo.z() < z3 && lastpo.perp() < 100.) {
424  //UserNtuples->fillg528(eta,1.);
425  //UserNtuples->fillg549(XV,1.);
426  //UserNtuples->fillg554(YV,1.);
427  //UserNtuples->fillg559(phigrad,1.);
428  }
429  if (lastpo.z() < z4 && lastpo.perp() < 100.) {
430  //UserNtuples->fillg529(eta,1.);
431  //UserNtuples->fillg550(XV,1.);
432  //UserNtuples->fillg555(YV,1.);
433  //UserNtuples->fillg560(phigrad,1.);
434  //UserNtuples->fillg531(lastpo.x(),1.);
435  //UserNtuples->fillg532(lastpo.y(),1.);
436  //UserNtuples->fillg533(lastpo.z(),1.);
437  }
438  }
439 }
440 
441 // ====================================================
442 
443 //=================================================================== each STEP
444 void FP420Test::update(const G4Step* aStep) {
445  // ==========================================================================
446 
447  if (verbosity > 2) {
448  G4int stepnumber = aStep->GetTrack()->GetCurrentStepNumber();
449  std::cout << "FP420Test:update Step number = " << stepnumber << std::endl;
450  }
451  // track on aStep: !
452  G4Track* theTrack = aStep->GetTrack();
453  TrackInformation* trkInfo = dynamic_cast<TrackInformation*>(theTrack->GetUserInformation());
454  if (trkInfo == nullptr) {
455  std::cout << "FP420Test on aStep: No trk info !!!! abort " << std::endl;
456  }
457  G4int id = theTrack->GetTrackID();
458  G4String particleType = theTrack->GetDefinition()->GetParticleName(); // !!!
459  G4int parentID = theTrack->GetParentID(); // !!!
460  G4TrackStatus trackstatus = theTrack->GetTrackStatus(); // !!!
461  G4double tracklength = theTrack->GetTrackLength(); // Accumulated track length
462  G4ThreeVector trackmom = theTrack->GetMomentum();
463  G4double entot = theTrack->GetTotalEnergy(); // !!! deposited on step
464  G4int curstepnumber = theTrack->GetCurrentStepNumber();
465  // const G4ThreeVector& vert_pos = theTrack->GetVertexPosition(); // vertex ,where this track was created
466  // const G4ThreeVector& vert_mom = theTrack->GetVertexMomentumDirection();
467 
468  // double costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+vert_mom.y()*vert_mom.y()+vert_mom.z()*vert_mom.z());
469  // double theta = acos(min(max(costheta,double(-1.)),double(1.)));
470  // float eta = -log(tan(theta/2));
471  // double phi = -1000.;
472  // if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x());
473  // if (phi < 0.) phi += twopi;
474  // double phigrad = phi*360./twopi;
475 
476  //G4double trackvel = theTrack->GetVelocity();
477 
478  //std::cout << " trackstatus= " << trackstatus << " entot= " << entot << std::endl;
479 
480  // step points: !
481  G4double stepl = aStep->GetStepLength();
482  G4double EnerDeposit = aStep->GetTotalEnergyDeposit();
483  // pointers: !
484  //G4VPhysicalVolume* physvol = theTrack->GetVolume();
485  //G4VPhysicalVolume* nextphysvol = theTrack->GetNextVolume();
486  //G4Material* materialtr = theTrack->GetMaterial();
487  //G4Material* nextmaterialtr = theTrack->GetNextMaterial();
488 
489  // preStep
490  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
491  const G4ThreeVector& preposition = preStepPoint->GetPosition();
492  G4ThreeVector prelocalpoint = theTrack->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(preposition);
493  G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume();
494  const G4String& prename = currentPV->GetName();
495 
496  const G4VTouchable* pre_touch = preStepPoint->GetTouchable();
497  int pre_levels = detLevels(pre_touch);
498 
499  G4String name1[20];
500  int copyno1[20];
501  for (int i = 0; i < 20; ++i) {
502  name1[i] = "";
503  copyno1[i] = 0;
504  }
505  if (pre_levels > 0) {
506  detectorLevel(pre_touch, pre_levels, copyno1, name1);
507  }
508 
509  // G4LogicalVolume* lv = currentPV->GetLogicalVolume();
510  // G4Material* mat = lv->GetMaterial();
511  // std::string prenameVolume;
512  // prenameVolume.assign(prename,0,20);
513 
514  // G4double prebeta = preStepPoint->GetBeta();
515  // G4double precharge = preStepPoint->GetCharge();
516  // G4double prerad = mat->GetRadlen();
517 
518  // std::cout << " EneryDeposited = " << EnerDeposit << std::endl;
519  // std::cout << " prevolume = " << prename << std::endl;
521  // std::cout << " preposition = " << preposition << std::endl;
522  /*
523  // postStep
524  G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
525  G4ThreeVector posposition = postStepPoint->GetPosition();
526  G4ThreeVector poslocalpoint = theTrack->GetTouchable()->GetHistory()->
527  GetTopTransform().TransformPoint(posposition);
528  G4VPhysicalVolume* poscurrentPV = postStepPoint->GetPhysicalVolume();
529  G4String posname = poscurrentPV->GetName();
530 // G4LogicalVolume* poslv = poscurrentPV->GetLogicalVolume();
531 // G4Material* posmat = poslv->GetMaterial();
532 // std::string posnameVolume;
533 // posnameVolume.assign(posname,0,20);
534 
535 #ifdef ddebug
536  std::cout << "============posStep: information:============" << std::endl;
537  std::cout << " posposition = " << posposition
538  << " poslocalpoint = " << poslocalpoint
539  << " posvolume = " << posname
540  // << " posnameVolume = " << posnameVolume
541  << std::endl;
542 
543  std::cout << " ==========================================" << std::endl;
544 #endif
545 
546 
547 */
548 
549  // //
550  // //
551 
552  // 24.01.2006:
553  // tr : id parentID trackstatus tracklength curstepnumber entot vert_pos
554  // st : stepl EnerDeposit
555  // prestep: preposition prevolume = SBSTm SIDETL SIDETR name= SISTATION copy= 1,2,3 name= SIPLANE copy= 1..5..10
556 
557  // gen_track:
558  // id=1 parentID=1 trackstatus=0,2 tracklength(accumulated) curstepnumber entot vert_pos
559  if (id == 1) {
560  // on 1-st step:
561  if (curstepnumber == 1) {
562  entot0 = entot;
563  //UserNtuples->fillg519(entot0,1.);
564  }
565 
566  // on every step:
567 
568  // for Copper:
569  if (prename == "SBST") {
570  SumStepc += stepl;
571  // =========
572  }
573  // for ststeel:
574  // if(prename == "SBSTs") {
575  if (prename == "SBSTs") {
576  SumStepl += stepl;
577  // =========
578  }
579  // =========
580  // =========
581 
582  // exclude last track point if it is in SD (MI was started their)
583  if (trackstatus != 2) {
584  // for SD: Si Det.: SISTATION:SIPLANE:(SIDETL+BOUNDDET +SIDETR + CERAMDET)
585  if (prename == "SIDETL" || prename == "SIDETR") {
586  if (prename == "SIDETL") {
587  //UserNtuples->fillg569(EnerDeposit,1.);
588  }
589  if (prename == "SIDETR") {
590  //UserNtuples->fillg570(EnerDeposit,1.);
591  }
592 
593  // double numbcont = 10.*(copyno1[2]-1)+copyno1[3];
594 
595  // =========
596  // double xx = preposition.x();
597  // double yy = preposition.y();
598  // double zz = preposition.z();
599  // =========
600  //UserNtuples->fillg580(theta,1.);
601  //UserNtuples->fillg07(phigrad,1.);
602  // double xPrimAtZhit = vert_pos.x() + (zz-vert_pos.z())*tan(theta)*cos(phi);
603  // double yPrimAtZhit = vert_pos.y() + (zz-vert_pos.z())*tan(theta)*sin(phi);
604  // =========
605  // double dx = xPrimAtZhit - xx;
606  // double dy = yPrimAtZhit - yy;
607  // =========
608 
609  // //UserNtuples->fillp212(numbcont,dx,1.);
610  // //UserNtuples->fillp213(numbcont,dy,1.);
611  // =========
612 
613  // last step of current SD Volume:
614  G4String posname = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetName();
615  if ((prename == "SIDETL" && posname != "SIDETL") || (prename == "SIDETR" && posname != "SIDETR")) {
616  if (name1[2] == "SISTATION") {
617  //UserNtuples->fillg539(copyno1[2],1.);
618  }
619  if (name1[3] == "SIPLANE") {
620  //UserNtuples->fillg540(copyno1[3],1.);
621  }
622 
623  if (prename == "SIDETL") {
624  //UserNtuples->fillg541(EnerDeposit,1.);
625  //UserNtuples->fillg561(numbcont,1.);
626  if (copyno1[2] < 2) {
627  //UserNtuples->fillg571(dx,1.);
628  } else if (copyno1[2] < 3) {
629  //UserNtuples->fillg563(dx,1.);
630  if (copyno1[3] < 2) {
631  } else if (copyno1[3] < 3) {
632  //UserNtuples->fillg572(dx,1.);
633  } else if (copyno1[3] < 4) {
634  //UserNtuples->fillg573(dx,1.);
635  } else if (copyno1[3] < 5) {
636  //UserNtuples->fillg574(dx,1.);
637  } else if (copyno1[3] < 6) {
638  //UserNtuples->fillg575(dx,1.);
639  } else if (copyno1[3] < 7) {
640  //UserNtuples->fillg576(dx,1.);
641  } else if (copyno1[3] < 8) {
642  //UserNtuples->fillg577(dx,1.);
643  } else if (copyno1[3] < 9) {
644  //UserNtuples->fillg578(dx,1.);
645  } else if (copyno1[3] < 10) {
646  //UserNtuples->fillg579(dx,1.);
647  }
648  } else if (copyno1[2] < 4) {
649  //UserNtuples->fillg565(dx,1.);
650  } else if (copyno1[2] < 5) {
651  //UserNtuples->fillg567(dx,1.);
652  }
653  }
654  if (prename == "SIDETR") {
655  //UserNtuples->fillg542(EnerDeposit,1.);
656  //UserNtuples->fillg562(numbcont,1.);
657  if (copyno1[2] < 2) {
658  //UserNtuples->fillg581(dy,1.);
659  } else if (copyno1[2] < 3) {
660  //UserNtuples->fillg564(dy,1.);
661  if (copyno1[3] < 2) {
662  } else if (copyno1[3] < 3) {
663  //UserNtuples->fillg582(dy,1.);
664  } else if (copyno1[3] < 4) {
665  //UserNtuples->fillg583(dy,1.);
666  } else if (copyno1[3] < 5) {
667  //UserNtuples->fillg584(dy,1.);
668  } else if (copyno1[3] < 6) {
669  //UserNtuples->fillg585(dy,1.);
670  } else if (copyno1[3] < 7) {
671  //UserNtuples->fillg586(dy,1.);
672  } else if (copyno1[3] < 8) {
673  //UserNtuples->fillg587(dy,1.);
674  } else if (copyno1[3] < 9) {
675  //UserNtuples->fillg588(dy,1.);
676  } else if (copyno1[3] < 10) {
677  //UserNtuples->fillg589(dy,1.);
678  }
679  } else if (copyno1[2] < 4) {
680  //UserNtuples->fillg566(dy,1.);
681  } else if (copyno1[2] < 5) {
682  //UserNtuples->fillg568(dy,1.);
683  }
684  }
685  }
686  }
687  // end of prenames SIDETL // SIDETR
688  }
689  // end of trackstatus != 2
690 
691  // deposition of energy on steps along primary track
692  //UserNtuples->fillg500(EnerDeposit,1.);
693  // collect sum deposited energy on all steps along primary track
694  SumEnerDeposit += EnerDeposit;
695  // position of step for primary track:
696  //UserNtuples->fillg501(preposition.x(),1.);
697  //UserNtuples->fillg502(preposition.y(),1.);
698  //UserNtuples->fillg503(preposition.z(),1.);
699  //UserNtuples->fillg504(preposition.x(),EnerDeposit);
700  //UserNtuples->fillg505(preposition.y(),EnerDeposit);
701  //UserNtuples->fillg506(preposition.z(),EnerDeposit);
702  // 2D step coordinates weighted by energy deposited on step
703  // //UserNtuples->fillp201(preposition.x(),preposition.y(),EnerDeposit);
704  // //UserNtuples->fillp202(preposition.x(),preposition.z(),EnerDeposit);
705  // //UserNtuples->fillp203(preposition.y(),preposition.z(),EnerDeposit);
706  //UserNtuples->filld204(preposition.x(),preposition.y(),EnerDeposit);
707  //UserNtuples->filld205(preposition.x(),preposition.z(),EnerDeposit);
708  //UserNtuples->filld206(preposition.y(),preposition.z(),EnerDeposit);
709  //UserNtuples->filld223(preposition.y(),preposition.z(),EnerDeposit);
710  // last step of primary track
711  if (trackstatus == 2) {
712  // primary track length
713  // //UserNtuples->fillg508(tracklength,1.);
714  tracklength0 = tracklength;
715  // how many steps primary track consist
716  //UserNtuples->fillg509(curstepnumber,1.);
717  // tot. energy of primary track at the end of trajectory(before disappeare)
718  //UserNtuples->fillg510((entot0-entot),1.);
719  //UserNtuples->fillg520((entot0-entot),1.);
720  }
721  }
722  // end of primary track !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
723 
724  if (parentID == 1 && curstepnumber == 1) {
725  // particles deposit their energy along primary track
726  numofpart += 1;
727  // energy of radiated particle
728  //UserNtuples->fillg511(entot,1.);
729  // coordinates of radiated particle
730  //UserNtuples->fillg512(vert_pos.x(),1.);
731  //UserNtuples->fillg513(vert_pos.y(),1.);
732  //UserNtuples->fillg514(vert_pos.z(),1.);
733  //UserNtuples->fillg515(vert_pos.x(),entot);
734  //UserNtuples->fillg516(vert_pos.y(),entot);
735  //UserNtuples->fillg517(vert_pos.z(),entot);
736 
737  //UserNtuples->filld214(vert_pos.x(),vert_pos.y(),1.);
738  //UserNtuples->filld215(vert_pos.x(),vert_pos.z(),1.);
739  //UserNtuples->filld216(vert_pos.y(),vert_pos.z(),1.);
740  //UserNtuples->filld219(vert_pos.y(),vert_pos.z(),1.);
741  //UserNtuples->filld224(vert_pos.y(),vert_pos.z(),1.);
742  if (prename == "SBST") {
743  //UserNtuples->filld225(vert_pos.y(),vert_pos.z(),1.);
744  }
745  if (prename == "SBSTs") {
746  //UserNtuples->filld226(vert_pos.y(),vert_pos.z(),1.);
747  }
748  }
749 
750  // ==========================================================================
751 }
752 // ==========================================================================
753 // ==========================================================================
754 int FP420Test::detLevels(const G4VTouchable* touch) const {
755  //Return number of levels
756  if (touch)
757  return ((touch->GetHistoryDepth()) + 1);
758  else
759  return 0;
760 }
761 // ==========================================================================
762 
763 G4String FP420Test::detName(const G4VTouchable* touch, int level, int currentlevel) const {
764  //Go down to current level
765  if (level > 0 && level >= currentlevel) {
766  int ii = level - currentlevel;
767  return touch->GetVolume(ii)->GetName();
768  } else {
769  return "NotFound";
770  }
771 }
772 
773 void FP420Test::detectorLevel(const G4VTouchable* touch, int& level, int* copyno, G4String* name) const {
774  //Get name and copy numbers
775  if (level > 0) {
776  for (int ii = 0; ii < level; ii++) {
777  int i = level - ii - 1;
778  G4VPhysicalVolume* pv = touch->GetVolume(i);
779  if (pv != nullptr)
780  name[ii] = pv->GetName();
781  else
782  name[ii] = "Unknown";
783  copyno[ii] = touch->GetReplicaNumber(i);
784  }
785  }
786 }
787 // ==========================================================================
788 
789 //=================================================================== End Of Event
790 void FP420Test::update(const EndOfEvent* evt) {
791  // ==========================================================================
792 
793  if (verbosity > 1) {
794  iev = (*evt)()->GetEventID();
795  std::cout << "FP420Test:update EndOfEvent = " << iev << std::endl;
796  }
797  // Fill-in ntuple
798  fp420eventarray[ntfp420_evt] = (float)whichevent;
799 
800  //
801  int trackID = 0;
802  G4PrimaryParticle* thePrim = nullptr;
803 
804  // prim.vertex:
805  G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
806  if (nvertex != 1)
807  std::cout << "FP420Test: My warning: NumberOfPrimaryVertex != 1 --> = " << nvertex << std::endl;
808 
809  for (int i = 0; i < nvertex; i++) {
810  G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
811  if (avertex == nullptr)
812  std::cout << "FP420Test End Of Event ERR: pointer to vertex = 0" << std::endl;
813  G4int npart = avertex->GetNumberOfParticle();
814  if (npart != 1)
815  std::cout << "FP420Test: My warning: NumberOfPrimaryPart != 1 --> = " << npart << std::endl;
816  if (npart == 0)
817  std::cout << "FP420Test End Of Event ERR: no NumberOfParticle" << std::endl;
818 
819  // find just primary track: track pointer: thePrim
820  if (thePrim == nullptr)
821  thePrim = avertex->GetPrimary(trackID);
822 
823  if (thePrim != nullptr) {
824  // primary vertex:
825  G4double vx = 0., vy = 0., vz = 0.;
826  vx = avertex->GetX0();
827  vy = avertex->GetY0();
828  vz = avertex->GetZ0();
829  //UserNtuples->fillh01(vx);
830  //UserNtuples->fillh02(vy);
831  //UserNtuples->fillh03(vz);
832  TheHistManager->GetHisto("VtxX")->Fill(vx);
833  TheHistManager->GetHisto("VtxY")->Fill(vy);
834  TheHistManager->GetHisto("VtxZ")->Fill(vz);
835  }
836  }
837  // prim.vertex loop end
838 
839  //=========================== thePrim != 0 ================================================================================
840  if (thePrim != nullptr) {
841  // inline G4ParticleDefinition * GetG4code() const
842  // inline G4PrimaryParticle * GetNext() const
843  // inline G4PrimaryParticle * GetDaughter() const
844 
845  // prim.vertex
846  //int ivert = 0 ;
847  //G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(ivert);
848  //G4double vx=0.,vy=0.,vz=0.;
849  //vx = avertex->GetX0();
850  //vy = avertex->GetY0();
851  //vz = avertex->GetZ0();
852 
853  //
854  // number of secondary particles deposited their energy along primary track
855  //UserNtuples->fillg518(numofpart,1.);
856  if (lastpo.z() < z4 && lastpo.perp() < 100.) {
857  //UserNtuples->fillg536(numofpart,1.);
858  }
859  //
860 
861  // direction !!!
862  G4ThreeVector mom = thePrim->GetMomentum();
863 
864  double phi = atan2(mom.y(), mom.x());
865  if (phi < 0.)
866  phi += twopi;
867  double phigrad = phi * 180. / pi;
868 
869  double th = mom.theta();
870  double eta = -log(tan(th / 2));
871  // works OK:
872  // double costheta =mom.z()/sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z());
873  // double theta = acos(min(max(costheta,double(-1.)),double(1.)));
874 
875  //UserNtuples->fillh04(eta);
876  //UserNtuples->fillh05(phigrad);
877  //UserNtuples->fillh06(th);
878 
879  //UserNtuples->fillg13(lastpo.x(),1.);
880  //UserNtuples->fillg14(lastpo.y(),1.);
881  //UserNtuples->fillg15(lastpo.z(),1.);
882 
883  TheHistManager->GetHisto("PrimaryEta")->Fill(eta);
884  TheHistManager->GetHisto("PrimaryPhigrad")->Fill(phigrad);
885  TheHistManager->GetHisto("PrimaryTh")->Fill(th * 180. / pi);
886 
887  TheHistManager->GetHisto("PrimaryLastpoZ")->Fill(lastpo.z());
888  if (lastpo.z() < z4) {
889  TheHistManager->GetHisto("PrimaryLastpoX")->Fill(lastpo.x());
890  TheHistManager->GetHisto("PrimaryLastpoY")->Fill(lastpo.y());
891  }
892  if (numofpart > 4) {
893  TheHistManager->GetHisto("XLastpoNumofpart")->Fill(lastpo.x());
894  TheHistManager->GetHisto("YLastpoNumofpart")->Fill(lastpo.y());
895  }
896 
897  // ==========================================================================
898 
899  // hit map for FP420
900  // ==================================
901 
902  map<int, float, less<int> > themap;
903  map<int, float, less<int> > themap1;
904 
905  map<int, float, less<int> > themapxy;
906  map<int, float, less<int> > themapz;
907  // access to the G4 hit collections: -----> this work OK:
908 
909  // edm::LogInfo("FP420Test") << "1";
910  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
911  // edm::LogInfo("FP420Test") << "2";
912  if (verbosity > 0) {
913  std::cout << "FP420Test: accessed all HC" << std::endl;
914  ;
915  }
916  int CAFIid = G4SDManager::GetSDMpointer()->GetCollectionID("FP420SI");
917  // edm::LogInfo("FP420Test") << "3";
918  // std::cout << " CAFIid = " << CAFIid << std::endl;;
919 
920  FP420G4HitCollection* theCAFI = (FP420G4HitCollection*)allHC->GetHC(CAFIid);
921  // CaloG4HitCollection* theCAFI = (CaloG4HitCollection*) allHC->GetHC(CAFIid);
922  if (verbosity > 0) {
923  //std::cout << "FP420Test: theCAFI = " << theCAFI << std::endl;
924  std::cout << "FP420Test: theCAFI->entries = " << theCAFI->entries() << std::endl;
925  }
926  // edm::LogInfo("FP420Test") << "theCAFI->entries="<< theCAFI->entries();
927  TheHistManager->GetHisto("NumberOfHits")->Fill(theCAFI->entries());
928 
929  // access to the G4 hit collections ----> variant 2: give 0 hits
930  // FP420G4HitCollection * theCAFI;
931  // theCAFI = new FP420G4HitCollection();
932  // ==========================================================================
933  // Silicon Hit collection start
934  //0) if particle goes into flat beam pipe below detector:
935  int varia; // = 0 -all; =1 - MI; =2 - noMI
936  // Select MI or noMI over all 3 stations
937  // 1)MI:
938  // if particle goes through window into detector:
939  // lastpoint of track in lateral dir. outside the detector and inside in z
940  // lastpoint of track in lateral dir. outside the detector and inside in z
941  // for all except zzzmarta.xml
942  // if( lastpo.z()<z4 || abs(lastpo.x())> 5. || lastpo.y()< 10.2 || lastpo.y()>30.2 ) {
943  // for zzzmarta.xml
944  // if( lastpo.z()<z4 || abs(lastpo.x())> 10. || lastpo.y()< 3.2 || lastpo.y()>43.2 ) {
945  if (lastpo.z() < z4) {
946  // if( lastpo.z()<z4 && lastpo.perp()< 100. ) {
947  // if(lastpo.z()<z4 || lastpo.perp()> 45.) {
948  //UserNtuples->fillg66(theCAFI->entries(),1.);
949  varia = 1;
950  } else {
951  // 2) no MI start in detector:
952  //UserNtuples->fillg31(numofpart,1.);
953  //UserNtuples->fillg67(theCAFI->entries(),1.);
954  varia = 2;
955  } // no MI end:
956  int nhits = theCAFI->entries();
957  for (int j = 0; j < nhits; j++) {
958  FP420G4Hit* aHit = (*theCAFI)[j];
959  G4ThreeVector hitPoint = aHit->getEntry();
960  double zz = hitPoint.z();
961  TheHistManager->GetHisto("zHits")->Fill(zz);
962  if (tracklength0 > 8300.)
963  TheHistManager->GetHisto("zHitsTrLoLe")->Fill(zz);
964  }
965  // varia = 0;
966  // if( varia == 0) {
967  if (varia == 2) {
968  // .............
969  // number of hits < 50
970  // if(theCAFI->entries() <50) {
971  // if(theCAFI->entries() > 0) {
972  // if(theCAFI->entries() > -1) {
973  // .............
974  int nhit11 = 0, nhit12 = 0, nhit13 = 0;
975  double totallosenergy = 0.;
976  for (int j = 0; j < nhits; j++) {
977  FP420G4Hit* aHit = (*theCAFI)[j];
978 
979  G4ThreeVector hitEntryLocalPoint = aHit->getEntryLocalP();
980  G4ThreeVector hitExitLocalPoint = aHit->getExitLocalP();
981  G4ThreeVector hitPoint = aHit->getEntry();
982  // double elmenergy = aHit->getEM();
983  // double hadrenergy = aHit->getHadr();
984  // double incidentEnergyHit = aHit->getIncidentEnergy();
985  int trackIDhit = aHit->getTrackID();
986  unsigned int unitID = aHit->getUnitID();
987  // double timeslice = aHit->getTimeSlice();
988  // int timesliceID = aHit->getTimeSliceID();
989  // double depenergy = aHit->getEnergyDeposit();
990  // float pabs = aHit->getPabs();
991  // float tof = aHit->getTof();
992  double losenergy = aHit->getEnergyLoss();
993  // int particletype = aHit->getParticleType();
994  // float thetaEntry = aHit->getThetaAtEntry();
995  // float phiEntry = aHit->getPhiAtEntry();
996  // float xEntry = aHit->getX();
997  // float yEntry = aHit->getY();
998  // float zEntry = aHit->getZ();
999  // int parentID = aHit->getParentId();
1000  // float vxget = aHit->getVx();
1001  // float vyget = aHit->getVy();
1002  // float vzget = aHit->getVz();
1003 
1004  // double th_hit = hitPoint.theta();
1005  // double eta_hit = -log(tan(th_hit/2));
1006  // double phi_hit = hitPoint.phi();
1007  // if (phi_hit < 0.) phi_hit += twopi;
1008  // double phigrad_hit = phi_hit*180./pi;
1009  //UserNtuples->fillg60(eta_hit,losenergy);
1010  //UserNtuples->fillg61(eta_hit,1.);
1011  //UserNtuples->fillg62(phigrad_hit,losenergy);
1012  //UserNtuples->fillg63(phigrad_hit,1.);
1013 
1014  // double xx = hitPoint.x();
1015  // double yy = hitPoint.y();
1016  double zz = hitPoint.z();
1017 
1018  TheHistManager->GetHisto("zHitsnoMI")->Fill(zz);
1019 
1020  if (verbosity > 2) {
1021  std::cout << "FP420Test:zHits = " << zz << std::endl;
1022  }
1023  // double rr = hitPoint.perp();
1024  /*
1025  if(aHit->getTrackID() == 1) {
1026  emu += aHit->getEnergyDeposit();} else {
1027  erest += aHit->getEnergyDeposit();}
1028  */
1029 
1030  //collect lost in Si en.of hits in every plane and put it into themap[]
1031  //UserNtuples->fillg30(losenergy,1.);
1032  themap[unitID] += losenergy;
1033  totallosenergy += losenergy;
1034 
1035  int det, zside, sector, zmodule;
1036  // CaloNumberingPacker::unpackCastorIndex(unitID, det, zside, sector, zmodule);
1037  FP420NumberingScheme::unpackFP420Index(unitID, det, zside, sector, zmodule);
1038  int justlayer = FP420NumberingScheme::unpackLayerIndex(rn00, zside); // 1,2
1039  if (justlayer < 1 || justlayer > 2) {
1040  std::cout << "FP420Test:WRONG justlayer= " << justlayer << std::endl;
1041  }
1042  // zside=1,2 ; zmodule=1,10 ; sector=1,3
1043  //UserNtuples->fillg44(float(sector),1.);
1044  //UserNtuples->fillg45(float(zmodule),1.);
1045  //UserNtuples->fillg46(float(zside),1.);
1046  // int sScale = 20;
1047  // intindex is a continues numbering of FP420
1048  //int zScale = 2; unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside; //intindex=1-30:X,Y,X,Y,X,Y...
1049  // int zScale = 10; unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule; //intindex=1-30:XXXXXXXXXX,YYYYYYYYYY,...
1050  //UserNtuples->fillg40(float(intindex),1.);
1051  //UserNtuples->fillg48(float(intindex),losenergy);
1052  //
1053  //=======================================
1054  // G4ThreeVector middle = (hitExitLocalPoint+hitEntryLocalPoint)/2.;
1055  G4ThreeVector middle = (hitExitLocalPoint - hitEntryLocalPoint) / 2.;
1056  themapz[unitID] = hitPoint.z() + fabs(middle.z());
1057  if (verbosity > 2) {
1058  std::cout << "1111111111111111111111111111111111111111111111111111111111111111111111111 " << std::endl;
1059  std::cout << "FP420Test: det, zside, sector, zmodule = " << det << zside << sector << zmodule << std::endl;
1060  std::cout << "FP420Test: justlayer = " << justlayer << std::endl;
1061  std::cout << "FP420Test: hitExitLocalPoint = " << hitExitLocalPoint << std::endl;
1062  std::cout << "FP420Test: hitEntryLocalPoint = " << hitEntryLocalPoint << std::endl;
1063  std::cout << "FP420Test: middle= " << middle << std::endl;
1064  std::cout << "FP420Test: hitPoint.z()-419000.= " << hitPoint.z() - 419000. << std::endl;
1065 
1066  std::cout << "FP420Test:zHits-419000. = " << themapz[unitID] - 419000. << std::endl;
1067  }
1068  //=======================================
1069  // Y
1070  if (zside == 1) {
1071  //UserNtuples->fillg24(losenergy,1.);
1072  if (losenergy > 0.00003) {
1073  themap1[unitID] += 1.;
1074  }
1075  }
1076  //X
1077  if (zside == 2) {
1078  //UserNtuples->fillg25(losenergy,1.);
1079  if (losenergy > 0.00005) {
1080  themap1[unitID] += 1.;
1081  }
1082  }
1083  // }
1084  //
1085  if (sector == 1) {
1086  nhit11 += 1;
1087  //UserNtuples->fillg33(rr,1.);
1088  //UserNtuples->fillg11(yy,1.);
1089  }
1090  if (sector == 2) {
1091  nhit12 += 1;
1092  //UserNtuples->fillg34(rr,1.);
1093  //UserNtuples->fillg86(yy,1.);
1094  }
1095  if (sector == 3) {
1096  nhit13 += 1;
1097  //UserNtuples->fillg35(rr,1.);
1098  //UserNtuples->fillg87(yy,1.);
1099  }
1100  //UserNtuples->fillg10(xx,1.);
1101  //UserNtuples->fillg12(zz,1.);
1102  //UserNtuples->fillg32(rr,1.);
1103 
1104  // =========
1105  // double xPrimAtZhit = vx + (zz-vz)*tan(th)*cos(phi);
1106  // double yPrimAtZhit = vy + (zz-vz)*tan(th)*sin(phi);
1107 
1108  // double dx = xPrimAtZhit - xx;
1109  // double dy = yPrimAtZhit - yy;
1110 
1111  // Select SD hits
1112  // if(rr<120.) {
1113  // Select MI or noMI over all 3 stations
1114  if (lastpo.z() < z4 && lastpo.perp() < 120.) {
1115  // MIonly:
1116  //UserNtuples->fillg16(lastpo.z(),1.);
1117  //UserNtuples->fillg18(zz,1.);
1118  // Station I
1119  if (zz < z2) {
1120  //UserNtuples->fillg54(dx,1.);
1121  //UserNtuples->fillg55(dy,1.);
1122  }
1123  // Station II
1124  if (zz < z3 && zz > z2) {
1125  //UserNtuples->fillg50(dx,1.);
1126  //UserNtuples->fillg51(dy,1.);
1127  }
1128  // Station III
1129  if (zz < z4 && zz > z3) {
1130  //UserNtuples->fillg64(dx,1.);
1131  //UserNtuples->fillg65(dy,1.);
1132  //UserNtuples->filld209(xx,yy,1.);
1133  }
1134  } else {
1135  // no MIonly:
1136  //UserNtuples->fillg17(lastpo.z(),1.);
1137  //UserNtuples->fillg19(zz,1.);
1138  //UserNtuples->fillg74(incidentEnergyHit,1.);
1139  //UserNtuples->fillg75(float(trackIDhit),1.);
1140  // Station I
1141  if (zz < z2) {
1142  //UserNtuples->fillg56(dx,1.);
1143  //UserNtuples->fillg57(dy,1.);
1144  //UserNtuples->fillg20(numofpart,1.);
1145  //UserNtuples->fillg21(SumEnerDeposit,1.);
1146  if (zside == 1) {
1147  //UserNtuples->fillg26(losenergy,1.);
1148  }
1149  if (zside == 2) {
1150  //UserNtuples->fillg76(losenergy,1.);
1151  }
1152  if (trackIDhit == 1) {
1153  //UserNtuples->fillg70(dx,1.);
1154  //UserNtuples->fillg71(incidentEnergyHit,1.);
1155  //UserNtuples->fillg79(losenergy,1.);
1156  } else {
1157  //UserNtuples->fillg82(dx,1.);
1158  }
1159  }
1160  // Station II
1161  if (zz < z3 && zz > z2) {
1162  //UserNtuples->fillg52(dx,1.);
1163  //UserNtuples->fillg53(dy,1.);
1164  //UserNtuples->fillg22(numofpart,1.);
1165  //UserNtuples->fillg23(SumEnerDeposit,1.);
1166  //UserNtuples->fillg80(incidentEnergyHit,1.);
1167  //UserNtuples->fillg81(float(trackIDhit),1.);
1168  if (zside == 1) {
1169  //UserNtuples->fillg27(losenergy,1.);
1170  }
1171  if (zside == 2) {
1172  //UserNtuples->fillg77(losenergy,1.);
1173  }
1174  if (trackIDhit == 1) {
1175  //UserNtuples->fillg72(dx,1.);
1176  //UserNtuples->fillg73(incidentEnergyHit,1.);
1177  } else {
1178  //UserNtuples->fillg83(dx,1.);
1179  }
1180  }
1181  // Station III
1182  if (zz < z4 && zz > z3) {
1183  //UserNtuples->fillg68(dx,1.);
1184  //UserNtuples->fillg69(dy,1.);
1185  //UserNtuples->filld210(xx,yy,1.);
1186  //UserNtuples->fillg22(numofpart,1.);
1187  //UserNtuples->fillg23(SumEnerDeposit,1.);
1188  if (zside == 1) {
1189  //UserNtuples->fillg28(losenergy,1.);
1190  }
1191  if (zside == 2) {
1192  //UserNtuples->fillg78(losenergy,1.);
1193  }
1194  }
1195  } // MIonly or noMIonly ENDED
1196  // }
1197 
1198  // !!!!!!!!!!!!!
1199 
1200  } // for loop on all hits ENDED ENDED ENDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1201 
1202  // !!!!!!!!!!!!!
1203 
1204  //======================================================================================================SUMHIT
1205  //UserNtuples->fillg29(totallosenergy,1.);
1206  //UserNtuples->fillg36(nhit11,1.);
1207  //UserNtuples->fillg37(nhit12,1.);
1208  //UserNtuples->fillg38(nhit13,1.);
1209  //======================================================================================================SUMHIT
1210  // int rn00=3;//test only with 2 sensors in superlayer, not 4
1211  // int rn00=rn0;//always
1212  if (verbosity > 2) {
1213  std::cout << "22222222222222222222222222222222222222222222222222222222222222222222222222 " << std::endl;
1214  }
1215  int det = 1;
1216  int allplacesforsensors = 7;
1217  for (int sector = 1; sector < sn0; sector++) {
1218  for (int zmodule = 1; zmodule < pn0; zmodule++) {
1219  for (int zsideinorder = 1; zsideinorder < allplacesforsensors; zsideinorder++) {
1220  int zside = FP420NumberingScheme::realzside(rn00, zsideinorder); //1,3,5,2,4,6
1221  if (verbosity > 2) {
1222  std::cout << "FP420Test: sector= " << sector << " zmodule= " << zmodule
1223  << " zsideinorder= " << zsideinorder << " zside= " << zside << std::endl;
1224  }
1225  if (zside != 0) {
1226  int justlayer = FP420NumberingScheme::unpackLayerIndex(rn00, zside); // 1,2
1227  if (justlayer < 1 || justlayer > 2) {
1228  std::cout << "FP420Test:WRONG justlayer= " << justlayer << std::endl;
1229  }
1230  int copyinlayer = FP420NumberingScheme::unpackCopyIndex(rn00, zside); // 1,2,3
1231  if (copyinlayer < 1 || copyinlayer > 3) {
1232  std::cout << "FP420Test:WRONG copyinlayer= " << copyinlayer << std::endl;
1233  }
1234  int orientation = FP420NumberingScheme::unpackOrientation(rn00, zside); // Front: = 1; Back: = 2
1235  if (orientation < 1 || orientation > 2) {
1236  std::cout << "FP420Test:WRONG orientation= " << orientation << std::endl;
1237  }
1238 
1239  // iu is a continues numbering of planes(!) over two arm FP420 set up
1240  int detfixed =
1241  1; // use this treatment for each set up arm, hence no sense to do it defferently for +FP420 and -FP420;
1242  // and ...[ii] massives have prepared in such a way
1243  unsigned int ii =
1244  FP420NumberingScheme::packMYIndex(rn00, pn0, sn0, detfixed, justlayer, sector, zmodule) - 1;
1245  // ii = 0-19 --> 20 items
1246  if (verbosity > 2) {
1247  std::cout << "FP420Test: justlayer = " << justlayer << " copyinlayer = " << copyinlayer
1248  << " orientation = " << orientation << " ii= " << ii << std::endl;
1249  }
1250  double zdiststat = 0.;
1251  if (sn0 < 4) {
1252  if (sector == 2)
1253  zdiststat = zD3;
1254  } else {
1255  if (sector == 2)
1256  zdiststat = zD2;
1257  if (sector == 3)
1258  zdiststat = zD3;
1259  }
1260  double kplane = -(pn0 - 1) / 2 - 0.5 + (zmodule - 1); //-3.5 +0...5 = -3.5,-2.5,-1.5,+2.5,+1.5
1261  double zcurrent = zinibeg + z420 + (ZSiStep - ZSiPlane) / 2 + kplane * ZSiStep + zdiststat;
1262  //double zcurrent = zinibeg +(ZSiStep-ZSiPlane)/2 + kplane*ZSiStep + (sector-1)*zUnit;
1263  if (verbosity > 2) {
1264  std::cout << "FP420Test: Leftzcurrent-419000. = " << zcurrent - 419000. << std::endl;
1265  std::cout << "FP420Test: ZGapLDet = " << ZGapLDet << std::endl;
1266  }
1267  if (justlayer == 1) {
1268  if (orientation == 1)
1269  zcurrent += (ZGapLDet + ZSiDet / 2);
1270  if (orientation == 2)
1271  zcurrent += zBlade - (ZGapLDet + ZSiDet / 2);
1272  }
1273  if (justlayer == 2) {
1274  if (orientation == 1)
1275  zcurrent += (ZGapLDet + ZSiDet / 2) + zBlade + gapBlade;
1276  if (orientation == 2)
1277  zcurrent += 2 * zBlade + gapBlade - (ZGapLDet + ZSiDet / 2);
1278  }
1279  // .
1280  //
1281  if (det == 2)
1282  zcurrent = -zcurrent;
1283  //
1284  if (verbosity > 2) {
1285  std::cout << "FP420Test: zcurrent-419000. = " << zcurrent - 419000. << std::endl;
1286  }
1287  //================================== end of for loops in continuius number iu:
1288  } //if(zside!=0
1289  } // for superlayer
1290  } // for zmodule
1291  } // for sector
1292 
1293  if (verbosity > 2) {
1294  std::cout << "----------------------------------------------------------------------------- " << std::endl;
1295  }
1296 
1297  //======================================================================================================CHECK
1298  if (totallosenergy == 0.0) {
1299  std::cout << "FP420Test: number of hits = " << theCAFI->entries() << std::endl;
1300  for (int j = 0; j < nhits; j++) {
1301  FP420G4Hit* aHit = (*theCAFI)[j];
1302  double losenergy = aHit->getEnergyLoss();
1303  std::cout << " j hits = " << j << "losenergy = " << losenergy << std::endl;
1304  }
1305  }
1306  //======================================================================================================CHECK
1307 
1308  //====================================================================================================== HIT START
1309 
1310  // FIBRE Hit collected analysis
1311  double totalEnergy = 0.;
1312  int nhitsX = 0, nhitsY = 0, nsumhit = 0;
1313  for (int sector = 1; sector < 4; sector++) {
1314  int nhitsecX = 0, nhitsecY = 0;
1315  for (int zmodule = 1; zmodule < 11; zmodule++) {
1316  for (int zside = 1; zside < 3; zside++) {
1317  int det = 1;
1318  // int nhit = 0;
1319  // int sScale = 20;
1320  int index = FP420NumberingScheme::packFP420Index(det, zside, sector, zmodule);
1321  double theTotalEnergy = themap[index];
1322  // X planes
1323  if (zside < 2) {
1324  //UserNtuples->fillg47(theTotalEnergy,1.);
1325  if (theTotalEnergy > 0.00003) {
1326  nhitsX += 1;
1327  // nhitsecX += themap1[index];
1328  // nhit=1;
1329  }
1330  }
1331  // Y planes
1332  else {
1333  //UserNtuples->fillg49(theTotalEnergy,1.);
1334  if (theTotalEnergy > 0.00005) {
1335  nhitsY += 1;
1336  // nhitsecY += themap1[index];
1337  // nhit=1;
1338  }
1339  }
1340  // intindex is a continues numbering of FP420
1341  // int zScale=2; unsigned int intindex = sScale*(sector - 1)+zScale*(zmodule - 1)+zside;
1342  // int zScale=10; unsigned int intindex = sScale*(sector - 1)+zScale*(zside - 1)+zmodule;
1343  //UserNtuples->fillg41(float(intindex),theTotalEnergy);
1344  //UserNtuples->fillg42(float(intindex),1.);
1345  //UserNtuples->fillp208(float(intindex),float(nhit),1.);
1346  //UserNtuples->fillp211(float(intindex),float(themap1[index]),1.);
1347  totalEnergy += themap[index];
1348  } // for
1349  } // for
1350  //UserNtuples->fillg39(nhitsecY,1.);
1351  if (nhitsecX > 10 || nhitsecY > 10) {
1352  nsumhit += 1;
1353  //UserNtuples->fillp213(float(sector),float(1.),1.);
1354  } else { //UserNtuples->fillp213(float(sector),float(0.),1.);
1355  }
1356  } // for
1357  //====================================================================================================== HIT END
1358 
1359  //====================================================================================================== HIT ALL
1360  //UserNtuples->fillg43(totalEnergy,1.);
1361  //UserNtuples->fillg58(nhitsX,1.);
1362  //UserNtuples->fillg59(nhitsY,1.);
1363  // if( nsumhit !=0 ) { //UserNtuples->fillp212(vy,float(1.),1.);
1364  if (nsumhit >= 2) { //UserNtuples->fillp212(vy,float(1.),1.);
1365  } else { //UserNtuples->fillp212(vy,float(0.),1.);
1366  }
1367 
1368  //====================================================================================================== HIT ALL
1369 
1370  //====================================================================================================== number of hits
1371  // .............
1372  // } // number of hits < 50
1373  // .............
1374  } // MI or no MI or all - end
1375 
1376  } // primary end
1377  //=========================== thePrim != 0 end ===
1378  // ==========================================================================
1379  if (verbosity > 0) {
1380  std::cout << "FP420Test: END OF Event " << (*evt)()->GetEventID() << std::endl;
1381  }
1382 }
1383 
1384 // ==========================================================================
personalPlayback.level
level
Definition: personalPlayback.py:22
HIPAlignmentAlgorithm_cfi.verbosity
verbosity
Definition: HIPAlignmentAlgorithm_cfi.py:7
iev
const HitContainer *__restrict__ const TkSoA *__restrict__ const Quality *__restrict__ const CAHitNtupletGeneratorKernelsGPU::HitToTuple *__restrict__ uint32_t int iev
Definition: CAHitNtupletGeneratorKernelsImpl.h:544
FP420G4Hit::getEnergyLoss
float getEnergyLoss() const
Definition: FP420G4Hit.cc:140
Fp420AnalysisHistManager::BookHistos
void BookHistos()
Definition: FP420Test.cc:176
runGCPTkAlMap.title
string title
Definition: runGCPTkAlMap.py:94
mps_fire.i
i
Definition: mps_fire.py:428
geometryCSVtoXML.zz
zz
Definition: geometryCSVtoXML.py:19
MessageLogger.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
FP420Track_cfi.z420
z420
Definition: FP420Track_cfi.py:10
FP420G4HitCollection.h
FP420Track_cfi.zD2
zD2
Definition: FP420Track_cfi.py:12
ntfp420_evt
Definition: FP420Test.cc:49
ecaldqm::zside
int zside(DetId const &)
Definition: EcalDQMCommonUtils.cc:189
FP420Test::~FP420Test
~FP420Test() override
Definition: FP420Test.cc:119
edm
HLT enums.
Definition: AlignableModifier.h:19
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
FP420NumberingScheme::unpackOrientation
static int unpackOrientation(int rn0, int zside)
Definition: FP420NumberingScheme.h:88
FP420G4HitCollection
G4THitsCollection< FP420G4Hit > FP420G4HitCollection
Definition: FP420G4HitCollection.h:14
gather_cfg.cout
cout
Definition: gather_cfg.py:144
FP420Test.h
FP420NumberingScheme.h
FP420Test::detectorLevel
void detectorLevel(const G4VTouchable *, int &, int *, G4String *) const
Definition: FP420Test.cc:773
HcalTBWriter_cfi.Number
Number
Definition: HcalTBWriter_cfi.py:19
EndOfEvent.h
Fp420AnalysisHistManager::GetHisto2
TH2F * GetHisto2(Int_t Number)
Definition: FP420Test.cc:257
npart
double npart
Definition: HydjetWrapper.h:46
EndOfTrack
Definition: EndOfTrack.h:6
FP420NumberingScheme::unpackFP420Index
static void unpackFP420Index(const unsigned int &idx, int &det, int &zside, int &station, int &superplane)
Definition: FP420NumberingScheme.cc:163
testProducerWithPsetDescEmpty_cfi.z2
z2
Definition: testProducerWithPsetDescEmpty_cfi.py:41
FP420Track_cfi.zD3
zD3
Definition: FP420Track_cfi.py:11
FP420NumberingScheme::packFP420Index
static unsigned int packFP420Index(int det, int zside, int station, int superplane)
Definition: FP420NumberingScheme.cc:143
FP420G4Hit::getUnitID
unsigned int getUnitID() const
Definition: FP420G4Hit.cc:124
FP420NumberingScheme::unpackCopyIndex
static int unpackCopyIndex(int rn0, int zside)
Definition: FP420NumberingScheme.h:69
FP420G4Hit::getEntry
G4ThreeVector getEntry() const
Definition: FP420G4Hit.cc:103
PVValHelper::eta
Definition: PVValidationHelpers.h:70
BeginOfTrack
Definition: BeginOfTrack.h:6
BeginOfJob
Definition: BeginOfJob.h:8
Fp420AnalysisHistManager::WriteToFile
void WriteToFile(const TString &fOutputFile, const TString &fRecreateFile)
Definition: FP420Test.cc:200
nhits
Definition: HIMultiTrackSelector.h:42
EndOfEvent
Definition: EndOfEvent.h:6
FP420NumberingScheme::packMYIndex
static unsigned packMYIndex(int rn0, int pn0, int sn0, int det, int zside, int sector, int zmodule)
Definition: FP420NumberingScheme.h:37
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
Fp420AnalysisHistManager::HistInit
void HistInit(const char *name, const char *title, Int_t nbinsx, Axis_t xlow, Axis_t xup)
Definition: FP420Test.cc:214
Fp420AnalysisHistManager::~Fp420AnalysisHistManager
~Fp420AnalysisHistManager() override
Definition: FP420Test.cc:161
edm::ParameterSet
Definition: ParameterSet.h:47
Event.h
FP420NumberingScheme::realzside
static int realzside(int rn0, int zsideinorder)
Definition: FP420NumberingScheme.h:99
FP420NumberingScheme::unpackLayerIndex
static int unpackLayerIndex(int rn0, int zside)
Definition: FP420NumberingScheme.h:56
FP420G4Hit::getEntryLocalP
G4ThreeVector getEntryLocalP() const
Definition: FP420G4Hit.cc:106
BeginOfEvent.h
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
TrackInformation
Definition: TrackInformation.h:12
Fp420AnalysisHistManager::StoreWeights
void StoreWeights()
Definition: FP420Test.cc:289
FrontierConditions_GlobalTag_cff.file
file
Definition: FrontierConditions_GlobalTag_cff.py:13
MetAnalyzer.pv
def pv(vc)
Definition: MetAnalyzer.py:7
FP420G4Hit::getExitLocalP
G4ThreeVector getExitLocalP() const
Definition: FP420G4Hit.cc:109
BeginOfEvent
Definition: BeginOfEvent.h:6
BeginOfRun
Definition: BeginOfRun.h:6
EndOfRun
Definition: EndOfRun.h:6
Fp420AnalysisHistManager::GetHisto
TH1F * GetHisto(Int_t Number)
Definition: FP420Test.cc:241
TrackInformation.h
Fp420AnalysisHistManager::Fp420AnalysisHistManager
Fp420AnalysisHistManager(const TString &managername)
Definition: FP420Test.cc:144
Fp420AnalysisHistManager
Definition: FP420Test.h:58
FP420G4Hit
Definition: FP420G4Hit.h:18
std
Definition: JetResolutionObject.h:76
writedatasetfile.run
run
Definition: writedatasetfile.py:27
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
pi
const Double_t pi
Definition: trackSplitPlot.h:36
FP420Test::FP420Test
FP420Test(const edm::ParameterSet &p)
Definition: FP420Test.cc:51
PbPb_ZMuSkimMuonDPG_cff.particleType
particleType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:27
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
FP420Test::detName
G4String detName(const G4VTouchable *, int, int) const
Definition: FP420Test.cc:763
FP420G4Hit::getTrackID
unsigned int getTrackID() const
Definition: FP420G4Hit.cc:121
FP420Test::update
void update(const BeginOfJob *run) override
This routine will be called when the appropriate signal arrives.
Definition: FP420Test.cc:304
cuy.ii
ii
Definition: cuy.py:590
ntfp420_elements
ntfp420_elements
Definition: FP420Test.cc:49
TrackWithHistory.h
FP420Test::detLevels
int detLevels(const G4VTouchable *) const
Definition: FP420Test.cc:754
Generated for CMSSW Reference Manual by doxygen 1.8.16