[Epic] update macro

Projects/Epic_sim_proto/root/simulation/ReadGefEpicSimulatedTree.C

@@ -84,10 +84,6 @@ void run()
   TCanvas * can1D_QvsT = new TCanvas("1D_QvsT","One Event",0,0,2500,2500);
 #endif
 
-
-
-
-
   // === ============================================================================
   // === Loop over the entries
   unsigned int nentries = t->GetEntries();
@@ -116,14 +112,71 @@ void run()
 
 
     if(m_Data->GetMultiplicity()==0) continue;
+
     else if(m_Data->GetMultiplicity()!=1) {
-      cout << "Particle Multiplicity = " << m_Data->GetMultiplicity() << endl;
-      time_window_ns = 0;
-      // FIXME : LOOP
-      // loop over the particle multiplicity
+#if DEBUG
+      cout << "CASE OF 2FF : m_Data->GetMultiplicity() = " << m_Data->GetMultiplicity() << endl;
+#endif
+      double time_window_ns_max = 0;
+      for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
+        int    nsteps  = m_Data->GetNumberOfSteps(pmult);
+        double tfirst  = m_Data->GetTimeCreationElectronsPerStep(pmult).at(0);
+        double tlast   = m_Data->GetTimeCreationElectronsPerStep(pmult).at(nsteps-1);
+        time_window_ns = max(tfirst+time_KtoA_ns,tlast+(114-nsteps)*time_bin_width_ns);
+        if(time_window_ns>time_window_ns_max) time_window_ns_max = time_window_ns;
+      }
+      nbins_QvsT     = ceil( time_window_ns_max / time_bin_width_ns );
+      nbins_convol   = nbins_QvsT + 2*nbins_gauss;
+      vector<double> influence(nbins_convol,0);
+#if DEBUG
+      TH1D * h1_QvsT = new TH1D("influence_vs_time_1D","influence_vs_time_1D",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
+      h1_QvsT->GetXaxis()->SetTitle("Absolute time [ps] 200ps/bin"); 
+      h1_QvsT->SetLineColor(kBlue);
+      h1_QvsT->SetDirectory(0);
+      TH1D * h1_convolution_QvsT = new TH1D("QvsT_convol","QvsT_convol",nbins_convol,-2*range_gauss,nbins_QvsT*1000*time_bin_width_ns+2*range_gauss);
+      h1_convolution_QvsT->SetLineColor(kRed);
+#endif
       for(unsigned short pmult=0; pmult<m_Data->GetMultiplicity(); pmult++){
+        int    nsteps  = m_Data->GetNumberOfSteps(pmult);
+        for(int smult=0; smult<nsteps; smult++){
+          double  t_Step = m_Data->GetTimeCreationElectronsPerStep(pmult).at(smult) ; // [ns]
+          int     n_Step = m_Data->GetNumElectronsPerStep(pmult).at(smult);
+          int     b_Step = t_Step/time_bin_width_ns; 
+          for(int step=smult; step<114; step++){
+            if(nbins_QvsT < b_Step+step-smult){
+              cout << "size of the vector: " << nbins_QvsT 
+                   << "  but index to fill vector is " << b_Step+step-smult << endl;
+            }
+            else {
+              influence[b_Step+step-smult+nbins_gauss] += n_Step * 22. / 2500. ;
+            }
+          }// end of loop over the steps up to anodes  
+        }// end of loop over the steps
+      }// end of loop over the particle multiplicity (if 2 FF)
+      double Qmax = 0;
+      for(int i=0; i<nbins_QvsT+nbins_gauss; i++){
+        double convolution = 0;
+        for (int j = 0; j < nbins_gauss; j++) convolution += influence[i+j] * gaussian_distribution[j]; 
+#if DEBUG
+        h1_QvsT->SetBinContent(i+1,influence[i]);
+        h1_QvsT->SetBinError(i+1,0);
+        h1_convolution_QvsT->SetBinContent(i+1+0.5*nbins_gauss,convolution);
+#endif
+        if (convolution>Qmax) Qmax=convolution;
       }
+      h1_Qmax->Fill(Qmax);
+      influence.clear();
+#if DEBUG
+      can1D_QvsT->cd(); 
+      h1_QvsT->Draw();
+      h1_convolution_QvsT->Draw("same");
+      can1D_QvsT->Update();
+      can1D_QvsT->WaitPrimitive();
+      h1_QvsT->Reset();
+      h1_convolution_QvsT->Reset();
+#endif
     }// end of if else mult>1
+    
     else{
       int    nsteps  = m_Data->GetNumberOfSteps(0);
       double tfirst  = m_Data->GetTimeCreationElectronsPerStep(0).at(0);
@@ -207,7 +260,7 @@ void run()
   can_Qmax->cd();
   h1_Qmax->Draw();
 
-//  TFile * fsave = new TFile(Form("/media/audrey/DATA1/EPIC/simulation/out_K11_alpha240Pu_window%ius_%i.root",AlphaDecay_TimeWindow,fileno),"recreate");
+//  TFile * fsave = new TFile(Form("/media/audrey/DATA1/EPIC/simulation/out_K11_GEF240Pu.root",AlphaDecay_TimeWindow,fileno),"recreate");
 //  h1_Qmax->Write();
 //  fsave->Close();