TreeTest.C 5.89 KB
Newer Older
1
void TreeEnergyWrite(){
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
2
3


4
5
6
7
  int Ecol1;
  double Ecol2,Ecol3,Ecol4;
  string line_E;
   
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
8
  
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
  TFile* Efile = new TFile("Energy_table.root","recreate");
  TTree* ETree = new TTree("ETree","Energy tree");

  ETree->Branch("Ecol1",&Ecol1,"Ecol1/I"); // ECol1 = Energy values = {200,300.....14000}
  ETree->Branch("Ecol4",&Ecol4,"Ecol4/D"); // Voltage values 

  std::ifstream Energy_file("Energy_table.dat");
  if(Energy_file.is_open()){
    while(getline(Energy_file,line_E)){
      istringstream ssE(line_E);
      ssE >> Ecol1 >> Ecol2 >> Ecol3 >> Ecol4;
      ETree->Fill();
    }
  }

  ETree->Write();
  Efile->Close();
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
26
  
27
28
}

29
30
31
void TreeTPWrite(){
 
  int energy, col, row, pc;
32
  double voltage,tot, tot_var, toa, toa_var, time_ref,toa_cut,mean,RMS;
33
34
35
36
  string line_TP;

  TFile* Efile = new TFile("Energy_table.root");
  TTree* ETree = (TTree*)Efile->Get("ETree");
37
38
  ETree->SetBranchAddress("Ecol1",&energy);
  ETree->SetBranchAddress("Ecol4",&voltage);
39

Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
40
  
41
42
43
44
45
  TFile* TPfile = new TFile("TP.root","update");  
  TTree* TPtree = new TTree("TPtree"," TP tree");
  TPtree->Branch("col",&col,"col/I");
  TPtree->Branch("row",&row,"row/I");
  TPtree->Branch("pc",&pc,"pc/I");
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
46
47
48
49
  TPtree->Branch("tot",&tot,"tot/D");
  TPtree->Branch("tot_var",&tot_var,"tot_var/D");
  TPtree->Branch("toa",&toa,"toa/D");
  TPtree->Branch("toa_var",&toa_var,"toa_var/D");
50
51
  TPtree->Branch("time_ref",&time_ref,"time_ref/D");
  auto voltageBranch = TPtree->Branch("voltage",&voltage,"voltage/D");// make a new branch and connect it to the same variable as the branch address of the energy tree
52
  auto Branch_toacut = TPtree->Branch("toa_cut",&toa_cut,"toa_cut/D");
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
53

54
55
56
57
58
59
  for(int i = 0; i < ETree->GetEntries(); i++){
    ETree->GetEntry(i);
    // std:: cout << energy << " " << voltage << std::endl;
    std::string TPfilename = "tp_dat_" + std::to_string(energy) + ".dat";
    //std:: cout << " The TPfilenames are : " <<  TPfilename.c_str() <<std::endl;
    std::ifstream TPdatafile(TPfilename.c_str());
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
60

61
62
    voltageBranch->Fill();
    
63
64
65
66
    if(TPdatafile.is_open()){

      std::cout << " The data file " << TPfilename.c_str() << " is opened ..." <<std::endl;

67
      
68
69
70
      // reading through individial tp_dat_ files

      while(getline(TPdatafile,line_TP)){
71
72
73
74
75
76
	  if(line_TP.rfind("#", 0) == 0) {
	    continue;
	  }

	  istringstream ssTP(line_TP);
	  ssTP >> col >> row >> pc >> tot >> tot_var >> toa >> toa_var >> time_ref;
77
78
79
	  TPtree->Fill();
      } //While loop getline
    }//if statement TPdatafile open
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

    TPtree->Draw("toa",Form("voltage==(%f)",voltage),"goff");
    mean = TMath::Mean(TPtree->GetSelectedRows(), TPtree->GetV1());
    RMS = TMath::RMS(TPtree->GetSelectedRows(), TPtree->GetV1());
    // std::cout <<" mean " << mean << " RMS " << RMS <<std::endl;

    /* for (int isigma=0; isigma < 15; isigma++){
      std::cout << " cut at " << isigma << " times the RMS leaves : " 
      << 100.0*TPtree->GetEntries(Form("toa<(%f) && voltage==(%f)",mean+isigma*RMS,voltage))/TPtree->GetEntries(Form("voltage==(%f)",voltage)) 
      << " percentage of entries " << std::endl;
      }*/

    toa_cut = mean + 6*RMS;
    Branch_toacut->Fill();
    
95
96
97
98
  }//For loop GetEntries

  //Save only the new version of the tree 

99
  TPtree->Write("",TObject::kOverwrite);
100
  
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  TPfile->Close();
}


void timewalkfit(){


  // gROOT->Reset();
  
  int col, row, pc ;
  double tot, tot_var, toa, toa_var;

  TFile* TPfile = new TFile("TP.root");
  TTree* TPtree = (TTree*)TPfile->Get("TPtree");
  TPtree->SetBranchAddress("col",&col);
  TPtree->SetBranchAddress("row",&row);
  TPtree->SetBranchAddress("pc",&pc);
  TPtree->SetBranchAddress("tot",&tot);
  TPtree->SetBranchAddress("tot_var",&tot_var);
  TPtree->SetBranchAddress("toa",&toa);
  TPtree->SetBranchAddress("toa_var",&toa_var);


  for(int i = 0; i < TPtree->GetEntries(); i++){

    TPtree->GetEntry(i);
    std::cout << col << " " << row << std::endl; //Correct values of col and row printed

    auto N = TPtree->Draw("voltage:toa:sqrt(toa_var)",Form("col==(%i) && row==(%i) && toa>0 && toa<toa_cut", col, row),"goff");
    // cout << N <<std::endl;
    TGraphErrors *graph = new TGraphErrors(N,TPtree->GetV1(),TPtree->GetV2(),0,TPtree->GetV3());

    // std::cout << col << " " << row << std::endl; // Always printing 255 255
    
    auto fitmin = TMath::MinElement(graph->GetN(),graph->GetX());
    // std::cout << "fitmin" << fitmin << std::endl;
    
    auto fitmax = TMath::MaxElement(graph->GetN(),graph->GetX());
    // std::cout << "fitmax " << fitmax << std::endl;
    
    TF1 *fitfuncTOA = new TF1("fitfuncTOA","(([2]/(x-[0]))+ [1])",fitmin,fitmax);
    fitfuncTOA->SetParNames("asymtote","offset","curvature");
    
    int iymax=TMath::LocMax(graph->GetN(),graph->GetY());
    // cout << "iymax " << iymax << std::endl;
    auto ttoa = graph->GetPointX(iymax);
    // cout << "ttoa" << " " << ttoa << std::endl;;
    fitfuncTOA->SetParameter(0,ttoa);
    
    auto xmax = TMath::MaxElement(graph->GetN(),graph->GetX());
    auto dtoa = graph->Eval(xmax);
    // cout << "dtoa" << " " << dtoa <<std::endl;
    fitfuncTOA->SetParameter(1,dtoa);
    
    auto testx = 0.5*fitmax;
    auto ctoa_y = graph->Eval(testx);
    auto ctoa = (ctoa_y - dtoa)*(testx - ttoa);
    // cout << "ctoa" << " " << ctoa <<std::endl;
    fitfuncTOA->SetParameter(2,ctoa);
    
    // gStyle->SetOptFit(1111);
    // fitfuncTOA->SetLineColor(kRed);
    graph->Fit(fitfuncTOA,"R");
    /* graph->GetXaxis()->SetTitle("Volatge[mV]");
       graph->GetYaxis()->SetTitle ("time difference[ns]");
       graph->SetTitle(" ");
       graph->SetMarkerStyle(7);*/
    // graph->Draw("AP");
    
    /* std::ofstream fitpars;
       fitpars.open("fit_parameters.txt",std::ofstream::app);
       fitpars << col << "\t" << row << "\t"  << fitfuncTOA->GetParameter(2)
       << setw(10) << fitfuncTOA->GetParameter(0) << setw(10) << fitfuncTOA->GetParameter(1) << std::endl;*/ //NO chi2/ndf for the time-being!
   
  }

177
178
  TPfile->Close();
}
Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
179

180
void TreeTest(){ 
181
182

  TreeEnergyWrite();
183
  TreeTPWrite();
184
185
  timewalkfit();

Keerthi Nakkalil's avatar
Keerthi Nakkalil committed
186
}