prova.cxx

#ifndef prova_cxx
#define prova_cxx

#include "prova.h"
#include <TH2.h>
#include <TH1.h>
#include <TLine.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <iostream>
#include <vector>
#include "TROOT.h"
#include <TLegend.h>
#include <fstream>
#include <algorithm>
#include <map>
#include <TGraph.h>
#include <TDirectory.h>
#include <math.h>
#include <TMultiGraph.h>
#include <TObject.h>
#include <TTree.h>
#include <string.h>
#include <TString.h>
#include <TProfile.h>
#include "cluster.h"

using namespace std;

void prova::Loop()
{
  // In a ROOT session, you can do:
  //   root> .L prova.C
  //   root> prova t
  //   root> t.GetEntry(12); // Fill t data members with entry number 12
  //   root> t.Show();       // Show values of entry 12
  //   root> t.Show(16);     // Read and show values of entry 16
  //   root> t.Loop();       // Loop on all entries

  // This is the loop skeleton where:
  //   jentry is the global entry number in the chain
  //   ientry is the entry number in the current Tree
  // Note that the argument to GetEntry must be:
  //   jentry for TChain::GetEntry
  //   ientry for TTree::GetEntry and TBranch::GetEntry

  // To read only selected branches, Insert statements like:
  // METHOD1:
  //   fChain->SetBranchStatus("*",0);  // disable all branches
  //   fChain->SetBranchStatus("branchname",1);  // activate branchname
  // METHOD2: replace line
  //   fChain->GetEntry(jentry);       //read all branches
  //by  b_branchname->GetEntry(ientry); //read only this branch
  //   gROOT−>Reset();
  //   gROOT−>SetStyle("Plain");
  //   gStyle − > SetOptStat ( 1 1 1 1 1 1 ) ;
  //   gStyle − > SetOptFit ( 1 1 1 1 ) ;
  //   gStyle − > SetPalette ( 1 ) ;
  //   gStyle − > SetOptTitle ( 0 ) ;
  if (fChain == nullptr)
    return;

  gErrorIgnoreLevel = kWarning;

  Long64_t nentries = fChain->GetEntriesFast(); 
  Long64_t nbytes = 0, nb = 0;
  //prendo il nome del run e lo uso per chiamarci i file su cui salvo gli output
  TString namefile =  ((TChain*)(fChain))->GetFile()->GetName();
  TString namef = namefile;
  namef.ReplaceAll("data/","");
  namef.Remove(namef.Index(".root"));
  TFile* f = new TFile(("/home/leonardo/Scrivania/Roba/TBReco/root_output/class_test_" + namef + ".root"), "RECREATE");
  _outdirname_lay = "LAYERS";
  f->mkdir(_outdirname_lay);

  //QUI CI VANNO LE DEFINIZIONI DEGLI HISTO PRIMA DI LOOPPARE SUGLI EVENTI
  //HISTO CAMERE
  vector<TH1F*> h_strip_chamb;
  vector<TH1F*> h_max_q_chamb;
  vector<TProfile*> profile_strip_max_q;
  vector<TH2F*> h2_strip_max_q_chamb;
  vector<TH1F*> h_number_strip_chamb;
  vector<TH1F*> h_number_cluster_chamb;
  vector<TH1F*> h_centroid_chamb;
  vector<TH1F*> h_sum_q_cluster_chamb;
  vector<TH2F*> h2_centroid_sum_q;
  //vector<TGraph*> centroid_reco_vs_real;
  TH1* evt_vs_layer = new TH1F("evt_vs_layer","Event vs layer turned on; Evt; Number of Layer",4,0.5,4.5);
  for(int j = 0; j<4; j++)
  {
    h_strip_chamb.push_back(new TH1F(("h_strip_chamb_"+ to_string(j+1)).c_str(),("Histo strip chamb "+to_string(j+1)+";Strip;Entries").c_str(),530,500,1030));
    h_max_q_chamb.push_back(new TH1F(("h_max_q_chamb_"+ to_string(j+1)).c_str(),("Distribution charge max chamb "+to_string(j+1)+";Max_q;Entries").c_str(),150,0.0,2500));
    profile_strip_max_q.push_back(new TProfile(("profile_strip_max_q"+ to_string(j+1)).c_str(),("profile_strip_max_q"+to_string(j+1)+":Strip;Profile Max_q").c_str(),530,500,1030,0.0,2500));
    h2_strip_max_q_chamb.push_back(new TH2F(("h2_strip_max_q_chamb_"+ to_string(j+1)).c_str(),("Histo strip vs max_q chamb "+to_string(j+1)+";Strip;Max_q").c_str(),530,500,1030,150,0.0,2500));
    h_number_strip_chamb.push_back(new TH1F(("h_number_strip_chamb_"+ to_string(j+1)).c_str(),("Histo molteplicita' strip per evento chamb "+to_string(j+1)+";Number of strip;Event").c_str(),200,0.0,200));
    h_number_cluster_chamb.push_back(new TH1F(("h_number_cluster_chamb_"+ to_string(j+1)).c_str(),("Histo molteplicita' cluster per evento chamb "+to_string(j+1)+";Number of cluster;Event").c_str(),50,0.0,50));
    h_centroid_chamb.push_back(new TH1F(("h_centroid_chamb_"+ to_string(j+1)).c_str(),("Histo centroid chamb "+to_string(j+1)+";Centroid [mm];Entries").c_str(),200,0.0,2200));
    h_sum_q_cluster_chamb.push_back(new TH1F(("h_sum_q_cluster_chamb_"+ to_string(j+1)).c_str(),("Distribution sum of charge chamb "+to_string(j+1)+";Sum_q;Event").c_str(),500,0.0,10000));
    h2_centroid_sum_q.push_back(new TH2F(("h2_centroid_sum_q"+ to_string(j+1)).c_str(),("Centroid vs sum of charge chamb "+to_string(j+1)+";Centroid [mm];Sum_q").c_str(),200, 0.0, 2200, 500,0.0,10000));
    h2_extrapolate_pos.push_back(new TH2F(("h2_extrapolate_pos"+ to_string(j+1)).c_str(),("Extrapolate cluster position vs real position"+to_string(j+1)+";Extrapolate[mm];Real[mm]").c_str(),200,0.0,2200,200,0.0,2200));
    h2_extrapolate_pos_fit.push_back(new TH2F(("h2_extrapolate_pos_fit"+ to_string(j+1)).c_str(),("Extrapolate cluster position vs real position by fit"+to_string(j+1)+";Extrapolate[mm];Real[mm]").c_str(),200,0.0,2200,200,0.0,2200));
    h_residues_pos_fit.push_back(new TH1F(("h_residues_pos_hist"+ to_string(j+1)).c_str(),("Residues of extrapolate cluster position vs real position by fit"+to_string(j+1)+";Residues[mm];Entries").c_str(),200,-100,100));
    h_residues_pos.push_back(new TH1F(("h_residues_pos"+ to_string(j+1)).c_str(),("Residues of extrapolate cluster position vs real position"+to_string(j+1)+";Residues[mm];Entries").c_str(),200,-100,100));
  }
  gROOT->SetStyle("Plain");
  double percentual = 0;
  int just_one = 0;
  for(Long64_t jentry=0; jentry<nentries;jentry++)
  {
    percentual = (double)jentry/((double)nentries);
    if((fmod(percentual*100,5) < 0.9) && (just_one == 0))
    { 
      cout<<"["<<(int)(percentual*100)<<"%]" <<endl;
      just_one ++;
    }
    else
      if(fmod(percentual*100,5) >= 0.9)
        just_one = 0;
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0)
      break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    n_event++;
    Inizializ();
    Int_t laytot=0;
    //******MAPPE
    Mapping(ch1, _ch_p1);
    Mapping(ch2, _ch_p2);
    mapINFNcards();
    Mapping(ch3, _ch_p3);
    Mapping(ch4, _ch_p4);
    //HISTO DI TUTTA LA CAMERA
    Histo_Layer(h_strip_chamb,h_max_q_chamb,profile_strip_max_q,h2_strip_max_q_chamb,h_number_strip_chamb,ch1,0);
    Histo_Layer(h_strip_chamb,h_max_q_chamb,profile_strip_max_q,h2_strip_max_q_chamb,h_number_strip_chamb,ch2,1);
    Histo_Layer(h_strip_chamb,h_max_q_chamb,profile_strip_max_q,h2_strip_max_q_chamb,h_number_strip_chamb,ch3,2);
    Histo_Layer(h_strip_chamb,h_max_q_chamb,profile_strip_max_q,h2_strip_max_q_chamb,h_number_strip_chamb,ch4,3);
    //******CLUSTER
    Cluster t1(ch1, _ch_p1);
    Cluster t2(ch2, _ch_p2);
    Cluster t3(ch3, _ch_p3);
    Cluster t4(ch4, _ch_p4);
    //******EVENT DISPLAY
    gStyle->SetPalette(1);  

    evt_vs_layer->Fill(laytot);
    h_number_cluster_chamb[0]->Fill(t1.GetN_Cluster());
    h_number_cluster_chamb[1]->Fill(t2.GetN_Cluster());
    h_number_cluster_chamb[2]->Fill(t3.GetN_Cluster());
    h_number_cluster_chamb[3]->Fill(t4.GetN_Cluster());

    if(jentry < 100)
    {
      cout<<"isGoodCluster dice che il cluster è ";
      if(t1.isGoodCluster(100,_ch_p1))
        cout<<"Buono"<<endl;
      else
        cout<<"Cattivo"<<endl;
    }

    vector<double> centroid_1 = t1.GetCentroid();
    vector<double> centroid_2 = t2.GetCentroid();
    vector<double> centroid_3 = t3.GetCentroid();
    vector<double> centroid_4 = t4.GetCentroid();
    vector<string> chamb_1 = t1.GetChamb();
    vector<string> chamb_2 = t2.GetChamb();
    vector<string> chamb_3 = t3.GetChamb();
    vector<string> chamb_4 = t4.GetChamb();
    vector<int>  q_sum_1 = t1.GetSum_q();
    vector<int>  q_sum_2 = t2.GetSum_q();
    vector<int>  q_sum_3 = t3.GetSum_q();
    vector<int>  q_sum_4 = t4.GetSum_q();

    for(unsigned int k = 0; k<chamb_1.size(); k++)
    {
      if(centroid_1.at(k)!=0)
      {
        h_centroid_chamb[0]->Fill(centroid_1.at(k));
        h_sum_q_cluster_chamb[0]->Fill(q_sum_1.at(k));
        h2_centroid_sum_q[0]->Fill(centroid_1.at(k),q_sum_1.at(k));
      }
    }
    for(unsigned int k = 0; k<chamb_2.size(); k++)
    {
      if(centroid_2.at(k)!=0)
      {
        h_centroid_chamb[1]->Fill(centroid_2.at(k));
        h_sum_q_cluster_chamb[1]->Fill(q_sum_2.at(k));
        h2_centroid_sum_q[1]->Fill(centroid_2.at(k),q_sum_2.at(k));
      }
    }
    for(unsigned int k = 0; k<chamb_3.size(); k++)
    {
      if(centroid_3.at(k)!=0)
      {
        h_centroid_chamb[2]->Fill(centroid_3.at(k));
        h_sum_q_cluster_chamb[2]->Fill(q_sum_3.at(k));
        h2_centroid_sum_q[2]->Fill(centroid_3.at(k),q_sum_3.at(k));
      }
    }
    for(unsigned int k = 0; k<chamb_4.size(); k++)
    {
      if(centroid_4.at(k)!=0)
      {
        h_centroid_chamb[3]->Fill(centroid_4.at(k));
        h_sum_q_cluster_chamb[3]->Fill(q_sum_4.at(k));
        h2_centroid_sum_q[3]->Fill(centroid_4.at(k),q_sum_4.at(k));
      }
    }
  }  
  for(UInt_t i = 0; i<h_strip_chamb.size(); i++)
  {
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_strip_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_max_q_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h2_strip_max_q_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(profile_strip_max_q.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_number_strip_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_number_cluster_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_centroid_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h_sum_q_cluster_chamb.at(i));
    f->GetDirectory(_outdirname_lay)->WriteTObject(h2_centroid_sum_q.at(i));
  }
  cout<<"HISTO SALVATI IN "<<"/home/leonardo/Scrivania/Roba/TBReco/root_output/class_test_" + namef + ".root"<<endl;

  h_strip_chamb.clear();
  h_max_q_chamb.clear();
  h2_strip_max_q_chamb.clear();
  h_number_strip_chamb.clear();
  h_number_cluster_chamb.clear();
  h_centroid_chamb.clear();
  h_sum_q_cluster_chamb.clear();
  h2_centroid_sum_q.clear();
}

void prova::Inizializ()
{
  n_clu_1 = 0;
  n_clu_2 = 0;
  n_clu_3 = 0;
  n_clu_4 = 0;
  for(int i = 0; i<20; i++)
    pcb_n_cluster[i] = 0;
  strip.clear();
  cluster_max_q.clear();
  chamb.clear();
  centroid.clear();
  q_sum.clear();
  n_cluster.clear();
  ch1.clear();
  ch2.clear();
  ch3.clear();
  ch4.clear();
}

void prova::Output_strip(ofstream &myfile)
{
  myfile.width(12);
  myfile<<"CHAMBER|";
  myfile.width(12);
  myfile<<"CLUSTER STRIP|";
  myfile.width(11);
  myfile<<"MAX_Q|";
  myfile.width(12);
  myfile<<"CENTROID|";
  myfile.width(12);
  myfile<<"SUM_Q|";
  myfile.width(13);
  myfile<<"N_CLUSTER|\n";
  for(UInt_t i = 0; i<chamb.size(); i++)
  {
    myfile.width(11);
    myfile<<chamb.at(i)<<"|";
    myfile.width(13);
    myfile<<strip.at(i)<<"|";
    myfile.width(10);
    myfile<<cluster_max_q.at(i)<<"|";
    myfile.width(11);
    myfile<<centroid.at(i)<<"|";
    myfile.width(11);
    myfile<<q_sum.at(i)<<"|";
    myfile.width(11);
    myfile<<n_cluster.at(i)<<"|\n";
  }
}

void prova::Output_ped(ofstream &myfile)
{
  myfile.width(12);
  myfile<<"mmChamber|";
  myfile.width(12);
  myfile<<"mmLayer|";
  myfile.width(11);
  myfile<<"mmReadout|";
  myfile.width(12);
  myfile<<"mmStrip|";
  myfile.width(12);
  myfile<<"ped_mean|";
  myfile.width(13);
  myfile<<"ped_stdev|\n";
  for(UInt_t i = 0; i<ped_chamb->size(); i++)
  {
    myfile.width(11);
    myfile<<ped_chamb->at(i)<<"|";
    myfile.width(11);
    myfile<<ped_layer->at(i)<<"|";
    myfile.width(10);
    myfile<<ped_readout->at(i)<<"|";
    myfile.width(11);
    myfile<<ped_Strip->at(i)<<"|";
    myfile.width(11);
    myfile<<ped_mean->at(i)<<"|";
    myfile.width(12);
    myfile<<ped_stdev->at(i)<<"|\n";
  }
}

void prova::Mapping(map<int,int> &map, string &name)
{
  for(UInt_t i = 0; i<mmChamber->size(); i++)
  {
    if(mmChamber->at(i) == name)
    {
      /*
      for(UInt_t j = 0; j<ped_chamb->size(); j++)
      {
        if((ped_chamb->at(j) == name) && (mmStrip->at(i) == ped_Strip->at(j)) && (max_q->at(i)/ped_stdev->at(j)>5))
        {
      */
      map.insert(pair<int,int>(mmStrip->at(i),max_q->at(i)));
        //}
      //}
    }
  }
}
void prova::mapINFNcards()
{
  vector<int> prima, seconda;
  ifstream file;
  file.open("/home/leonardo/Scrivania/Roba/TBReco/input_file/accoppiamento_strip.txt");
  if (file.is_open())
  {  }
  else
    cout<<"ERROR OPENING FILE"<<endl;
  int x,y,z;
  string a,b,c;
  while(file>> a >> x >> b >> y >> c >> z)
  {
    prima.push_back(x+512);
    seconda.push_back(y+512);
  }
  //  for(UInt_t i = 0; i<prima.size(); i++)
  //  {
  //    cout<<"prima "<<prima.at(i)<<"   seconda "<<seconda.at(i)<<endl;
  //  }
  for(UInt_t i = 0; i<prima.size(); i++)
  {
    for(const auto& it: ch2)
    {
      if(seconda.at(i) == it.first)
      {
        ch2.insert(pair<int,int>(prima.at(i), it.second));
      }
    }
  }
}

void prova::Histo_Layer(vector<TH1F*>& h_strip, vector<TH1F*>& h_max_q, vector<TProfile*>& p_strip_q, vector<TH2F*>& h2_strip_q, vector<TH1F*>& h_molt,  map<int,int>& map, int i)
{
  for(const auto& it: map)
  {
    h_strip[i]->Fill(it.first);
    h_max_q[i]->Fill(it.second);
    p_strip_q[i]->Fill(it.first,it.second);
    h2_strip_q[i]->Fill(it.first,it.second);
  }
  h_molt[i]->Fill(map.size());
}

void prova::Efficiency_1()
{
  double par[2];
  double x_axis[5] = { 1., 12.562, 33.08, 38.86, 44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, super_eta = 0;
  if(n_clu_2 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && (DAQCount(_ch_p2)==true) && (DAQCount(_ch_p3)==true) && (DAQCount(_ch_p4)==true))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
    }
    super_eta /= (2*cos(angle));
    parameters(temp_centroid_2, super_eta, x_axis[1], x_axis[3], par);
    if(par[1]<=0.404 && par[1]>= -0.404)
    {
      den_eff[0] = den_eff[0] + 1;
      for(int i = 0; i<5; i++)
      {
        max_road[i] = par[0]+5+par[1]*(x_axis[i]);
        min_road[i] = par[0]-5+par[1]*(x_axis[i]);
      }
      for(UInt_t i = 0; i<centroid.size(); i++)
      {
        if((chamb.at(i) == "P1_M01") && (centroid.at(i) != 0))
          h2_extrapolate_pos[0]->Fill(par[0]+par[1],centroid.at(i));
        if((chamb.at(i) == "P1_M01") && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[0] && centroid.at(i)>=min_road[0]))
        {
          num_eff[0] = num_eff[0] + 1;
          h_residues_pos[0]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[0]));
          break;
        }
      }
    }
    else
    {
      return;
    }
  }
  else
  {
    return;
  }
}

void prova::Efficiency_2()
{
  double par[2];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, super_eta = 0;
  if(n_clu_1 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && (DAQCount(_ch_p1)==true) && (DAQCount(_ch_p3)==true) && (DAQCount(_ch_p4)==true))
  {
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
    }
    super_eta /= (2*cos(angle));
    parameters(temp_centroid_1, super_eta, x_axis[0], x_axis[3], par);
    if(par[1]<=0.404 && par[1] >= -0.404)
    {
      den_eff[1] = den_eff[1] + 1;
      for(int i = 0; i<5; i++)
      {
        max_road[i] = par[0]+5+par[1]*(x_axis[i]);
        min_road[i] = par[0]-5+par[1]*(x_axis[i]);
      }
      for(UInt_t i = 0; i<centroid.size(); i++)
      {
        if((chamb.at(i) == "P2_M01") && (centroid.at(i) != 0))
          h2_extrapolate_pos[1]->Fill(par[0]+par[1]*x_axis[1],centroid.at(i));
        if((chamb.at(i) == "P2_M01") && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[1] && centroid.at(i)>=min_road[1]))
        {
          num_eff[1] = num_eff[1] + 1;
          h_residues_pos[1]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[1]));
          break;
        }
      }
    }
    else
    {
      return;
    }
  }
  else
    return;
}

void prova::Efficiency_stereo()
{
  double par[2];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2;
  if(n_clu_1 == 1 && n_clu_2  == 1 && (DAQCount(_ch_p1)==true) && (DAQCount(_ch_p2)==true))
  { 
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
      {
        temp_centroid_1 = centroid.at(k);
      }
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
      {
        temp_centroid_2 = centroid.at(k);
      }
    }
    //double start_point = (temp_centroid_1+temp_centroid_2)/2;
    parameters(temp_centroid_1, temp_centroid_2, x_axis[0], x_axis[1], par);
    if(par[1]<=0.404 && par[1] >= -0.404)
    {
      den_eff[2] = den_eff[2] + 1;
      den_eff[3] = den_eff[3] + 1;
      for(int i = 0; i<5; i++)
      {
        //max_road[i] = par[0]+Stereo_Road(start_point)+par[1]*(x_axis[i]);
        //min_road[i] = par[0]-Stereo_Road(start_point)+par[1]*(x_axis[i]);
        max_road[i] = par[0]+5+par[1]*(x_axis[i]);
        min_road[i] = par[0]-5+par[1]*(x_axis[i]);
      }
      //cout<<(temp_centroid_1+temp_centroid_2)/2<<endl;
      //cout<<Stereo_Road(start_point)<<endl;
      for(UInt_t i = 0; i<centroid.size(); i++)
      {
        if((chamb.at(i) == "P3_M01") && (centroid.at(i) != 0))
          h2_extrapolate_pos[2]->Fill(par[0]+par[1]*x_axis[2],centroid.at(i));
        if((chamb.at(i) == "P3_M01") && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[2] && centroid.at(i)>=min_road[2]))
        {
          num_eff[2] = num_eff[2] + 1;
          h_residues_pos[2]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[2]));
          break;
        }
      }
      for(UInt_t i = 0; i<centroid.size(); i++)
      {
        if((chamb.at(i) == "P4_M01") && (centroid.at(i) != 0))
          h2_extrapolate_pos[3]->Fill(par[0]+par[1]*x_axis[4],centroid.at(i));
        if((chamb.at(i) == "P4_M01") && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[4] && centroid.at(i)>=min_road[4]))
        {
          num_eff[3] = num_eff[3] + 1.;
          h_residues_pos[3]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[4]));
          break;
        }
      }
    }
    else
    {
      return;
    }
  }
  else
    return;
}

void prova::parameters(double &y1, double &y2, double &x1, double &x2, double par[2])
{
  par[0] = (x2*y1-x1*y2)/(x2-x1);
  par[1] = (y2-y1)/(x2-x1);
}

void prova::Efficiency_Fit_1()
{
  double par[3];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, temp_centroid_3, temp_centroid_4;
  if(n_clu_2 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && DAQCount(_ch_p2) && DAQCount(_ch_p3) && DAQCount(_ch_p4))
  { 
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        temp_centroid_3 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        temp_centroid_4 = centroid.at(k);
    }
  }
  else
    return;
  FitParameters(temp_centroid_2, temp_centroid_3, temp_centroid_4, x_axis[1], x_axis[2], x_axis[4], par);
  /*
  if(n_event <= 100000)
  {
    TCanvas *c = new TCanvas("c",("event_display " + to_string(par[2])).c_str());
    //cout<<"CHI QUADRO = "<<par[2]<<endl;
    TGraph *g = new TGraph();
    int k = 0;
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p1) && (centroid.at(i) != 0))
      {
        //cout<<centroid.at(i)<<endl;
        g->SetPoint(k,x_axis[0],centroid.at(i));
        k++;
      }
    }
    if(k == 0)
      g->RemovePoint(0);
    g->SetPoint(k+1,x_axis[1],temp_centroid_2);
    g->SetPoint(k+2,x_axis[2],temp_centroid_3);
    g->SetPoint(k+3,x_axis[4],temp_centroid_4);
    TF1 *f = new TF1("f", "[1]*x + [0]", x_axis[0], x_axis[4]);
    f->SetParameter(0, par[0]);
    f->SetParameter(1, par[1]);
    g->SetMarkerStyle(20);
    g->SetMarkerColor(kRed);
    g->SetTitle(("chi quadro = "+to_string(par[2])).c_str());
    g->Draw("AP");
    f->Draw("same");
    c->SaveAs(("/home/leonardo/Scrivania/Roba/TBReco/Event_display/event_display_" + to_string(par[2]) + ".pdf").c_str());
    //c->Print(("event_display_" + to_string(par[2]) + ".pdf").c_str(),"pdf");
  }
  */
  if((par[1]<=0.404 && par[1] >= -0.404))
  {
    den_fit_eff[0] = den_fit_eff[0] + 1;
    for(int i = 0; i<5; i++)
    {
      max_road[i] = par[0]+13+par[1]*(x_axis[i]);
      min_road[i] = par[0]-13+par[1]*(x_axis[i]);
    }
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p1) && (centroid.at(i) != 0))
        h2_extrapolate_pos_fit[0]->Fill(par[0]+par[1]*x_axis[0],centroid.at(i));
      if((chamb.at(i) == _ch_p1) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[0] && centroid.at(i)>=min_road[0]))
      {
        num_fit_eff[0] = num_fit_eff[0] + 1;
        h_residues_pos_fit[0]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[0]));
        break;
      }
    }  
  }
  else
    return;
}

void prova::Efficiency_Fit_2()
{
  double par[3];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, temp_centroid_3, temp_centroid_4;
  if(n_clu_1 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && DAQCount(_ch_p1) && DAQCount(_ch_p3) && DAQCount(_ch_p4))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        temp_centroid_3 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        temp_centroid_4 = centroid.at(k);
    }
  }
  else
    return;
  FitParameters(temp_centroid_1, temp_centroid_3, temp_centroid_4, x_axis[0], x_axis[2], x_axis[4], par);
  if((par[1]<=0.404 && par[1] >= -0.404))
  {
    den_fit_eff[1] = den_fit_eff[1] + 1;
    for(int i = 0; i<5; i++)
    {
      max_road[i] = par[0]+13+par[1]*(x_axis[i]);
      min_road[i] = par[0]-13+par[1]*(x_axis[i]);
    }
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p2) && (centroid.at(i) != 0))
        h2_extrapolate_pos_fit[1]->Fill(par[0]+par[1]*x_axis[1],centroid.at(i));
      if((chamb.at(i) == _ch_p2) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[1] && centroid.at(i)>=min_road[1]))
      {
        num_fit_eff[1] = num_fit_eff[1] + 1;
        h_residues_pos_fit[1]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[1]));
        break;
      }
    }
  }
  else
    return;
}

void prova::Efficiency_Fit_3()
{
  double par[3];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, temp_centroid_3, temp_centroid_4;
  if(n_clu_1 == 1 && n_clu_2 == 1 && n_clu_4 == 1 && DAQCount(_ch_p1) && DAQCount(_ch_p2) && DAQCount(_ch_p4))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        temp_centroid_4 = centroid.at(k);
    }
  }
  else
    return;
  FitParameters(temp_centroid_1, temp_centroid_2, temp_centroid_4, x_axis[0], x_axis[1], x_axis[4], par);
  if((par[1]<=0.404 && par[1] >= -0.404))
  {
    den_fit_eff[2] = den_fit_eff[2] + 1;
    for(int i = 0; i<5; i++)
    {
      max_road[i] = par[0]+13+par[1]*(x_axis[i]);
      min_road[i] = par[0]-13+par[1]*(x_axis[i]);
    }
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p3) && (centroid.at(i) != 0))
        h2_extrapolate_pos_fit[2]->Fill(par[0]+par[1]*x_axis[2],centroid.at(i));
      if((chamb.at(i) == _ch_p3) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[2] && centroid.at(i)>=min_road[2]))
      {
        num_fit_eff[2] = num_fit_eff[2] + 1;
        h_residues_pos_fit[2]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[2]));
        break;
      }
    }  
  }
  else
    return;
}

void prova::Efficiency_Fit_4()
{
  double par[3];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, temp_centroid_3, temp_centroid_4;
  if(n_clu_1 == 1 && n_clu_2 == 1 &&  n_clu_3 == 1 && DAQCount(_ch_p1) && DAQCount(_ch_p2) && DAQCount(_ch_p3))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        temp_centroid_3 = centroid.at(k);
    }
  }
  else
    return;
  FitParameters(temp_centroid_1, temp_centroid_2, temp_centroid_3, x_axis[0], x_axis[1], x_axis[2], par);
  if((par[1]<=0.404 && par[1] >= -0.404))
  {
    den_fit_eff[3] = den_fit_eff[3] + 1;
    for(int i = 0; i<5; i++)
    {
      max_road[i] = par[0]+13+par[1]*(x_axis[i]);
      min_road[i] = par[0]-13+par[1]*(x_axis[i]);
    }
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p4) && (centroid.at(i) != 0))
        h2_extrapolate_pos_fit[3]->Fill(par[0]+par[1]*x_axis[4],centroid.at(i));
      if((chamb.at(i) == _ch_p4) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[4] && centroid.at(i)>=min_road[4]))
      {
        num_fit_eff[3] = num_fit_eff[3] + 1;
        h_residues_pos_fit[3]->Fill(centroid.at(i)-(par[0]+par[1]*x_axis[4]));
        break;
      }
    }
  }
  else
    return;
}

void prova::FitParameters(double &y1, double &y2, double &y3, double &x1, double &x2, double &x3, double par[3])
{
  TGraph *g = new TGraph();
  g->SetPoint(0,x1, y1);
  g->SetPoint(1,x2, y2);
  g->SetPoint(2,x3, y3);
  TF1 *f = new TF1("f", "[1]*x+[0]", x1, x3);
  g->Fit(f, "QR");
  par[0] = f->GetParameter(0);
  par[1] = f->GetParameter(1);
  par[2] = f->GetChisquare();
}

void prova::Efficiency_Position_1()
{
  double par[2];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, super_eta = 0;
  if(n_clu_2 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && DAQCount(_ch_p2) && DAQCount(_ch_p3) && DAQCount(_ch_p4))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
    }
    super_eta /= (2*cos(angle));
  }
  else
    return;
  parameters(temp_centroid_2, super_eta, x_axis[1], x_axis[3], par);
  if(par[1]<=0.404 && par[1]>= -0.404)
  {
    den_eff_pos[0]->Fill(par[0]+par[1]*x_axis[0]);
    max_road[0] = par[0]+5+par[1]*x_axis[0];
    min_road[0] = par[0]-5+par[1]*x_axis[0];
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p1) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[0] && centroid.at(i)>=min_road[0]))
      {
        num_eff_pos[0]->Fill(par[0]+par[1]*x_axis[0]);
        break;
      }
    }
  }
  else
    return;
}

void prova::Efficiency_Position_2()
{
  double par[2];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[5], max_road[5];
  double temp_centroid_1, temp_centroid_2, super_eta = 0;
  if(n_clu_1 == 1 && n_clu_3 == 1 && n_clu_4 == 1 && DAQCount(_ch_p1) && DAQCount(_ch_p3) && DAQCount(_ch_p4))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 3) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 4) && (centroid.at(k)!=0))
        super_eta += centroid.at(k);
    }
    super_eta /= (2*cos(angle));
  }
  else
    return;
  parameters(temp_centroid_1, super_eta, x_axis[0], x_axis[3], par);
  if(par[1]<=0.404 && par[1]>= -0.404)
  {
    den_eff_pos[1]->Fill(par[0]+par[1]*x_axis[1]);
    max_road[1] = par[0]+5+par[1]*x_axis[1];
    min_road[1] = par[0]-5+par[1]*x_axis[1];
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p2) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[1] && centroid.at(i)>=min_road[1]))
      {
        num_eff_pos[1]->Fill(par[0]+par[1]*x_axis[1]);
        break;
      }
    }
  }
  else
    return;
}

void prova::Efficiency_Position_Stereo()
{
  double par[2];
  double x_axis[5] = {1.,12.562,33.08,38.86,44.64};
  double min_road[4], max_road[4];
  double temp_centroid_1, temp_centroid_2;
  if(n_clu_1 == 1 && n_clu_2 == 1 && DAQCount(_ch_p1) && DAQCount(_ch_p2))
  {
    eff_evt++;
    for(UInt_t k =0; k<centroid.size(); k++)
    {
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 1) && (centroid.at(k)!=0))
        temp_centroid_1 = centroid.at(k);
      if((((int)(atoi((chamb.at(k)).substr(1,1).c_str()))) == 2) && (centroid.at(k)!=0))
        temp_centroid_2 = centroid.at(k);
    }
  }
  else
    return;
  //double start_point = (temp_centroid_1+temp_centroid_2)/2;
  parameters(temp_centroid_1, temp_centroid_2, x_axis[0], x_axis[1], par);
  if(par[1]<=0.404 && par[1]>= -0.404)
  {
    den_eff_pos[2]->Fill(par[0]+par[1]*x_axis[2]);
    den_eff_pos[3]->Fill(par[0]+par[1]*x_axis[4]);
    //max_road[2] = par[0]+Stereo_Road(start_point)+par[1]*x_axis[2];
    //min_road[2] = par[0]-Stereo_Road(start_point)+par[1]*x_axis[2];
    //max_road[3] = par[0]+Stereo_Road(start_point)+par[1]*x_axis[4];
    //min_road[3] = par[0]-Stereo_Road(start_point)+par[1]*x_axis[4];
    max_road[2] = par[0]+5+par[1]*x_axis[2];
    min_road[2] = par[0]-5+par[1]*x_axis[2];
    max_road[3] = par[0]+5+par[1]*x_axis[4];
    min_road[3] = par[0]-5+par[1]*x_axis[4];
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p3) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[2] && centroid.at(i)>=min_road[2]))
      {
        num_eff_pos[2]->Fill(par[0]+par[1]*x_axis[2]);
        break;
      }
    }
    for(UInt_t i = 0; i<centroid.size(); i++)
    {
      if((chamb.at(i) == _ch_p4) && (centroid.at(i)!= 0) && (centroid.at(i)<=max_road[3] && centroid.at(i)>=min_road[3]))
      {
        num_eff_pos[3]->Fill(par[0]+par[1]*x_axis[4]);
        break;
      }
    }
  }
  else
    return;
}

bool prova::DAQCount(string &name)
{
  for(UInt_t i = 0; i<chamb.size(); i++)
  {
    if(chamb.at(i) == name && cluster_max_q.at(i)<100)
    {
      n_volte++;
      return false;
    }
  }
  return true;
}

int prova::Stereo_Road(double &start_point)
{
  double base = (major_base/2)-start_point*tan(phi);
  return base*tan(angle)+sigma_road;
  //return 5;
}

prova::prova(TTree *tree) : fChain(0)
{
  // if parameter tree is not specified (or zero), connect the file
  // used to generate this class and read the Tree.
  if(tree == nullptr)
  {
    TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("data/run177.root");
    if (!f || !f->IsOpen())
    {
      f = new TFile("data/run177.root");
    }
    f->GetObject("apv_raw", tree);
  }
  Init(tree);
}

prova::~prova()
{
  if (!fChain)
    return;
  delete fChain->GetCurrentFile();
}

Int_t prova::GetEntry(Long64_t entry)
{
  // Read contents of entry.
  if (!fChain)
    return 0;
  return fChain->GetEntry(entry);
}

Long64_t prova::LoadTree(Long64_t entry)
{
  // Set the environment to read one entry
  if (!fChain)
    return -5;
  Long64_t centry = fChain->LoadTree(entry);
  if (centry < 0)
    return centry;
  if (fChain->GetTreeNumber() != fCurrent)
  {
    fCurrent = fChain->GetTreeNumber();
    Notify();
  }
  return centry;
}

void prova::Init(TTree *tree)
{
  // The Init() function is called when the selector needs to initialize
  // a new tree or chain. Typically here the branch addresses and branch
  // pointers of the tree will be set.
  // It is normally not necessary to make changes to the generated
  // code, but the routine can be extended by the user if needed.
  // Init() will be called many times when running on PROOF
  // (once per file to be processed).

  // Set object pointer
  srsFec = 0;
  srsChip = 0;
  srsChan = 0;
  mmChamber = 0;
  mmLayer = 0;
  mmReadout = 0;
  mmStrip = 0;
  raw_q = 0;
  max_q = 0;
  t_max_q = 0;
   
  // Set branch addresses and branch pointers
  if (!tree)
    return;
  fChain = tree;
  fCurrent = -1;
  fChain->SetMakeClass(1);

  fChain->SetBranchAddress("evt", &evt, &b_evt);
  fChain->SetBranchAddress("error", &error, &b_error);
  fChain->SetBranchAddress("daqTimeSec", &daqTimeSec, &b_daqTimeSec);
  fChain->SetBranchAddress("daqTimeMicroSec", &daqTimeMicroSec, &b_daqTimeMicroSec);
  fChain->SetBranchAddress("srsTimeStamp", &srsTimeStamp, &b_srsTimeStamp);
  fChain->SetBranchAddress("srsTrigger", &srsTrigger, &b_srsTrigger);
  fChain->SetBranchAddress("srsFec", &srsFec, &b_srsFec);
  fChain->SetBranchAddress("srsChip", &srsChip, &b_srsChip);
  fChain->SetBranchAddress("srsChan", &srsChan, &b_srsChan);
  fChain->SetBranchAddress("mmChamber", &mmChamber, &b_mmChamber);
  fChain->SetBranchAddress("mmLayer", &mmLayer, &b_mmLayer);
  fChain->SetBranchAddress("mmReadout", &mmReadout, &b_mmReadout);
  fChain->SetBranchAddress("mmStrip", &mmStrip, &b_mmStrip);
  fChain->SetBranchAddress("raw_q", &raw_q, &b_raw_q);
  fChain->SetBranchAddress("max_q", &max_q, &b_max_q);
  fChain->SetBranchAddress("t_max_q", &t_max_q, &b_t_max_q);
  /*
  h_sum_q_cluster_ch1 = new TH1F("h_sum_q_cluster_ch1","Distribution sum of charge in cluster P1_M01;sum_q_P1_M01;Entry",150,0.0,5000.);
  h_sum_q_cluster_ch2 = new TH1F("h_sum_q_cluster_ch2","Distribution sum of charge in cluster P2_M01;sum_q_P2_M01;Entry",150,0.0,5000.);
  h_sum_q_cluster_ch3 = new TH1F("h_sum_q_cluster_ch3","Distribution sum of charge in cluster P3_M01;sum_q_P3_M01;Entry",150,0.0,5000.);
  h_sum_q_cluster_ch4 = new TH1F("h_sum_q_cluster_ch4","Distribution sum of charge in cluster P4_M01;sum_q_P4_M01;Entry",150,0.0,5000.);
  h_sum_q_centroid_ch1 = new TH2F("h_sum_q_centroid_ch1","Histo sum of charge vs centroid P1_M01;sum_q;Centroid",150,0.0,5000.,512,500.,1124.);
  h_sum_q_centroid_ch2 = new TH2F("h_sum_q_centroid_ch2","Histo sum of charge vs centroid P2_M01;sum_q;Centroid",150,0.0,5000.,512,500.,1124.);
  h_sum_q_centroid_ch3 = new TH2F("h_sum_q_centroid_ch3","Histo sum of charge vs centroid P3_M01;sum_q;Centroid",150,0.0,5000.,512,500.,1124.);
  h_sum_q_centroid_ch4 = new TH2F("h_sum_q_centroid_ch4","Histo sum of charge vs centroid P4_M01;sum_q;Centroid",150,0.0,5000.,512,500.,1124.);
  */
  Notify();
}

Bool_t prova::Notify()
{
  // The Notify() function is called when a new file is opened. This
  // can be either for a new TTree in a TChain or when when a new TTree
  // is started when using PROOF. It is normally not necessary to make changes
  // to the generated code, but the routine can be extended by the
  // user if needed. The return value is currently not used.

  return kTRUE;
}

void prova::Show(Long64_t entry)
{
  // Print contents of entry.
  // If entry is not specified, print current entry
  if (!fChain)
    return;
  fChain->Show(entry);
}

Int_t prova::Cut(Long64_t entry)
{
  // This function may be called from Loop.
  // returns  1 if entry is accepted.
  // returns -1 otherwise.
  return 1;
}
#endif