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