/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/Microroc/Calibration/Fit_scurve_chip.C

Go to the documentation of this file.

{gROOT->Reset();
  gROOT->SetStyle("Video");
  gStyle->SetOptStat(0);



  bool display = 1;
  bool write = 1;
  bool write_in_tfile = 0;
  bool plot = 1;

  TString dac = "2";

  TString dir = "/lapp_data/LC/Detecteurs/MicroMegas/data/MICROROC/MR1/Root_files/Calibration_Dans_M2/";
  TString file = "scurve_noise_offset7_dac0_asu1314";
  TString t_file = dir + file + ".root";

  cout<<t_file<<endl;




  if (write_in_tfile){TFile * tf = new TFile(t_file,"UPDATE");}
  else{TFile * tf = new TFile(t_file);}


  
  int threshold = 0;
  int chipid = 0;
  int hardid = 0;
  int value = 0;
  int content = 0;
  
  int tmin = 0;
  int tmed = 0;
  int tmax = 0;
  int delta = 0;

  TH1I * hscurve[64];

  TH1F * hinflex = new TH1F("hinflex","",1000,0,500);
  TH1F * hsigma = new TH1F("hsigma","",50,-2,2);

  TH1F * hinflex_err = new TH1F("hinflex","",100,0,10);
  TH1F * hsigma_err = new TH1F("hsigma","",100,0,2);

  TH1F * hc_inflex = new TH1F("hc_inflex","",64,0,64);
  TH1F * hc_sigma = new TH1F("hc_sigma","",64,0,64);

  TH1F * hc_inflex_err = new TH1F("hc_inflex_err","",64,0,64);
  TH1F * hc_sigma_err = new TH1F("hc_sigma_err","",64,0,64);

  TH1F * hc_inflex_sigma = new TH1F("hc_inflex_sigma","",64,0,64);

  TString name;

  for (int j=0;j<64;j++){

    name  = Form("hscurve_%i",j);
    TH1I * h = (TH1I*)tf.Get(name);
    hscurve[j] = h;}



 


  TString func  = "[0]/(1+exp((x-[1])/[2]))";
  TF1 * fermi[64];
  float scurve_max = 0;
  float scurve_inflex = 0;
  float inflex = 0;
  float sigma = 0;
  float inflex_err = 0;
  float sigma_err = 0;
  float content0,content1,content2;
  float chi2 = 0;
  float ndf = 0;


  
  for (int j=0;j<64;j++){

    name  = Form("fermi_%i",j);
       
    fermi[j] = new TF1(name,func,0,1024);
    fermi[j]->SetLineColor(j+1);
    fermi[j]->SetLineWidth(0.5);
      
    /*If SCurve exist*/
    if (hscurve[j]->GetEntries()!=0){

      /*Find threshold min tmin*/
      for (int b=0;b<1024;b++){
        content = hscurve[j]->GetBinContent(b+1);
        if (content != 0){tmin = b;b = 1024;scurve_max = content;}}

      /*Find threshold at which scurve = 0 tmax*/
      for (int b=1024;b>tmin;b--){
        if (hscurve[j]->GetBinContent(b+1) != 0){tmax = b+1;b = tmin;}}
        
      /*Find threshold at which scurve reaches the max*/
      for (int b=tmax;b>=0;b--){
        if (hscurve[j]->GetBinContent(b+1) == scurve_max){tmed = b+1;b = -1;}}

      scurve_inflex = 0.5*(tmax+tmed);
      delta = tmax - tmed;
      
      /*If SCurve is a step function, do not fit parameters*/
      if (delta==1){

        inflex = scurve_inflex;
        sigma = 0.;
        inflex_err = 0.;
        sigma_err = 0.;}

      else{

        /*Set parameters*/
        fermi[j]->FixParameter(0,scurve_max);
        fermi[j]->SetParameter(1,scurve_inflex);
        fermi[j]->SetParameter(2,delta/10.);
        
        /*Perform fit*/
        hscurve[j]->Fit(fermi[j],"q","",0,tmax+50);

        /*Get parameters*/
        inflex = fermi[j]->GetParameter(1);
        sigma = sqrt(pow(fermi[j]->GetParameter(2),2));

        inflex_err = sqrt(pow(fermi[j]->GetParError(1),2));
        sigma_err = sqrt(pow(fermi[j]->GetParError(2),2));

        chi2 = fermi[j]->GetChisquare();
        ndf = fermi[j]->GetNDF();}


      hc_inflex_sigma->SetBinContent(j+1,inflex);
      hc_inflex_sigma->SetBinError(j+1,sigma);

      hc_inflex->SetBinContent(j+1,inflex);
      hc_sigma->SetBinContent(j+1,sigma);
      
      hc_inflex_err->SetBinContent(j+1,inflex_err);
      hc_sigma_err->SetBinContent(j+1,sigma_err);
      
      hinflex->Fill(inflex);
      hsigma->Fill(sigma);

      if (display){cout<<"Channel\t"<<j<<"\t"<<inflex<<"\t width\t"<<sigma<<"\r"<<flush;}}}




  
  if (write){

    TString outfile = dir + "parameters_" + file + ".dat";
    ofstream out(outfile.Data());

    for (int j=0;j<64;j++){

      out<<j<<" "<<hc_inflex->GetBinContent(j+1)<<" "<<hc_inflex_err->GetBinContent(j+1)<<" "<<hc_sigma->GetBinContent(j+1)<<" "<<hc_sigma_err->GetBinContent(j+1)<<endl;}}



  if (plot){

    TString adj;
    if (dac=="0"){adj = "bas";}
    if (dac=="1"){adj = "intermediaire";}
    if (dac=="2"){adj = "haut";}



    if (0){

    TCanvas * c0 = new TCanvas("c0","",10,10,500,500);
    TString title = "MICROROC pedestals - Seuil " + adj + ";threshold (DAC)";

    TH2F * frame = new TH2F("frame",title,100,0,1024,100,0,600);

    frame->GetXaxis()->SetRangeUser(140,200);
    frame->Draw();

    for (int j=0;j<64;j++){fermi[j]->Draw("same");}}}

  

  }

---