Please enter the commit message for your changes. Lines starting

NPSimulation/Detectors/TACTIC/GARFDRIFT.h

@@ -8,28 +8,34 @@
 #include "Garfield/ViewSignal.hh"
 #include "Garfield/ViewMedium.hh"
 		 
-void GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector delta_pos, double R, int Pad_start, double ID,
-	  double ScoLength, double SegLength) {
+double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector delta_pos, double R, int Pad_start, double ID,
+	       double ScoLength, double SegLength, double event) {
 
 
   const double rWire = 1.2;
   const double rTube = 5.;
   const double lTube = 25.19;
+  //TH1F *hist = new TH1F("","",1000,0,10000); // 10 ns bins
   //For R6 = 2, R7 = 1
   //const double vWire = -375.; //For V_total = 750V
   //const double vWire = -645.3312; //For V_total = 1200V //SHOULD BE 583.741 !?
   //const double vWire = -535.0959; //For V_total = 1100V
   //const double vWire = -875.6116; //For V_total = 1800V 
-  const double vWire =  -681.0312; //For V_total = 1400V
+  //const double vWire =  -681.0312; //For V_total = 1400V
+  //const double vWire = -730.; 
+  const double vWire = -1000.;
   const double vTube = 0.;
   double x_start, y_start, z_start;
   ofstream file;
-  double production_bias = 1.e-01;
+  double production_bias = 0.001;
   double t_end, x_end, y_end, z_end;
   double Pad_end;
   //double w_value = 35.; //for HeCO2 --guess
-  double w_value = 26.31; //for argon
+  //double w_value = 26.31; //for argon
+  double w_value = 41.1;// for He
   int electrons  = (int)(energy/w_value*production_bias);//*production_bias; //Need to carry excess energy over later.
+  double energy_excess = 0.;
+  vector<double> pad_vec, time_vec; 
   
   Garfield::MediumMagboltz* gas = new Garfield::MediumMagboltz();
   Garfield::ViewMedium* mediumView = new Garfield::ViewMedium();
@@ -42,10 +48,15 @@ void GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector
   cout << "\n" << electrons << "\n" << endl;
   
   if(R>rWire && electrons>0) {
-
-    gas->LoadGasFile("ar_90_ch4_10.gas"); //atomic fraction in this case (P10).
+    
+    energy_excess = (energy/w_value*production_bias - electrons)*w_value/production_bias;
+    
+    //gas->LoadGasFile("he_90_co2_10_200mbar.gas");
+    //gas->LoadGasFile("ar_90_ch4_10.gas"); //atomic fraction in this case (P10).
     //gas->LoadGasFile("he_90_co2_10_500mbar.gas"); //mass fraction in this case
-   
+    //gas->LoadGasFile("he_90_co2_10_1bar.gas");
+    gas->LoadGasFile("he_90_co2_10_100mbar.gas");
+    
     gas->Initialise(true);
     
     mediumView->SetMedium(gas);
@@ -69,14 +80,18 @@ void GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector
       drift->GetDriftLinePoint(np-1, x_end, y_end, z_end, t_end); //np-1 is last DriftLineEndPoint
       
       Pad_end = (int)((z_end + ScoLength / 2.) / SegLength ) + 1; //new Pad number 
-      
-      file.open("signal.dat", std::ios::app);
-      file <<  Pad_end << "\t" << t_end << endl;
-      file.close();
-      
+
+      if(pow(pow(x_end,2)+pow(y_end,2),0.5)>4.9) { //reached the anode
+	file.open("signal.dat", std::ios::app);
+	file << event << "\t" << ID << "\t" << Pad_start << "\t" <<  Pad_end << "\t" << t_end << endl;
+	file.close();
+      }
+    
     }
   }
-
+  
   delete gas; delete mediumView; delete geo; delete tube; delete cmp; delete sensor; delete drift; //otherwise these are held in memory and cause a killed: 9 crash.
+
+  return energy_excess;
   
 }

NPSimulation/Detectors/TACTIC/TACTIC.cc

@@ -404,12 +404,13 @@ void TACTIC::InitializeScorers() {
   */
     LightFilter->addIon(m_Reaction.GetParticle3()->GetZ(),m_Reaction.GetParticle3()->GetA());
     HeavyFilter->addIon(m_Reaction.GetParticle4()->GetZ(),m_Reaction.GetParticle4()->GetA());
-    LightScorer->SetFilter(LightFilter);
-    HeavyScorer->SetFilter(HeavyFilter);
     if(m_Reaction.GetParticle1()->GetZ() == m_Reaction.GetParticle4()->GetZ()) BeamFilter->add("geantino");
     else BeamFilter->addIon(m_Reaction.GetParticle1()->GetZ(),m_Reaction.GetParticle1()->GetA());
-    BeamScorer->SetFilter(BeamFilter);
     
+    LightScorer->SetFilter(LightFilter);
+    HeavyScorer->SetFilter(HeavyFilter);
+    BeamScorer->SetFilter(BeamFilter);
+
     m_Scorer->RegisterPrimitive(LightScorer);
     m_Scorer->RegisterPrimitive(HeavyScorer);
     m_Scorer->RegisterPrimitive(BeamScorer);

NPSimulation/Detectors/TACTIC/TACTICScorer.cc

@@ -29,7 +29,7 @@ Gas_Scorer::~Gas_Scorer(){}
 
 G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
 
-  G4double* Infos = new G4double[14];
+  G4double* Infos = new G4double[15];
   m_Position  = aStep->GetPreStepPoint()->GetPosition();
 
   Infos[0] = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
@@ -60,18 +60,30 @@ G4bool Gas_Scorer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
     aStep->GetTrack()->SetTrackStatus(fStopAndKill);
     return 0;
   }
-            
-#ifdef USE_Garfield
-  G4ThreeVector delta_Position = aStep->GetDeltaPosition();
-
-  GARFDRIFT(Infos[5]/eV, Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0], Infos[2], Infos[11]);
-#endif
-  
+    
   map<G4int, G4double**>::iterator it;
   it= EvtMap->GetMap()->find(m_Index);
   if(it!=EvtMap->GetMap()->end()){
     G4double* dummy = *(it->second);
     if(Infos[1]==dummy[1]) Infos[5]+=dummy[5]; //accumulate ionisation energy deposit to get total accross pad
+
+          
+#ifdef USE_Garfield
+    G4ThreeVector delta_Position = aStep->GetDeltaPosition();
+    G4double excess;
+
+    if(aStep->IsFirstStepInVolume() == 1) excess = 0.;
+    if(aStep->IsFirstStepInVolume() == 0) excess = dummy[14]/eV;
+    
+    if(excess > 1.e06) excess = 0.; //If dummy[12] is not defined returns random value > 1.e06 ? 
+  
+    if(Infos[8] > 12.) Infos[14] = GARFDRIFT((Infos[5]/eV+excess), Infos[3], m_Position/cm, delta_Position/cm, Infos[8]/cm, Infos[6], Infos[2], m_ScorerLength/cm, m_SegmentLength/cm, Infos[0])*eV;
+#endif
+    if(Infos[8] < 12.) Infos[14] = 0;
+    
+    
+    if(Infos[14]>0.) cout << "Infos[14] " << Infos[14]/eV << endl;
+    
     delete dummy;
   }