/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2010/Testbeam_june_10/Gassiplex/Telescope/Pedestals_study.C

Go to the documentation of this file.

{
  gROOT->Reset();
  gStyle->SetPalette(1);
  gSystem->Load("/home/chefdevi/micromegas_framework_03032011/trunk/libMicro.so");


  TString file = "p0624_1.root";
  TString filename = "/home/chefdevi/lapp_data/LC/Detecteurs/MicroMegas/data/TB2010/SPS_H4_june_2010/Root_files/pro/"+file;






  /*Variables*/

  float nEvent = 0;
  int nchannel = 0;
  int x_channel[4][96*4];
  int y_channel[4][96*4];
  float ped_mean[4][96*4];
  float ped_rms[4][96*4];
  float ped_mean_err[4][96*4];
  float ped_rms_err[4][96*4];

  TH1F * hped[4][96*4];
  TH1F * hc_ped[4];
  TH2F * hxy_noise[4];
  TH1F * hnoise[4];
  TH1F * hc_noise[4];
  TString title,name;

  for (int i=0;i<4;i++){

    title = Form("hnoise_%i",i+1);
    name = Form("Noise in chamber %i;noise (fC)",i+1);
    hnoise[i] = new TH1F(title,name,1000,0,10);

    title = Form("hc_noise_%i",i+1);
    name = Form("Noise per channel in chamber %i;channel;noise (fC)",i+1);
    hc_noise[i] = new TH1F(title,name,96*4,0,96*4);

    title = Form("hxy_noise_%i",i+1);
    name = Form("Noise per channel in chamber %i;y (cm);x (cm);noise (fC)",i+1);
    hxy_noise[i] = new TH2F(title,name,32,0,32,32,0,32);

    title = Form("hc_ped_%i",i+1);
    name = Form("Pedestal per channel in chamber %i;pedestal (ADC)",i+1);
    hc_ped[i] = new TH1F(title,name,96*4,0,96*4);

    for (int j=0;j<96*4;j++){

      ped_mean[i][j] = 0;
      ped_rms[i][j] = 0;
      ped_mean_err[i][j] = 0;
      ped_rms_err[i][j] = 0;
      x_channel[i][j] = 0;
      y_channel[i][j] = 0;
      title = Form("hped_%i_%i",i+1,j); 
      name = Form("Pedestal position distribution in chamber %i channel %i;pedetsal (ADC)",i+1,j); 
      hped[i][j] = new TH1F(title,name,1024,0,1024);}}




  /*Read pedestal file*/

  TFile *f = new TFile(filename);
  cout<<"reading File : "<<filename<<endl;
  TIter nextkey(f.GetListOfKeys());
  TKey *key;

  while (key = (TKey*)nextkey()) {
    
    TTree *tree = (TTree*)key->ReadObj();                
    MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

    nEvent = tree.GetEntries();
    cout <<"nEvent="<<nEvent<<endl;

    MTEvent *evt = new MTEvent();
    TBranch *branch= tree->GetBranch("MTEvent");
    branch->SetAddress(&evt);

    for (int i=0;i<nEvent;i++){  

      tree.GetEntry(i);
      nchannel = evt->GetNchannel();
      if (i%100==0){cout<<i/(1.*nEvent)*100.<<" % \r"<<flush;}

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

        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
        hped[channel->GetChamberId()-1][channel->GetHardId()]->Fill(channel->GetAnalogValue());
        hc_ped[channel->GetChamberId()-1]->SetBinContent(channel->GetHardId()+1,channel->GetAnalogValue());

        if (i==1){
          x_channel[channel->GetChamberId()-1][channel->GetHardId()] = channel->GetX();
          y_channel[channel->GetChamberId()-1][channel->GetHardId()] = channel->GetY();}}}}





  /*Read preamp gain file*/

  ifstream in("gassiplex_preamp_gains.dat");
  float gain[4][96*4];
  float g = 0;

  for (int i=0;i<4;i++){
    for (int j=0;j<96*4;j++){
      gain[i][j] = 0;}}

  while(in>>j>>g){
    if (j<96){gain[0][j] = g;}
    else{if (j>=96 && j<96*2){gain[1][j-96] = g;}
      else{if (j>=96*2 && j<96*3){gain[2][j-96*2] = g;}
        else{gain[3][j-96*3] = g;}}}}





  /*Calculate noise*/

  TF1 * fped = new TF1("fped","gaus",0,1024);

  for (int i=0;i<4;i++){
    for (int j=0;j<96*4;j++){

      if (hped[i][j]->GetEntries()!=0){

        fped->SetParameters(hped[i][j]->GetEntries(),hped[i][j]->GetMean(),hped[i][j]->GetRMS());
        hped[i][j]->Fit(fped,"q");
        ped_mean[i][j] = fped->GetParameter(1);
        ped_mean_err[i][j] = fped->GetParError(1);
        ped_rms[i][j] = fped->GetParameter(2);
        ped_rms_err[i][j] = fped->GetParError(2);

        if (gain[i][j]!=0){
          hnoise[i]->Fill(ped_rms[i][j]/gain[i][j]);
          hc_noise[i]->SetBinContent(j+1,ped_rms[i][j]/gain[i][j]);
          hxy_noise[i]->SetBinContent(y_channel[i][j]+1,x_channel[i][j]+1,ped_rms[i][j]/gain[i][j]);}

        cout<<i<<" "<<x_channel[i][j]<<" "<<y_channel[i][j]<<" "<<ped_mean[i][j]<<" "<<ped_rms[i][j]<<" \r "<<flush;}}}

  for (int i=0;i<4;i++){

    fped->SetParameters(hnoise[i]->GetEntries(),hnoise[i]->GetMean(),hnoise[i]->GetRMS());
    hnoise[i]->Fit(fped,"q");}




  /*Calculate and write threshold in ADC units*/

  TString thr_file;
  float thr_adc = 0;
  float thr_fc = 0;

  int Nthr0 = 0;
  int Nthr = 50;
  TH1F * hthr[4][Nthr];
  TH1F * hc_thr[4][Nthr];

  for (i=0;i<4;i++){file.Chop();}

  TString ped_file = file+"dat";
  ofstream out_ped(ped_file.Data());

  for (int i=0;i<4;i++){
    for (int j=0;j<96*4;j++){
        if (gain[i][j]!=0){
          out_ped<<j<<" "<<x_channel[i][j]<<" "<<y_channel[i][j]<<" "<<i+1<<" "<<ped_mean[i][j]<<" "<<ped_mean_err[i][j]<<" "<<ped_rms[i][j]<<" "<<ped_rms_err[i][j]<<" "<<gain[i][j]<<endl;}}}



  /*Plots*/

  if (0){

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

    TCanvas * c0 = new TCanvas("c0","",10,10,500,500);

    for (i=0;i<4;i++){
      title = Form("chb. %i",i+1);
      hc_noise[i]->SetTitle(title);
      hc_noise[i]->SetLineWidth(3);
      hc_noise[i]->GetYaxis()->SetRangeUser(0,1.2);
      if (i<3){hc_noise[i]->GetXaxis()->SetRangeUser(0,96);}


      hc_noise[i]->Draw();
      title_eps = Form("Plots/noise_chb_%i.eps",i+1);
      title_png = Form("Plots/noise_chb_%i.png",i+1);
      c0->Print(title_eps);
      c0->Print(title_png);}



  }




}

---