full_toa.C 3.62 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
std::vector<double> energy;
std::vector<double> TOA;
std::vector<double> TOA_ERR;
      
void timewalkfit(int pixel_col,int pixel_row,std::vector<double> &voltage, std::vector<double> &time,std::vector<double> &timeerr){

  TGraphErrors* graph = new TGraphErrors(time.size(), &(voltage[0]), &(time[0]), 0, &(timeerr[0]));

  TF1 *fitfuncTOA = new TF1("fitfuncTOA","(([0]/(x-[1]))+ [2])",25.0,700.0);
  fitfuncTOA->SetParNames("curvature","asymtote","offset");

  auto ymax = TMath::MaxElement(graph->GetN(),graph->GetY());
  /*int iymax=TMath::LocMax(graph->GetN(),graph->GetY()); -> Fitting did not work here properlyy!!
    auto ttoa = graph->GetPointX(iymax);*/ 
  auto ttoa = graph->Eval(ymax);
  // cout << "ttoa" << " " << ttoa << std::endl;;
  fitfuncTOA->SetParameter(1,ttoa);

  auto xmax = TMath::MaxElement(graph->GetN(),graph->GetX());
  auto dtoa = graph->Eval(xmax);
  // cout << "dtoa" << " " << dtoa <<std::endl;
  fitfuncTOA->SetParameter(2,dtoa);

  auto ctoa_y = graph->Eval(400);
  auto ctoa = (ctoa_y - dtoa)*(400 - ttoa);
  // cout << "ctoa" << " " << ctoa <<std::endl;
  fitfuncTOA->SetParameter(0,ctoa);

  // fitfuncTOA->SetLineColor(kRed);
  graph->Fit(fitfuncTOA,"R");
  // graph->Draw("AP");


  std::ofstream fitpars;
  fitpars.open("fit_parameters_new.txt",std::ofstream::app);
  fitpars << pixel_col << "\t" << pixel_row << "\t"  << fitfuncTOA->GetParameter(0)
	  << setw(10) << fitfuncTOA->GetParameter(1) << setw(10) << fitfuncTOA->GetParameter(2)
	  << setw(10) << (fitfuncTOA->GetChisquare()/fitfuncTOA->GetNDF()) << std::endl;

}  

 
void createfiles(int c, int r){

  std::ifstream Energy_file("Energy_table.dat");
  std::vector<double>Energy_e = {200,300,400,500,600,700,800,900,1000,1100,1200,
                                1300,1400,1500,1700,2200,2700,3200,4000,5000,6000,
				7000,8000,9000,10000,11000,12000,13000,14000};
  if (Energy_file.is_open()){
    std::cout << "Energy_table.dat is opened for pixel" << std::to_string(r) << "_" << std::to_string(c) << std::endl;

    std::string line_1,line_2;
    int Ecol1,Ecol2,Ecol3, TPcol1,TPcol2,TPcol3;
    int ctr_E = 0, ctr_TOA = 0;
    double Ecol4=0, tot,tot_var,toa=0,toa_var =0,toa_ref ;
   
    
    while(getline(Energy_file,line_1)){
      istringstream ss1(line_1);
      ss1 >> Ecol1 >> Ecol2 >> Ecol3 >> Ecol4;
      energy.push_back(Ecol4);
      // std::cout << "Elements " << Ecol4[ctr_E] << " ctr " << ctr_E << std::endl;
      ctr_E++;
    }
    
    for(int i= 0; i < Energy_e.size(); i++){
      stringstream Sfile;
      Sfile << Energy_e.at(i) << ".dat";
      // std::cout << "Files : " << Sfile.str() << std::endl;
      std::vector<std::string>File_names = {std::string("tp_dat_" + Sfile.str())};
      
      for (int j = 0; j < File_names.size(); j++ ) {

	std::ifstream Data_file(File_names.at(j));

	if(Data_file.is_open()){
	  //  std::cout << "The data file " << File_names.at(j) << " is opened ..." << std::endl;
	  while(getline(Data_file,line_2)){
	    // cout << " looking at each line ... " <<std::endl;
	    if(line_2.rfind("#", 0) == 0) {
	      continue;
	    }

	    istringstream ss2(line_2);
	    ss2 >> TPcol1 >> TPcol2 >> TPcol3 >> tot >> tot_var >> toa >> toa_var;
	    if (TPcol1==c){
	      if (TPcol2==r){
		
		TOA.push_back(toa);
		auto toa_err = sqrt(toa_var);
		TOA_ERR.push_back(toa_err);
		ctr_TOA++;
	      }
	    }
	  }
	}
	else {
	  std::cout << "Unable to open " <<File_names.at(j) << " file ..." << std::endl;
	}
      }
    }
  }

  timewalkfit(c,r,energy,TOA,TOA_ERR);
}


void full_toa(){
  int col,row;

  for(row =0; row <256; row ++){
    for(col =0; col <256; col++){
      createfiles(col,row);
    }
  }
  std::cout << " Done ..." <<std ::endl;
}