B1EventAction.cc 31 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
/// \file B1EventAction.cc
/// \brief Implementation of the B1EventAction class

#include "B1EventAction.hh"
#include "B1RunAction.hh"
32
#include "B1Hits.hh"
33
#include "B1DetectorConstruction.hh"
34 35 36

#include "G4Event.hh"
#include "G4RunManager.hh"
37 38 39 40 41 42
#include "G4EventManager.hh"
#include "G4HCofThisEvent.hh"
#include "G4VHitsCollection.hh"
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"
43
#include "B1Analysis.hh"
44

45 46
#include <iostream>
#include <fstream>
47 48 49

using namespace std;

50 51 52 53 54 55
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B1EventAction::B1EventAction(B1RunAction* runAction)
: G4UserEventAction(),
  fRunAction(runAction),
  fEdep(0.)
56
{}
57 58 59 60 61 62 63 64

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B1EventAction::~B1EventAction()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

65 66
void B1EventAction::BeginOfEventAction(const G4Event* event)
{
67
  fEdep = 0.;
68

69 70


71 72 73 74
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

75 76
void B1EventAction::EndOfEventAction(const G4Event* event)
{
77 78
  // accumulate statistics in run action
  fRunAction->AddEdep(fEdep);
79

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
  G4SDManager * SDman = G4SDManager::GetSDMpointer();
  G4String detName;

  auto Col1 = SDman->GetCollectionID("SD1");
  auto Col2 = SDman->GetCollectionID("SD2");
  auto Col3 = SDman->GetCollectionID("SD3");
  auto Col4 = SDman->GetCollectionID("SD4");
  auto Col5 = SDman->GetCollectionID("SD5");
  auto Col6 = SDman->GetCollectionID("SD6");
  auto Col7 = SDman->GetCollectionID("SD7");
  auto Col8 = SDman->GetCollectionID("SD8");
  auto Col9 = SDman->GetCollectionID("SD9");
  auto Col10 = SDman->GetCollectionID("SD10");
  auto Col11 = SDman->GetCollectionID("SD11");
  auto Col12 = SDman->GetCollectionID("SD12");
  auto Col13 = SDman->GetCollectionID("SD13");
  auto Col14 = SDman->GetCollectionID("SD14");
  auto Col15 = SDman->GetCollectionID("SD15");
  auto Col16 = SDman->GetCollectionID("SD16");
  auto Col17 = SDman->GetCollectionID("SD17");
  auto Col18 = SDman->GetCollectionID("SD18");
101 102 103 104 105 106 107 108 109 110 111 112 113 114
  auto Col19 = SDman->GetCollectionID("SD19");
  auto Col20 = SDman->GetCollectionID("SD20");
  auto Col21 = SDman->GetCollectionID("SD21");
  auto Col22 = SDman->GetCollectionID("SD22");
  auto Col23 = SDman->GetCollectionID("SD23");
  auto Col24 = SDman->GetCollectionID("SD24");
  auto Col25 = SDman->GetCollectionID("SD25");
  auto Col26 = SDman->GetCollectionID("SD26");
  auto Col27 = SDman->GetCollectionID("SD27");
  auto Col28 = SDman->GetCollectionID("SD28");
  auto Col29 = SDman->GetCollectionID("SD29");
  auto Col30 = SDman->GetCollectionID("SD30");
  auto Col31 = SDman->GetCollectionID("SD31");
  auto Col32 = SDman->GetCollectionID("SD32");
115

116 117

  G4HCofThisEvent* HCE = event->GetHCofThisEvent();
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
  B1HitsCollection* HitsCol1 = 0;
  B1HitsCollection* HitsCol2 = 0;
  B1HitsCollection* HitsCol3 = 0;
  B1HitsCollection* HitsCol4 = 0;
  B1HitsCollection* HitsCol5 = 0;
  B1HitsCollection* HitsCol6 = 0;
  B1HitsCollection* HitsCol7 = 0;
  B1HitsCollection* HitsCol8 = 0;
  B1HitsCollection* HitsCol9 = 0;
  B1HitsCollection* HitsCol10 = 0;
  B1HitsCollection* HitsCol11 = 0;
  B1HitsCollection* HitsCol12 = 0;
  B1HitsCollection* HitsCol13 = 0;
  B1HitsCollection* HitsCol14 = 0;
  B1HitsCollection* HitsCol15 = 0;
  B1HitsCollection* HitsCol16 = 0;
  B1HitsCollection* HitsCol17 = 0;
  B1HitsCollection* HitsCol18 = 0;
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
  B1HitsCollection* HitsCol19 = 0;
  B1HitsCollection* HitsCol20 = 0;
  B1HitsCollection* HitsCol21 = 0;
  B1HitsCollection* HitsCol22 = 0;
  B1HitsCollection* HitsCol23 = 0;
  B1HitsCollection* HitsCol24 = 0;
  B1HitsCollection* HitsCol25 = 0;
  B1HitsCollection* HitsCol26 = 0;
  B1HitsCollection* HitsCol27 = 0;
  B1HitsCollection* HitsCol28 = 0;
  B1HitsCollection* HitsCol29 = 0;
  B1HitsCollection* HitsCol30 = 0;
  B1HitsCollection* HitsCol31 = 0;
  B1HitsCollection* HitsCol32 = 0;

151 152

  if(HCE) {
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
    HitsCol1 = (B1HitsCollection*)(HCE->GetHC(Col1));
    HitsCol2 = (B1HitsCollection*)(HCE->GetHC(Col2));
    HitsCol3 = (B1HitsCollection*)(HCE->GetHC(Col3));
    HitsCol4 = (B1HitsCollection*)(HCE->GetHC(Col4));
    HitsCol5 = (B1HitsCollection*)(HCE->GetHC(Col5));
    HitsCol6 = (B1HitsCollection*)(HCE->GetHC(Col6));
    HitsCol7 = (B1HitsCollection*)(HCE->GetHC(Col7));
    HitsCol8 = (B1HitsCollection*)(HCE->GetHC(Col8));
    HitsCol9 = (B1HitsCollection*)(HCE->GetHC(Col9));
    HitsCol10 = (B1HitsCollection*)(HCE->GetHC(Col10));
    HitsCol11 = (B1HitsCollection*)(HCE->GetHC(Col11));
    HitsCol12 = (B1HitsCollection*)(HCE->GetHC(Col12));
    HitsCol13 = (B1HitsCollection*)(HCE->GetHC(Col13));
    HitsCol14 = (B1HitsCollection*)(HCE->GetHC(Col14));
    HitsCol15 = (B1HitsCollection*)(HCE->GetHC(Col15));
    HitsCol16 = (B1HitsCollection*)(HCE->GetHC(Col16));
    HitsCol17 = (B1HitsCollection*)(HCE->GetHC(Col17));
    HitsCol18 = (B1HitsCollection*)(HCE->GetHC(Col18));
171 172 173 174 175 176 177 178 179 180 181 182 183 184
    HitsCol19 = (B1HitsCollection*)(HCE->GetHC(Col19));
    HitsCol20 = (B1HitsCollection*)(HCE->GetHC(Col20));
    HitsCol21 = (B1HitsCollection*)(HCE->GetHC(Col21));
    HitsCol22 = (B1HitsCollection*)(HCE->GetHC(Col22));
    HitsCol23 = (B1HitsCollection*)(HCE->GetHC(Col23));
    HitsCol24 = (B1HitsCollection*)(HCE->GetHC(Col24));
    HitsCol25 = (B1HitsCollection*)(HCE->GetHC(Col25));
    HitsCol26 = (B1HitsCollection*)(HCE->GetHC(Col26));
    HitsCol27 = (B1HitsCollection*)(HCE->GetHC(Col27));
    HitsCol28 = (B1HitsCollection*)(HCE->GetHC(Col28));
    HitsCol29 = (B1HitsCollection*)(HCE->GetHC(Col29));
    HitsCol30 = (B1HitsCollection*)(HCE->GetHC(Col30));
    HitsCol31 = (B1HitsCollection*)(HCE->GetHC(Col31));
    HitsCol32 = (B1HitsCollection*)(HCE->GetHC(Col32));
185 186
  }

187 188 189 190
 //G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
  std::ofstream gammafile("gamma.txt");
  std::ofstream mu_plusfile("mu+.txt");
  std::ofstream mu_minusfile("mu-.txt");
191 192 193 194 195
  std::ofstream nu_e_file("nu_e.txt");
  std::ofstream nu_mu_file("nu_mu.txt");
  std::ofstream anti_nu_e_file("anti_nu_e.txt");
  std::ofstream anti_nu_mu_file("anti_nu_mu.txt");

196 197


198 199 200 201 202 203 204
  detectorConstruction = (B1DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction();
  G4cout << "tttttt" << G4endl;
  dec_pos = detectorConstruction->GetPositions();
  half_height = detectorConstruction->get_size();
  variation = detectorConstruction->get_delta();


205
  if(HitsCol1) {
206 207 208
    G4cout << "test" << G4endl;
    G4double pos_detec = SetPosition(1);
    G4cout << "RRRRRRR" << G4endl;
209 210
    int n_hit = HitsCol1->entries();
    G4cout << "My detector 1 has " << n_hit << "hits" << G4endl;
211 212
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
213 214 215 216 217 218 219
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol1)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 1);
    }
220
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
221 222
    fparticles.clear();
  }
223

224
  if(HitsCol2) {
225
    G4double pos_detec = SetPosition(2);
226 227 228 229
    int n_hit = HitsCol2->entries();
    G4cout << "My detector 2 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
230
      for(int i1 = 0; i1 < n_hit; i1++) {
231
      B1Hits* hit = (*HitsCol2)[i1];
232 233 234
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
235
      WriteHistogram(name, 2);
236
    }
237
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
238
    fparticles.clear();
239 240
  }

241
  if(HitsCol3) {
242
    G4double pos_detec = SetPosition(3);
243 244 245 246 247 248 249 250 251 252 253
    int n_hit = HitsCol3->entries();
    G4cout << "My detector 3 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol3)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 3);
    }
254
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
255 256 257 258
    fparticles.clear();
  }

  if(HitsCol4) {
259
    G4double pos_detec = SetPosition(4);
260 261 262 263 264 265 266 267 268 269 270
    int n_hit = HitsCol4->entries();
    G4cout << "My detector 4 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol4)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      //WriteHistogram(name, 4);
    }
271
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
272 273 274 275
    fparticles.clear();
  }

  if(HitsCol5) {
276
    G4double pos_detec = SetPosition(5);
277 278 279 280 281 282 283 284 285 286 287
    int n_hit = HitsCol5->entries();
    G4cout << "My detector 5 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol5)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 5);
    }
288
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
289 290 291 292
    fparticles.clear();
  }

  if(HitsCol6) {
293
    G4double pos_detec = SetPosition(6);
294 295 296 297 298 299 300 301 302 303 304
    int n_hit = HitsCol6->entries();
    G4cout << "My detector 6 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol6)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
305
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
306 307 308
    fparticles.clear();
  }
  if(HitsCol7) {
309
    G4double pos_detec = SetPosition(7);
310 311 312 313 314 315 316 317 318 319 320
    int n_hit = HitsCol7->entries();
    G4cout << "My detector 7 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol7)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 7);
    }
321
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
322 323 324
    fparticles.clear();
  }
  if(HitsCol8) {
325
    G4double pos_detec = SetPosition(8);
326 327 328 329 330 331 332 333 334 335 336
    int n_hit = HitsCol8->entries();
    G4cout << "My detector 8 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol8)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 8);
    }
337
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
338 339 340 341
    fparticles.clear();
  }

  if(HitsCol9) {
342
    G4double pos_detec = SetPosition(9);
343 344 345 346 347 348 349 350 351 352 353
    int n_hit = HitsCol9->entries();
    G4cout << "My detector 9 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol9)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 9);
    }
354
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
355 356 357 358
    fparticles.clear();
  }

  if(HitsCol10) {
359
    G4double pos_detec = SetPosition(10);
360 361 362 363 364 365 366 367 368 369 370
    int n_hit = HitsCol10->entries();
    G4cout << "My detector 10 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol10)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 10);
    }
371
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
372 373 374 375
    fparticles.clear();
  }

  if(HitsCol11) {
376
    G4double pos_detec = SetPosition(11);
377 378 379 380 381 382 383 384 385 386 387
    int n_hit = HitsCol11->entries();
    G4cout << "My detector 11 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol11)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 11);
    }
388
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
389 390 391
    fparticles.clear();
  }
  if(HitsCol12) {
392
    G4double pos_detec = SetPosition(12);
393 394 395 396 397 398 399 400 401 402 403
    int n_hit = HitsCol12->entries();
    G4cout << "My detector 12 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol12)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 12);
    }
404
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
405 406 407
    fparticles.clear();
  }
  if(HitsCol13) {
408
    G4double pos_detec = SetPosition(13);
409 410 411 412 413 414 415 416 417 418 419
    int n_hit = HitsCol13->entries();
    G4cout << "My detector 13 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol13)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 13);
    }
420
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
421 422 423
    fparticles.clear();
  }
  if(HitsCol14) {
424
    G4double pos_detec = SetPosition(14);
425
    int n_hit = HitsCol14->entries();
426
    G4cout << "My detector 14 has " << n_hit << "hits" << G4endl;
427 428 429 430 431 432 433 434 435
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol14)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 14);
    }
436
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
437 438 439 440
    fparticles.clear();
  }

  if(HitsCol15) {
441
    G4double pos_detec = SetPosition(15);
442 443 444 445 446 447 448 449 450 451 452
    int n_hit = HitsCol15->entries();
    G4cout << "My detector 15 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol15)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
    //  WriteHistogram(name, 15);
    }
453
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
454 455 456
    fparticles.clear();
  }
  if(HitsCol16) {
457
    G4double pos_detec = SetPosition(16);
458 459 460 461 462 463 464 465 466 467 468
    int n_hit = HitsCol16->entries();
    G4cout << "My detector 16 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol16)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
469
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
470 471 472
    fparticles.clear();
  }
  if(HitsCol17) {
473
    G4double pos_detec = SetPosition(17);
474 475 476 477 478 479 480 481 482 483 484
    int n_hit = HitsCol17->entries();
    G4cout << "My detector 17 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol17)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
485
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
486 487 488
    fparticles.clear();
  }
  if(HitsCol18) {
489
    G4double pos_detec = SetPosition(18);
490 491 492 493 494 495 496 497 498 499 500
    int n_hit = HitsCol18->entries();
    G4cout << "My detector 18 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol18)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
501
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
502 503
    fparticles.clear();
  }
504

505
  if(HitsCol19) {
506
    G4double pos_detec = SetPosition(19);
507 508 509 510 511 512 513 514 515 516 517
    int n_hit = HitsCol19->entries();
    G4cout << "My detector 19 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol19)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol20) {
    G4double pos_detec = SetPosition(20);
    int n_hit = HitsCol20->entries();
    G4cout << "My detector 20 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol20)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol21) {
    G4double pos_detec = SetPosition(21);
    int n_hit = HitsCol21->entries();
    G4cout << "My detector 21 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol21)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol22) {
    G4double pos_detec = SetPosition(22);
    int n_hit = HitsCol22->entries();
    G4cout << "My detector 22 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol22)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol23) {
    G4double pos_detec = SetPosition(23);
    int n_hit = HitsCol23->entries();
    G4cout << "My detector 23 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol23)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }


  if(HitsCol24) {
    G4double pos_detec = SetPosition(24);
    int n_hit = HitsCol24->entries();
    G4cout << "My detector 24 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol24)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol25) {
    G4double pos_detec = SetPosition(25);
    int n_hit = HitsCol25->entries();
    G4cout << "My detector 25 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol25)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol26) {
    G4double pos_detec = SetPosition(26);
    int n_hit = HitsCol26->entries();
    G4cout << "My detector 26 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol26)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol27) {
    G4double pos_detec = SetPosition(27);
    int n_hit = HitsCol27->entries();
    G4cout << "My detector 27 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol27)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol28) {
    G4double pos_detec = SetPosition(28);
    int n_hit = HitsCol28->entries();
    G4cout << "My detector 28 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol28)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol29) {
    G4double pos_detec = SetPosition(29);
    int n_hit = HitsCol29->entries();
    G4cout << "My detector 29 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol29)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol30) {
    G4double pos_detec = SetPosition(30);
    int n_hit = HitsCol30->entries();
    G4cout << "My detector 30 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol30)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol31) {
    G4double pos_detec = SetPosition(31);
    int n_hit = HitsCol31->entries();
    G4cout << "My detector 31 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol31)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
    fparticles.clear();
  }

  if(HitsCol32) {
    G4double pos_detec = SetPosition(32);
    int n_hit = HitsCol32->entries();
    G4cout << "My detector 32 has " << n_hit << "hits" << G4endl;
    B1Hits* hit = new B1Hits;
    std::map<const G4String, int> fparticles;
      for(int i1 = 0; i1 < n_hit; i1++) {
      B1Hits* hit = (*HitsCol32)[i1];
      const G4String name = hit->getParticleInTarget();
      //G4cout << name << G4endl;
      fparticles[name]++;
      WriteHistogram(name, 6);
    }
    PrintParticles(fparticles, gammafile, mu_plusfile, mu_minusfile, nu_e_file, nu_mu_file, anti_nu_e_file, anti_nu_mu_file, pos_detec);
741 742 743 744 745 746 747
    fparticles.clear();
  }


gammafile.close();
mu_plusfile.close();
mu_minusfile.close();
748 749 750 751
nu_e_file.close();
nu_mu_file.close();
anti_nu_e_file.close();
anti_nu_mu_file.close();
752 753
}

754 755
void B1EventAction::PrintParticles(std::map<const G4String, int>& container, std::ofstream& gamma, std::ofstream& mu_minus, std::ofstream& mu_plus,
                                   std::ofstream& nu_e, std::ofstream& nu_mu, std::ofstream& anti_nu_e, std::ofstream& anti_nu_mu, G4double position) {
756
  std::map<const G4String, int>::iterator it;
757
//  G4cout << "Número de párticulas identificadas no detetor: " << G4endl;
758
   G4bool isGamma, ismu_minus, ismu_plus, isnu_e, isnu_mu, isanti_nu_e, isanti_nu_mu = false;
759

760 761 762
    for(it = container.begin() ;it != container.end(); it ++)
    {
      G4cout << "N " << it->first << " : " << it->second << G4endl;
763
      //FUTURO -> substituir por switch;
764 765 766 767 768 769 770 771 772
      if(it->first == "gamma") {
        gamma << it->second << " " << position << "\n";
        isGamma = true;
      } else if (it->first == "mu+") {
        mu_plus << it->second << " " << position << "\n";
        ismu_plus= true;
      } else if (it->first == "mu-") {
        mu_minus << it->second << " " << position << "\n";
        ismu_minus = true;
773 774 775 776 777 778 779 780 781 782 783 784
      } else if (it->first == "nu_e") {
        nu_e << it->second << " " << position << "\n";
        isnu_e = true;
      } else if (it->first == "nu_mu") {
        nu_mu << it->second << " " << position << "\n";
        isnu_mu = true;
      } else if (it->first == "anti_nu_e") {
        anti_nu_e << it->second << " " << position << "\n";
        isanti_nu_e = true;
      } else if (it->first == "anti_nu_mu") {
        anti_nu_mu << it->second << " " << position << "\n";
        isanti_nu_mu = true;
785 786 787 788 789 790 791 792 793 794 795 796
      }
  }

  if (!isGamma) {
    gamma << "0 " << position << "\n";
  }
  if (!ismu_plus) {
    mu_plus << "0 " << position << "\n";
  }
  if (!ismu_minus) {
    mu_minus << "0 " << position << "\n";
  }
797 798 799 800 801 802 803 804 805 806 807 808
  if (!isnu_e) {
    nu_e << "0 " << position << "\n";
  }
  if (!isnu_mu) {
    nu_mu << "0 " << position << "\n";
  }
  if (!isanti_nu_e) {
    anti_nu_e << "0 " << position << "\n";
  }
  if (!isanti_nu_mu) {
    anti_nu_mu << "0 " << position << "\n";
  }
809 810
}

811 812 813 814 815 816 817 818 819 820 821 822 823
void B1EventAction::WriteHistogram(const G4String name, G4int Detec) {
  /*
   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
  if (name == "gamma") {
    analysisManager->FillNtupleDColumn(2, Detec);
  } else if (name == "mu+") {
    analysisManager->FillNtupleDColumn(1, Detec);
  } else if (name == "mu-") {
    analysisManager->FillNtupleDColumn(0, Detec);
  }
  analysisManager->AddNtupleRow();
*/
}
824

825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
G4double B1EventAction::SetPosition(G4int detec) {
  G4int detec_index = detec - 1;
  G4double first_value = dec_pos[0]/1000000;
  G4cout << first_value << G4endl;
  G4double value =  dec_pos[detec_index];
  G4double value_c = value/1000000;
  G4double half_height_real = half_height/1000000;
  G4double h_real = half_height_real*2;
  G4double real_delta = variation/1000000;
  G4double first_distance = half_height_real + first_value ;
  if (detec_index == 0) {
    return first_distance;
  }
  return (first_distance + detec_index*real_delta);
}


842
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......