Generic analysis

analysis/generic/ATLHECTBanalysis.C

@@ -0,0 +1,83 @@
+//**************************************************
+// \file ATLHECTBanalysis.C
+// \brief: Generic analysis of ATLHECTB data 
+// \author: Lorenzo Pezzotti (CERN EP-SFT-sim)
+//          @lopezzot
+// \start date: 17 August 2023
+//**************************************************
+
+//Usage: root 'ATLHECTBanalysis.C(folderpath/,version,g4version,physicslist)'
+//e.g. root -l 'ATLHECTBanalysis.C("Data1/","2.6","11.1.ref05","FTFP_BERT")'
+//
+//If hadronic data are generated with the Fluka interface add "true"
+//as the last macro argument
+
+//Includers from C++
+//
+#include <string>
+#include <array>
+
+//Includers from project files
+//
+#include "emanalysis.h"
+#include "pianalysis.h"
+#include "ecalibrate.h"
+#include "picalibrate.h"
+
+void ATLHECTBanalysis(const string folder,const string version,const string g4version,const string pl, const bool isFluka=false){
+
+    //This is a generic analysis therefore the legend for MC data
+    //must be taken as input
+    //e.g. "#splitline{ATLHECTB v2.6 }{Geant4.11.1.ref05 FTFP_BERT }"
+    //
+    const string emMClegend = "#splitline{ATLHECTB "+version+ "}{Geant4-"+g4version+" "+pl+ "}";
+    string piMClegend = emMClegend;
+    if(isFluka) piMClegend = "#splitline{ATLHECTB "+version+ "}{FLUKA}";
+
+    // Analysis of e- data
+    // energies 20, 40, 50, 80, 100, 119.1, 147.8 GeV
+    //
+    vector<double> emenergies = {20.,40.,50.,80.,100.,119.1,147.8};
+    vector<string> emfiles;
+    for ( unsigned int i=11; i<18; i++ ){
+        emfiles.push_back( "ATLHECTBout_Run"+std::to_string(i)+".root" );
+    }
+    //First call to emanalysis to extract emdata, using dummy 44.95 value for SF
+    //will be called again later with correct value
+    emoutput emdata = emanalysis( emenergies, emfiles, folder, emMClegend );
+    //Second call to emanalysis with correct SF value, final plots are
+    //now recreated and correct
+    emdata = emanalysis( emenergies, emfiles, folder, emMClegend, emdata.avgSF*10. );
+    emdata.print();
+
+    //Reconstrcuted energies for em runs
+    //For missing energy points (30, 60, 120, 180 and 200 GeV) using 0.99*beamenergy
+    //
+    vector<double> recemenergies = {emdata.recenergies.at(0),0.99*30.,emdata.recenergies.at(1),
+        emdata.recenergies.at(2),0.99*60.,emdata.recenergies.at(3),emdata.recenergies.at(4),
+        emdata.recenergies.at(5),emdata.recenergies.at(6),0.99*180.,0.99*200}; 
+
+    // Analysis of pi- data
+    // energies 20, 30, 40, 50, 60, 80, 100, 120, 150, 180, 200 GeV
+    //
+    vector<double> pienergies = {20.,30.,40.,50.,60.,80.,100.,120.,150.,180.,200.};
+    vector<string> pifiles;
+    for ( unsigned int i=0; i<11; i++ ){
+        pifiles.push_back( "ATLHECTBout_Run"+std::to_string(i)+".root" );
+    }
+    //Last argument is FLUKA(bool). If true assumes the g4-to-fluka interface
+    //was used to generate the data.
+    //
+    pianalysis( pienergies, pifiles, recemenergies, folder, emdata.avgSF*10., piMClegend, isFluka);
+
+    //Analysis to select channels for pi- analysis
+    //
+    //picalibrate(180., "ATLHECTBout_Run9.root");
+
+    //Analysis to select channels for e- analysis
+    //
+    //ecalibrate(147.8,"ATLHECTBout_Run17.root"); 
+
+}
+
+//**************************************************

analysis/generic/ecalibrate.h

@@ -0,0 +1,123 @@
+//**************************************************
+// \file ecalibrate.h
+// \brief: Analysis #1 of ATLHECTB v2.6 
+//         for e- channels 
+// \author: Lorenzo Pezzotti (CERN EP-SFT-sim)
+//          @lopezzot
+// \start date: 10 July 2023
+//**************************************************
+
+#ifndef ecalibrate_H
+#define ecalibrate_H
+
+void ecalibrate( const double& eenergy, const string& efile ){
+
+    //Initiate objects through all the analysis
+    //
+    cout<<"ATLHECTB analysis of channels to be selected with e- runs"<<endl;
+    cout<<"---> Analysis at energy(GeV) "<<eenergy<<endl;  
+
+    string filename = "Data1/"+efile;
+    TFile* file = TFile::Open( filename.c_str(), "READ" );
+    TTree* tree = (TTree*)file->Get( "ATLHECTBout" );
+
+    vector<double>* M1L1BelAr = NULL; 
+    tree->SetBranchAddress( "M1L1BirkeLayer", &M1L1BelAr );
+    vector<double>* M1L2BelAr = NULL; 
+    tree->SetBranchAddress( "M1L2BirkeLayer", &M1L2BelAr );
+    vector<double>* M1L3BelAr = NULL; 
+    tree->SetBranchAddress( "M1L3BirkeLayer", &M1L3BelAr );
+    vector<double>* M1L4BelAr = NULL; 
+    tree->SetBranchAddress( "M1L4BirkeLayer", &M1L4BelAr );
+    vector<double>* M2L1BelAr = NULL; 
+    tree->SetBranchAddress( "M2L1BirkeLayer", &M2L1BelAr );
+    vector<double>* M2L2BelAr = NULL; 
+    tree->SetBranchAddress( "M2L2BirkeLayer", &M2L2BelAr );
+    vector<double>* M2L3BelAr = NULL; 
+    tree->SetBranchAddress( "M2L3BirkeLayer", &M2L3BelAr );
+    vector<double>* M2L4BelAr = NULL; 
+    tree->SetBranchAddress( "M2L4BirkeLayer", &M2L4BelAr );
+    vector<double>* M3L1BelAr = NULL; 
+    tree->SetBranchAddress( "M3L1BirkeLayer", &M3L1BelAr );
+    vector<double>* M3L2BelAr = NULL; 
+    tree->SetBranchAddress( "M3L2BirkeLayer", &M3L2BelAr );
+    vector<double>* M3L3BelAr = NULL; 
+    tree->SetBranchAddress( "M3L3BirkeLayer", &M3L3BelAr );
+    vector<double>* M3L4BelAr = NULL; 
+    tree->SetBranchAddress( "M3L4BirkeLayer", &M3L4BelAr );
+
+    double M1L1avg[24]; memset( M1L1avg, 0., 24*sizeof(double));
+    double M2L1avg[24]; memset( M2L1avg, 0., 24*sizeof(double));
+    double M3L1avg[24]; memset( M3L1avg, 0., 24*sizeof(double));
+
+    double M1L2avg[23]; memset( M1L2avg, 0., 23*sizeof(double));
+    double M2L2avg[23]; memset( M2L2avg, 0., 23*sizeof(double));
+    double M3L2avg[23]; memset( M3L2avg, 0., 23*sizeof(double));
+
+    double M1L3avg[21]; memset( M1L3avg, 0., 21*sizeof(double));
+    double M2L3avg[21]; memset( M2L3avg, 0., 21*sizeof(double));
+    double M3L3avg[21]; memset( M3L3avg, 0., 21*sizeof(double));
+
+    double M1L4avg[20]; memset( M1L4avg, 0., 20*sizeof(double));
+    double M2L4avg[20]; memset( M2L4avg, 0., 20*sizeof(double));
+    double M3L4avg[20]; memset( M3L4avg, 0., 20*sizeof(double));
+
+    //loop over events
+    //
+    for (unsigned int evtNo = 0; evtNo<tree->GetEntries(); evtNo++){
+        tree->GetEntry(evtNo);
+
+        for (unsigned int i = 0; i<24; i++){
+            M2L1avg[i] += M2L1BelAr->at(i)/tree->GetEntries();
+            M1L1avg[i] += M1L1BelAr->at(i)/tree->GetEntries();
+            M3L1avg[i] += M3L1BelAr->at(i)/tree->GetEntries();
+        }
+        for (unsigned int i = 0; i<23; i++){
+            M2L2avg[i] += M2L2BelAr->at(i)/tree->GetEntries();
+            M1L2avg[i] += M1L2BelAr->at(i)/tree->GetEntries();
+            M3L2avg[i] += M3L2BelAr->at(i)/tree->GetEntries();
+        }
+        for (unsigned int i = 0; i<21; i++){
+            M2L3avg[i] += M2L3BelAr->at(i)/tree->GetEntries();
+            M1L3avg[i] += M1L3BelAr->at(i)/tree->GetEntries();
+            M3L3avg[i] += M3L3BelAr->at(i)/tree->GetEntries();
+        }
+        for (unsigned int i = 0; i<20; i++){
+            M2L4avg[i] += M2L4BelAr->at(i)/tree->GetEntries();
+            M1L4avg[i] += M1L4BelAr->at(i)/tree->GetEntries();
+            M3L4avg[i] += M3L4BelAr->at(i)/tree->GetEntries();
+        }
+    }
+
+    double ecut = 5.0;
+    int channels = 0;
+
+    cout<<"List of channels with avg signal above cut"<<endl;
+    for (unsigned int i = 0; i<24; i++){
+        if (M2L1avg[i]>ecut){cout<<"M2L1 "<<i<<" "<<M2L1avg[i]<<endl;channels=channels+1;}
+        if (M1L1avg[i]>ecut){cout<<"M1L1 "<<i<<" "<<M1L1avg[i]<<endl;channels=channels+1;}
+        if (M3L1avg[i]>ecut){cout<<"M3L1 "<<i<<" "<<M3L1avg[i]<<endl;channels=channels+1;}
+    }
+    for (unsigned int i = 0; i<M2L2BelAr->size(); i++){
+        if (M2L2avg[i]>ecut){cout<<"M2L2 "<<i<<" "<<M2L2avg[i]<<endl;channels=channels+1;}
+        if (M1L2avg[i]>ecut){cout<<"M1L2 "<<i<<" "<<M1L2avg[i]<<endl;channels=channels+1;}
+        if (M3L2avg[i]>ecut){cout<<"M3L2 "<<i<<" "<<M3L2avg[i]<<endl;channels=channels+1;}
+    }
+    for (unsigned int i = 0; i<M2L3BelAr->size(); i++){
+        if (M2L3avg[i]>ecut){cout<<"M2L3 "<<i<<" "<<M2L3avg[i]<<endl;channels=channels+1;}
+        if (M1L3avg[i]>ecut){cout<<"M1L3 "<<i<<" "<<M1L3avg[i]<<endl;channels=channels+1;}
+        if (M3L3avg[i]>ecut){cout<<"M3L3 "<<i<<" "<<M3L3avg[i]<<endl;channels=channels+1;}
+    }
+    for (unsigned int i = 0; i<M2L4BelAr->size(); i++){
+        if (M2L4avg[i]>ecut){cout<<"M2L4 "<<i<<" "<<M2L4avg[i]<<endl;channels=channels+1;}
+        if (M1L4avg[i]>ecut){cout<<"M1L4 "<<i<<" "<<M1L4avg[i]<<endl;channels=channels+1;}
+        if (M3L4avg[i]>ecut){cout<<"M3L4 "<<i<<" "<<M3L4avg[i]<<endl;channels=channels+1;}
+    }
+    cout<<"Number of channels above cut: "<<channels<<endl;
+    cout<<"take first 7 with highest average signal"<<endl;
+
+}
+
+#endif
+
+//**************************************************