Merge pull request #12 from lopezzot/dev

From dev to master for v2.1 release

ATLHECTB.cc

@@ -1,7 +1,8 @@
 //**************************************************
 // \file ATLHECTB.cc
 // \brief: main() of ATLHECTB project
-// \author: Lorenzo Pezzotti (CERN EP-SFT-sim) @lopezzot
+// \author: Lorenzo Pezzotti (CERN EP-SFT-sim) 
+//          @lopezzot
 // \start date: 11 May 2021
 //**************************************************
 
@@ -12,7 +13,12 @@
 
 //Includers from Geant4
 //
-#include "G4RunManagerFactory.hh"
+//#include "G4RunManagerFactory.hh" //G4RunManagerFactory is only available from 10.7
+#ifdef G4MULTITHREADED
+#include "G4MTRunManager.hh"
+#else
+#include "G4RunManager.hh"
+#endif
 #include "G4UImanager.hh"
 #include "G4UIcommand.hh"
 #include "Randomize.hh"
@@ -74,11 +80,22 @@ int main( int argc, char** argv ) {
         ui = new G4UIExecutive(argc, argv, session);
     }
 
-    //Construct run manager (default of multithreaded)
-    auto* runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
+    //Construct run manager (default type)
+    /*uncomment this part for G4RunManagerFactory usage (10.7 on)
+		auto* runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
     #ifdef G4MULTITHREADED
     if ( nthreads > 0 ) runManager->SetNumberOfThreads(nthreads);
     #endif
+    */
+    //Construct run manager (works also for versions before 10.7)
+    #ifdef G4MULTITHREADED
+    auto runManager = new G4MTRunManager;
+    if ( nthreads > 0 ) { 
+    		runManager->SetNumberOfThreads(nthreads);
+    }  
+    #else
+    auto runManager = new G4RunManager;
+    #endif
 
     //Set mandatory classes (DetConstruction, PhysicsList, ActionInitialization)
     //
@@ -127,7 +144,3 @@ int main( int argc, char** argv ) {
 }
 
 //**************************************************
-
-
-
-

README.md

@@ -128,6 +128,7 @@ The following are results deployed on Geant Val so far.
 
 | ATLHECTB         | Reproduce data | Reproduce analysis | Comments     |
 | -------------    | ----------     | -----------        | -----------  |
+| v2.1 <br /> Dataset #1 <br /> tag 2.1_1 (Geant4.10.06.p01 and Geant4.11.0.beta, ATLHECTB v2.1, FTFP_BERT, FTFP_BERT_ATL, QGSP_BERT, FTFP_BERT_INCLXX) <br /> Added on 21/10/2021 <br /> | python mc-config-generator.py submit -t ATLHECTB -d OUTPUT -v 10.6.p01 (11.0.beta) -q "testmatch" -r | python mc-config-generator.py parse -t ATLHECTB -d OUTPUT |  Results for two Geant4 versions, can be compared to tag 2.0_1 and 2.0_2.|
 | v2.0 <br /> Dataset #2 <br /> tag 2.0_2 (Geant4.10.07.p01, ATLHECTB v2.0, FTFP_BERT_ATL, QGSP_BERT, FTFP_BERT_INCLXX) <br /> Added on 18/10/2021 <br /> | python mc-config-generator.py submit -t ATLHECTB -d OUTPUT -v 10.7.p01 -q "testmatch" -r | python mc-config-generator.py parse -t ATLHECTB -d OUTPUT | Everything identical to 2.0_1 but changed physics lists in params.conf (change it yourself to reproduce data).|
 | v2.0 <br /> Dataset #1 <br /> tag 2.0_1 (Geant4.10.07.p01, ATLHECTB v2.0, FTFP_BERT) <br /> Added on 11/10/2021 <br /> | python mc-config-generator.py submit -t ATLHECTB -d OUTPUT -v 10.7.p01 -q "testmatch" -r | python mc-config-generator.py parse -t ATLHECTB -d OUTPUT | First results on Geant Val, using Geant4.10.07.p01, ATLHECTB v2.0, FTFP_BERT. Analysis coded in parser.py. JSON files for test-beam data are created with parser.py (end of file). |
 

geantval_scripts/ATLHECTB/params.conf

@@ -1,4 +1,4 @@
-!PHYSLIST=FTFP_BERT
+!PHYSLIST=FTFP_BERT, QGSP_BERT, FTFP_BERT_ATL, FTFP_INCLXX
 !CONST:ENERGY_UNIT=GeV
 PARTICLE | ENERGY | PHYSLIST | NEVENTS
 pi- | 20. | PHYSLIST | 50000

geantval_scripts/ATLHECTB/parser.py

@@ -40,10 +40,10 @@ def parse(self, jobs):
 	pijobs = [x for x in jobs if x["PARTICLE"]=="pi-"]
 	pienergies = [float(x["ENERGY"]) for x in pijobs]
 	pifiles = [os.path.join(x["path"],"ATLHECTBout_Run0.root") for x in pijobs]
-	print "Found "+str(len(jobs))+" runs in jobs with:"
-	print "--->"+ str(len(ectrjobs)) + " jobs with e-, energies (GeV): ", ectrenergies
+	print "Found "+str(len(jobs))+" runs in jobs:"
+	print "--->"+ str(len(ectrjobs)) + " jobs with e-, energies (GeV): " + str(ectrenergies) + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
 	#print "------> files: ", ectrfiles
-	print "--->"+ str(len(pijobs)) + " jobs with pi-, energies (GeV): ", pienergies
+	print "--->"+ str(len(pijobs)) + " jobs with pi-, energies (GeV): " + str(pienergies) + " ,physlist: " + str(set([x["PHYSLIST"] for x in pijobs]))
 	#print "------> files: ", pifiles
 
 	#e- analysis
@@ -63,9 +63,9 @@ def parse(self, jobs):
 		H1sampfraction.Fill(evt.elAr/(energy*1000)*100) #percent value
 	    sampfractions.append(H1sampfraction.GetMean())
 	    ersampfractions.append(H1sampfraction.GetMeanError())
-	print "--->e- sampling fraction: " + str(sampfractions)
-	print "--->e- avg sampling fraction: " + str(np.mean(sampfractions)) + "%"
-	#outfile = TFile.Open("OUTe-.root","RECREATE")
+	print "--->e- sampling fraction: " + str(sampfractions) + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
+	print "--->e- avg sampling fraction: " + str(np.mean(sampfractions)) + "%" + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
+	outfile = TFile.Open("OUTe-.root","RECREATE")
 	for energy in eenergies:
 	    #Find e- job with corresponding energy
 	    job = [x for x in ectrjobs if float(x["ENERGY"])==energy][0]
@@ -91,10 +91,10 @@ def parse(self, jobs):
 	    erres = (recenergy.GetFunction("gaus").GetParError(2)/
 		     recenergy.GetFunction("gaus").GetParameter(2))*res
 	    erresolutions.append(erres)
-	    #outfile.cd()
-	    #recenergy.Write()
-	print "--->e- sampling terms in resolution: " + str(resolutions) + " %GeV^{1/2}"
-	print "--->e- avg sampling term in resolution: " + str(np.mean(resolutions)) + " %GeV^{1/2}" 
+	    outfile.cd()
+	    recenergy.Write()
+	print "--->e- sampling terms in resolution: " + str(resolutions) + " %GeV^{1/2}" + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
+	print "--->e- avg sampling term in resolution: " + str(np.mean(resolutions)) + " %GeV^{1/2}" + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
 
         #Create JSON output files for e- energy resolution (graph)
         #
@@ -274,13 +274,14 @@ def parse(self, jobs):
 		addchannelF4 += 2.*evt.M3L4BirkeLayer[7] #M3L4 
 		addchannelF4 += 2.*evt.M3L4BirkeLayer[9] 
 		addchannelF4 += 2.*evt.M3L4BirkeLayer[11] 
-
-		response.Fill((addchannel/energy)/(10.*np.mean(sampfractions)))
-		recenergy.Fill(addchannel/(10.*np.mean(sampfractions)))
-		H1F1.Fill(addchannelF1/addchannel)
-		H1F2.Fill(addchannelF2/addchannel)
-		H1F3.Fill(addchannelF3/addchannel)
-		H1F4.Fill(addchannelF4/addchannel)
+
+		if addchannel>0.:
+		    response.Fill((addchannel/energy)/(10.*np.mean(sampfractions)))
+		    recenergy.Fill(addchannel/(10.*np.mean(sampfractions)))
+		    H1F1.Fill(addchannelF1/addchannel)
+		    H1F2.Fill(addchannelF2/addchannel)
+		    H1F3.Fill(addchannelF3/addchannel)
+		    H1F4.Fill(addchannelF4/addchannel)
 
 	    responses.append(response.GetMean())
 	    erresponses.append(response.GetMeanError())
@@ -303,7 +304,7 @@ def parse(self, jobs):
 	    residual = (depths[0]*H1F1.GetMean())**2+(depths[1]*H1F2.GetMean())**2+(depths[2]*H1F3.GetMean())**2+(depths[3]*H1F4.GetMean())**2
 	    sigmaL0.append(2.*(residual/4.)**0.5) 
 
-	print "--->pi- pi/e: "+ str(responses) 
+	print "--->pi- pi/e: "+ str(responses) + " ,physlist: " + str(set([x["PHYSLIST"] for x in ectrjobs]))
 	print "--->pi- avg pi/e: "+ str(np.mean(responses))
 	print "--->pi- resolutions: " + str(resolutions) + " %"
 	print "--->pi- avg resolution: " + str(np.mean(resolutions)) + " %"