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