/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2010/Testbeam_june_10/Gassiplex/Telescope/Track_finder_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 = "b0624_5_acq_HR2_24062010_0634_merged.root";
  TString ped_file = "p0624_1.dat";

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

  int event = 0;
  int nsigma = 2;




  /*Variables*/

  int nsuper = 0;
  int ntrack = 0;
  int nstraight = 0;
  int thr = 0;
  int adc = 0;
  int xpos = 0;
  int ypos = 0;
  int nchannel = 0;
  int xadc_max_all = 0;
  int yadc_max_all = 0;
  int content = 0;
  int nquiet = 0;
  int nhot = 0;

  float nEvent = 0;

  bool hit = 0;
  bool inside = 0;
  bool quiet = 0;
  bool aligned = 0;

  bool hot[4];
  int nhit[4];
  int xadc_max[4];
  int yadc_max[4];
  float ped_rms[4][96*4];
  float ped_mean[4][96*4];

  TH1F * hhot = new TH1F("hhot","Hot chambers",4,0,4);
  TH1F * hnhot = new TH1F("hnhot","Number of hot chambers per event",5,0,5);
  TH1F * hnhot_cumul = new TH1F("hnhot_cumul","Cumulative of number of hot chambers per event",5,0,5);

  TH1F * hx_raw = new TH1F("hx_raw","Raw hit pattern in X direction;x (cm)",32,0,32);
  TH1F * hy_raw = new TH1F("hy_raw","Raw hit pattern in Y direction;y (cm)",32,0,32);

  TH1F * hx_clean = new TH1F("hx_clean","Hit pattern in X direction after noise removal;x (cm)",32,0,32);
  TH1F * hy_clean = new TH1F("hy_clean","Hit pattern in Y direction noise removal;y (cm)",32,0,32);

  TH1F * hx_adc_all = new TH1F("hx_adc_all","ADC pattern in X direction;x (cm); ADC counts",32,0,32);
  TH1F * hy_adc_all = new TH1F("hy_adc_all","ADC pattern in Y direction;y (cm); ADC counts",32,0,32);

  TH1F * hthr[4];
  TH1F hadc_all[4];
  TH1F * hx_adc[4];
  TH1F * hy_adc[4];
  TH1F * hnhit_raw[4];
  TH1F * hnhit_clean[4];
  TH1F * hnhit_gold[4];
  TH1F * hnhit_raw_cumul[4];
  TH2F * hxy_raw[4];
  TH2F * hxy_clean[4];
  TString title,name;

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

    hot[i] = 0;
    nhit[i] = 0;
    xadc_max[i] = 0;
    yadc_max[i] = 0;

    title = Form("hthr_%i",i+1);
    name = Form("Threshold at 2#sigma in chamber %i;thr (fC)",i+1);
    hthr[i] = new TH1F(title,name,500,0,10);

    title = Form("hadc_all_%i",i+1);
    name = Form("ADC distribution in chamber %i;ADC counts",i+1);
    hadc_all[i] = new TH1F(title,name,1024,0,1024);

    title = Form("hx_adc_%i",i+1);
    name = Form("ADC pattern in chamber %i after noise removal;x (cm); ADC counts",i+1);
    hx_adc[i] = new TH1F(title,name,32,0,32);

    title = Form("hy_adc_%i",i+1);
    name = Form("ADC pattern in chamber %i after noise removal;y (cm);ADC counts",i+1);
    hy_adc[i] = new TH1F(title,name,32,0,32);
    
    title = Form("hnhit_clean_%i",i+1);
    name = Form("Number of hit in chamber %i after noise removal;nhit",i+1);
    hnhit_clean[i] = new TH1F(title,name,96*4,0,96*4);

    title = Form("hnhit_gold_%i",i+1);
    name = Form("Number of hit in chamber %i for golden events;nhit",i+1);
    hnhit_gold[i] = new TH1F(title,name,96*4,0,96*4);

    title = Form("hnhit_raw_%i",i+1);
    name = Form("Number of hit in chamber %i;nhit",i+1);
    hnhit_raw[i] = new TH1F(title,name,96*4,0,96*4);

    title = Form("hnhit_raw_cumul_%i",i+1);
    name = Form("Cumulative of number of hit in chamber %i;nhit",i+1);
    hnhit_raw_cumul[i] = new TH1F(title,name,96*4,0,96*4);

    title = Form("hxy_raw_%i",i+1);
    name = Form("Raw hit pattern in chamber %i;y (cm);x (cm)",i+1);
    hxy_raw[i] = new TH2F(title,name,32,0,32,32,0,32);

    title = Form("hxy_clean_%i",i+1);
    name = Form("Hit pattern in chamber %i after noise removal;y (cm);x (cm)",i+1);
    hxy_clean[i] = new TH2F(title,name,32,0,32,32,0,32);

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

      ped_mean[i][k] = 0;
      ped_rms[i][k] = 0;}}




  /*Read pedestal file*/

  int chb,chn,thr_adc,x,y;
  float pede,pede_err,noise,noise_err,gain;

  cout<<"Read file "<<ped_file<<endl;
  ifstream in(ped_file.Data());
  while(in>>chn>>x>>y>>chb>>pede>>pede_err>>noise>>noise_err>>gain){
    hthr[chb-1]->Fill(noise*2/gain);
    ped_mean[chb-1][chn] = pede;
    ped_rms[chb-1][chn] = noise;}






  /*Read data file and look for hits*/

  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++){
      //for (int i=event;i<event+1;i++){

      /*Reset histos and variables*/

      nhot = 0;
      hx_adc_all->Reset();
      hy_adc_all->Reset();
      for (int k=0;k<4;k++){
        hot[k] = 0;
        nhit[k] = 0;
        xadc_max[k] = 0;
        yadc_max[k] = 0;
        hx_adc[k]->Reset();
        hy_adc[k]->Reset();
        hxy_clean[k]->Reset();
        hxy_raw[k]->Reset();}

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

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

        /*Get channel info*/

        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
        chb = channel->GetChamberId()-1;

        if (chb!=4){

        xpos = channel->GetX();
        ypos = channel->GetY();

        /*Test for hit*/

        inside = (ypos>=6) && (ypos<=22) && (xpos>=3) && (xpos<=9);

        if (inside){

          chn = channel->GetHardId();
          adc = channel->GetAnalogValue();

          hadc_all[chb]->Fill(adc);
          thr = floor(ped_mean[chb][chn] + nsigma * ped_rms[chb][chn] + 1);
          hit = ((adc + ped_mean[chb][chn] - 20) >= thr);

          if (hit){

            hx_adc_all->Fill(xpos,adc);
            hy_adc_all->Fill(ypos,adc);
            hxy_raw[chb]->Fill(ypos,xpos,adc);}}}}

      for (int k=0;k<4;k++){hnhit_raw[k]->Fill(hxy_raw[k]->GetEntries());}

      /*Look events with less than 10 hits in each chambers*/

      quiet = ((hxy_raw[0]->GetEntries()<10) && (hxy_raw[1]->GetEntries()<10) && (hxy_raw[2]->GetEntries()<10) && (hxy_raw[3]->GetEntries()<10));

      if (quiet){

        nquiet++;

        /*Identify 3x3 region of track candidate*/

        xadc_max_all = hx_adc_all->GetMaximumBin()-1;
        yadc_max_all = hy_adc_all->GetMaximumBin()-1;

        /*Keep hits inside region*/

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

        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
        chb = channel->GetChamberId()-1;

        if (chb!=4){

          /*Get channel info*/

          xpos = channel->GetX();
          ypos = channel->GetY();

          /*Test for hit*/

          inside = ((fabs(xpos-xadc_max_all)<=1) && (fabs(ypos-yadc_max_all)<=1));

          if (inside){

            chb = channel->GetChamberId()-1;
            chn = channel->GetHardId();
            adc = channel->GetAnalogValue();

            thr = floor(ped_mean[chb][chn] + nsigma * ped_rms[chb][chn] + 1);
            hit = ((adc + ped_mean[chb][chn] - 20) >= thr);

            if (hit){

              hx_adc[chb]->Fill(xpos,adc);
              hy_adc[chb]->Fill(ypos,adc);

              hxy_clean[chb]->Fill(ypos,xpos,adc);
              nhit[chb]++;}}}}

        /*Identify hot chambers and max of ADC pattern*/

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

          if (nhit[k]!=0){

            hot[k] = 1;
            hhot->Fill(k);
            nhot++;

            xadc_max[k] = hx_adc[k]->GetMaximumBin()-1;
            yadc_max[k] = hy_adc[k]->GetMaximumBin()-1;}

          hnhit_clean[k]->Fill(nhit[k]);}

        hnhot->Fill(nhot);


        /*Measure gross telescope performance*/

        if (hot[0] && hot[1] && hot[2]){hnhit_gold[3]->Fill(nhit[3]);}
        if (hot[0] && hot[1] && hot[3]){hnhit_gold[2]->Fill(nhit[2]);}
        if (hot[0] && hot[2] && hot[3]){hnhit_gold[1]->Fill(nhit[1]);}
        if (hot[1] && hot[2] && hot[3]){hnhit_gold[0]->Fill(nhit[0]);}


        /*Select straight tracks*/

        if (hot[0] && hot[3] && (hot[1] || hot[2])){

          ntrack++;

          if ((xadc_max[0] == xadc_max[3]) && (yadc_max[0] == yadc_max[3])){

            nstraight++;

            if (nhit[0]==1 && nhit[3]==1){

              nsuper++;}}}}}}



  

  /*Calculate cumulatives*/

    hnhot_cumul->SetBinContent(5,hnhot->GetBinContent(5));

    for (int j=3;j>=0;j--){
      content = hnhot_cumul->GetBinContent(j+2) + hnhot->GetBinContent(j+1);
      hnhot_cumul->SetBinContent(j+1,content);}

    hnhot_cumul->Scale(1/(1.*hnhot->GetEntries()));

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

      hnhit_raw_cumul[k]->SetBinContent(1,hnhit_raw[k]->GetBinContent(1));

      for (int j=1;j<96;j++){

        content = hnhit_raw_cumul[k]->GetBinContent(j) + hnhit_raw[k]->GetBinContent(j+1);
        hnhit_raw_cumul[k]->SetBinContent(j+1,content);}
      
      hnhit_raw_cumul[k]->Scale(1/(1.*hnhit_raw[k]->GetEntries()));}






    cout<<endl;
    cout<<"Fraction of events with less than 10 hits in each chamber = "<<nquiet/nEvent*100.<<" %"<<endl;
    cout<<"  with 3 or 4 hot chambers = "<<hnhot->Integral(4,5)/nEvent*100.<<" %"<<endl;
    cout<<"    with chamber 1 and 4 hot = "<<ntrack/nEvent*100.<<" %"<<endl;
    cout<<"      with max in chamber 1 and 4 aligned = "<<nstraight/nEvent*100.<<" %"<<endl;
    cout<<"        with single hits in chamber 1 and 4 = "<<nsuper/nEvent*100.<<" %"<<endl;
    cout<<endl;


    float eff[4];
    float eff_err[4];
    float mult[4];

    float fnt = 0;
    float fn0 = 0;
    float fm = 0;
    float fnh = 0;

    cout<<"Telescope performance with golden events"<<endl;
    cout<<endl;

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

      fm = 0;
      fnt = hnhit_gold[i]->GetEntries();
      fn0 = hnhit_gold[i]->GetBinContent(1);
      fnh = hnhit_gold[i]->Integral(2,10);

      for (k=2;k<10;k++){fm += (k-1)*hnhit_gold[i]->GetBinContent(k)/fnh;}

      cout<<"Chamber "<<i+1<<endl;
      cout<<"eff  = "<<(1-fn0/fnt)*100.<<" +/- "<<100.*sqrt(fn0/fnt/fnt/fnt*(fnt-fn0))<<" %"<<endl;
      cout<<"mult = "<<fm<<endl;
    }

    cout<<endl;








}

---