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