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