Elad michael input file patch

Construct.dat

@@ -0,0 +1,4 @@
+1 3 3 1
+
+hw, hh, columns, rows
+hw and hh can be float values

ConvoluteLoop.C

@@ -1,7 +1,9 @@
 #include <TF1.h>
 #include <TMath.h>
 #include <TCanvas.h>
+#include <TStopwatch.h>
 #include <TPad.h>
+#include <TFile.h>
 #include <TH1.h>
 #include "Math/Integrator.h"
 #include "RLAfuncs.cxx"
@@ -27,7 +29,7 @@ double Func_1(double *t, double *par) // The detector-response function
 double Func_2(double *t, double *par) // The charge dispersion function
 {
   double result;
-  result = Qfuncs::PadRect(t[0], par[0], par[1], par[2], par[3]);
+  result = Qfuncs::PadRect(t[0], par[0], par[1], par[2], par[3], par[4]);
   //result = Qfuncs::TestTri(t[0]);
   return result;
 }
@@ -47,10 +49,14 @@ Other than that, this is exactly the convolution code that
 Klaus has posted on the skipper physics link.
 */
 
-void ConvoluteLoop(double * maspar, double * relativepos)
+void ConvoluteLoop(double * maspar, double * relativepos, char name[], char filename[], char title[], bool newFile)
 {
   //unpacking the master parameter array. These are set in the parent class.
-  double w = 0.25;
+  TStopwatch clock;
+  clock.Start();
+  double R = maspar[5]; //75 micrometers in PowerPoint
+  if(newFile==true){TFile f(filename,"RECREATE");}
+  TFile f(filename,"UPDATE");
   int n = int(maspar[0]+.1);
   double tR = maspar[1];
   double sL = maspar[2];
@@ -66,14 +72,12 @@ void ConvoluteLoop(double * maspar, double * relativepos)
   fQ = new TF1("fQ", Func_2, tmin, tmax, 5);// 5 params: xl, xh, yl, yh, w
   // To be adjusted
   fRLA->SetParameters(tR, sL, tr, tf);//(double t_R, sigmaL, t_r, t_f)=(50.,13.,40.,2000.) all in ns
-  fQ->SetParameters(xl, xh, yl, yh, w);//(xl, xh, yl, yh)=(-1, 1, -3, 3) all in mm
+  fQ->SetParameters(xl, xh, yl, yh, R);//(xl, xh, yl, yh)=(-1, 1, -3, 3) all in mm
   fRLA->SetNpx(n);
   fQ->SetNpx(n);
   double tcplus = 0;
   double dt = (tmax-tmin)/double(n);
-  char name[40];
-  sprintf(name,"(%g,%g) to (%g,%g))",xl,yl,xh,yh);
-  TH1D *cumConv = new TH1D("Signal",name, n, tmin, tmax);
+  TH1D *cumConv = new TH1D(name,title, n, tmin, tmax);
   cumConv->Reset();
   TH1D *cumInt = new TH1D("cumInt", "", n, tmin, tmax);
   cumInt->Reset();
@@ -101,10 +105,15 @@ void ConvoluteLoop(double * maspar, double * relativepos)
       }
   }
 
-      cumConv->SetBinContent(iT, -maxInt);//cumInt->GetBinContent(n));
+      cumConv->SetBinContent(iT, maxInt);//cumInt->GetBinContent(n));
     }
-  cumConv->SetBinContent(n, 1.0E-12);//cumInt->GetBinContent(n));
   cumConv->GetXaxis()->SetTitle("Time in ns");
-  cumConv->Draw();
-  cumConv->GetYaxis()->SetRangeUser(-3500,0);
+  cumConv->Write();
+  delete cumConv;
+  delete cumInt;
+  cumConv = NULL;
+  cumInt = NULL;
+  f.Close();
+  clock.Stop();
+  cout << " Overall time: "; clock.Print(); cout << endl;
 }

Input.dat

@@ -0,0 +1,5 @@
+7 3 .5 1 11 4 1000 -200 4
+
+X-start, Y-start, x-stepsize, y-stepsize, steps-in-x, steps-in-y, R-start, R-step, steps-in-R
+
+Number of steps include the 0 step! 

Qfuncs.cxx

@@ -17,15 +17,16 @@
 
 #include "Qfuncs.h"
 #include "TMath.h"
+#include "iostream"
 
 using namespace TMath;
 
 //______________________________________________________________________________
-double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, double w)
+double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, double R)
 {
   // Calculate dispersion of charge across 2-D resistive anode 
-  double N = 8000;// in q_e, typical COMPASS GEM value
-  double R = 2*530000;// Ohms per square
+  double N = 10000;// in q_e, 8000 is typical COMPASS GEM value
+  double w = .250;//530000;// Ohms per square
   double C = 0.21E-12;// 0.21pF/mm^2
   //double w = 250E-3;// in mm
   /*
@@ -35,14 +36,15 @@ double Qfuncs::PadRect(double tt, double xl, double xh, double yl, double yh, do
   double yl = par[2];//1;// in mm
   double yh = par[3];//7;// in mm
   */
-
   double Qpad;// = 1/Sqrt(4*t/(R*C)+2*Power(w,2));
+  double ttn=-999;
   if(tt>0)
     {
-      double S = Sqrt(4*(tt*Power(double(10),-9))/(R*C)+2*Power(w,2));
+      ttn=tt/double(1000000000);
+      double S = Sqrt(4*ttn/(R*C)+2*Power(w,2));
       double diffxErf = Erf(xh/S)-Erf(xl/S);
       double diffyErf = Erf(yh/S)-Erf(yl/S);      
-      Qpad = N/4 * ( diffxErf*diffyErf);
+      Qpad = N/double(4) * ( diffxErf*diffyErf );
     }
   else
     {

Qfuncs.h

@@ -33,7 +33,7 @@
 */
 
 namespace Qfuncs {
-  double PadRect(double tt=0, double xl=1, double xh=1, double yl=1, double yh=1, double w = 1);
+  double PadRect(double tt=0, double xl=1, double xh=1, double yl=1, double yh=1, double R = 1);
   double TestBox(double tt=0);
   double TestTri(double tt=0);
 }

aragorn.C

@@ -1,41 +1,45 @@
 #include "aragorn.h"
 #include "quink.C"
 #include <iostream>
-#include "TCanvas.h"
+#include <TCanvas.h>
+#include <TStopwatch.h>
+#include <TSystem.h>
+#include <fstream>
+#include <TFile.h>
 
 using namespace std;
 
-double matpar [5] = {30, 60, 13, 40, 2000}; //Parameters of the materiasl and the detector
+double matpar [6] = {40, 60, 13, 40, 2000, 0}; //Parameters of the materials and the detector
+TDatime now;
+
+float aragorn::hw;
+float aragorn::hh;
+int aragorn::rows;
+int aragorn::columns;
 
 aragorn *aragorn::__instance=0;
 
 aragorn::aragorn() {cout << "Their pace is quickening." << endl;}   //Just to confirm it runs correctly
 
-/* NONE OF THIS IS CORRECT ANYMORE, BUT I WILL FIX IT SOON
-The following is an example of how the quink array would look, with the quinks numbered as if they have already
-correctly been associated with each other for x and y equal to 4:
-                                              On top is electrode 0
-                                             ------------------------- 
-                                           | 42  43 44  45  46  47  48 |
-                                           | 35  36 37  38  39  40  41 |                   
-                                           | 28  29 30  31  32  33  34 | On the right is electrode 1
-                On the left is electrode 3 | 21  22 23  24  25  26  27 |        
-                                           | 14  15 16  17  18  19  20 |
-                                           | 7   8   9  10  11  12  13 |
-                                           | 0   1   2   3   4   5   6 |
-                                             -------------------------
-                                             On bottom is electrode 2
-*/
 void aragorn::Construct()
 {
+  ifstream in;
+  in.open("Construct.dat");
+  int c,r;
+  double a,b;
+  in>>a>>b>>c>>r;
+  Sethw(a);
+  Sethh(b);
+  Setcol(c);
+  Setrow(r);
   if(TheQuinks.size()>0){cout<<"This setup has already been constructed.";}
   else
   {
-   for(int j=hh;j<2*hh*rows;j=j+2*hh)                         
+   for(int j=1;j<=2*rows;j+=2)                         
    {                           //Creating all quinks, one row at at time and then moving down a row.
-     for(int i=hw;i<2*hw*columns;i=i+2*hw)
+     for(int i=1;i<=2*columns;i+=2)
      {
-       quink * newest =(new quink(double(i-hw),double(i+hw),double(j-hh),double(j+hh)));
+       quink * newest =(new quink(double(i*hw-hw),double(i*hw+hw),double(j*hh-hh),double(j*hh+hh)));
        TheQuinks.push_back(newest);
      }
    }
@@ -51,30 +55,170 @@ void aragorn::Report()
 
 void aragorn::SetChargeCloud(double a,double b) 
 {
-  //eventually this function will hopefully
+  // eventually this function will hopefully
   // take some x and y position and turn it
   // into a circle of charge with a dense
   // center and diffuse radial charge.
   // for today however, it does nothing.
 }
 
-void aragorn::Event()
+// This is just a function to make sure that the setup was constructed correctly
+void aragorn::CheckPad(int n)
+{
+  if(n>TheQuinks.size()-1){cout<<"This is not a valid pad number.";}
+  else if (n<0){cout<<"This is not a valid pad number";}
+  else{TheQuinks[n]->Report();}
+}
+
+void aragorn::Observe()
 {
-  double x1,y1; //x and y coordinate of charge cloud 
-
-  cout<<"Where should the charge cloud hit?"<<endl;
-  cin>>x1>>y1;
-  TCanvas* c1 = new TCanvas("c1", "Readout", 800, 600);
-  c1->Divide(columns,rows);
-  for(int i =1;i<=TheQuinks.size();i++)
+  int n = 7;
+  vector < vector<double> > values(n+1);
+  int step = 40/(n); //tmax divided by the number of time samples
+  TH1F * h1 = new TH1F();
+  char buffer [50];
+  TCanvas * c1 = new TCanvas("c1","Progression of Charge");
+  //TFile f("C:\\Users\\Overlord\\Desktop\\Programs\\Research\\ChargeDiv\\PRF\\PRFData3.root");
+  TFile f("C:\\root_v5.34.25\\macros\\SmallPads.root");
+  for (int i = 0; i <= n; ++i)
+  {
+    for (int k = 1; k < TheQuinks.size()+1; ++k)
+    {
+      sprintf(buffer,"Pad_li_%i",k);
+      h1 = (TH1F*)f.Get(buffer);
+      values[i].push_back(h1->GetBinContent(i*step+1)); //Value at a time i(index)*n(step size)
+      delete(h1);
+    }
+  }
+  c1->Divide(4,2); //The product of these need to be the n value +1;
+  double max;
+  double myval;
+  for (int i = 0; i <= n; ++i)
   {
-    c1->cd(i);
-    TheQuinks[i-1]->Display(matpar, x1, y1);
-    cout<<i<<" Pads have been calculated for.\n";
+    max = abs(values[i][12]);
+    c1->cd(i+1);
+    gPad->Divide(columns,rows,0,0);
+    for (int j = 0; j < values[i].size(); ++j)
+    {
+      gPad->cd(j+1);
+      myval = abs(values[i][j]);
+      if(myval>.95*max){gPad->SetFillColor(kRed);}
+      else if((myval<.95*max)&&(myval>=.5*max)){gPad->SetFillColor(kRed+1);}
+      else if((myval<.5*max)&&(myval>=.3*max)){gPad->SetFillColor(kRed+2);}
+      else if((myval<.3*max)&&(myval>=.1*max)){gPad->SetFillColor(kRed+3);}
+      else if((myval<.1*max)&&(myval>=.01*max)){gPad->SetFillColor(kRed+4);}
+      else {gPad->SetFillColor(kBlack);}
+      gPad->Draw();
+      c1->cd(i+1);
+    }
   }
 }
-// This is just a function to make sure that the setup was constructed correctly
-void aragorn::CheckPad(int n)
+
+void aragorn::SingleEvent()
+{
+  ifstream in;
+  in.open("Input.dat");
+  in.clear();
+  double x1,y1;
+  in>>x1>>y1;
+  in.close();
+  char name[40];
+  bool newFile = true;
+  char filename[40];
+  int timenow=now.Convert();
+  sprintf(filename,"Single_%i.root",timenow);
+  for (int j = 0; j<TheQuinks.size(); ++j)
+  {
+    if(j==1){newFile=false;}
+    sprintf(name,"Pad_li_%i",j+1);
+    TheQuinks[j]->ToFile(matpar, x1, y1, name, filename, newFile);
+    cout<<"Pad "<<j+1<<" Calculation Complete."<<endl;
+  }
+}
+
+void aragorn::MakePRF()
+{
+  ifstream in;
+  in.open("Input.dat");
+  double x1,y1,stepx,stepy,xtemp,ytemp;
+  int xev,yev;
+  in>>x1>>y1>>stepx>>stepy>>xev>>yev;
+  in.close();
+  char name[40];
+  char filename[40];
+  bool newFile = true;
+  int timenow=now.Convert();
+  sprintf(filename,"PRF_%i.root",timenow);
+  TStopwatch loopclock;
+  loopclock.Start();
+  for(int j=0;j<yev;++j)
+  {
+    ytemp=y1+j*stepy;
+    for(int i=0;i<xev;++i)
+    {
+      xtemp=x1+i*stepx;
+      for (int k = 0; k<TheQuinks.size(); ++k)
+      {
+        if(k==1){newFile=false;}
+        sprintf(name,"PRF_li_%i_lx_%i_ly_%i",k+1,i,j);
+        TheQuinks[k]->ToFile(matpar, xtemp, ytemp, name, filename, newFile);
+      }
+      cout<<i+1<<" x-event done."<<endl;
+    }
+  }
+  loopclock.Stop();
+  cout << endl; loopclock.Print(); cout << endl;
+}
+
+void aragorn::VaryResistance()
 {
-  TheQuinks[n]->Report();
+  ifstream in;
+  in.open("Input.dat");
+  double x1,y1,stepx,stepy,xtemp,ytemp,r1,stepr;
+  int xev,yev,rev;
+  in>>x1>>y1>>stepx>>stepy>>xev>>yev>>r1>>stepr>>rev;
+  in.close();
+  cout<<"R start "<<r1<<" step size "<<stepr<<" number of events "<<rev<<endl;
+  cin.get();
+  char name[40];
+  char filename[40];
+  bool newFile = true;
+  int timenow=now.Convert();
+  sprintf(filename,"ResistorsVariation_%i.root",timenow);
+  TStopwatch loopclock;
+  loopclock.Start();
+for (int r = 0; r <rev; ++r)
+{
+  matpar[5] = 1000*(r1+r*stepr);
+  cout<<" Resistance for this run is "<<matpar[5]<<endl;
+  for(int j=0;j<yev;++j)
+  {
+    ytemp=y1+j*stepy;
+    for(int i=0;i<xev;++i)
+    {
+      xtemp=x1+i*stepx;
+      for (int k = 0; k<TheQuinks.size(); ++k)
+      {
+        if(k==1){newFile=false;}
+        sprintf(name,"PRF_li_%i_lx_%i_ly_%i_lr_%i",k+1,i,j,matpar[5]/100000);
+        TheQuinks[k]->ToFile(matpar, xtemp, ytemp, name, filename, newFile);
+      }
+      cout<<i+1<<" x-event done."<<endl;
+    }
+  }
+}
+  loopclock.Stop();
+  cout << endl; loopclock.Print(); cout << endl;
+}
+
+void aragorn::Examine(vector< vector<double> > x)
+{
+  for(int i=0;i<x.size();i++)
+  {
+    for(int j=0;j<x[i].size();j++)
+    {
+      cout<<x[i][j]<<" ";
+    }
+    cout<<endl;
+  }
 }