/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/geomVerif.cpp

Go to the documentation of this file.

#include <map>
#include <stdlib.h>
#include <string>
#include <set>
#include <list>
#include "TKey.h"
#include <TRint.h>
#include <TKey.h>
#include <TTree.h>
#include <TH2I.h>
#include <TProfile.h>
#include <TFile.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TChain.h>
#include <iostream>
#include <sstream>
#include "tools/Log.hh"
#include "tools/Toolbox.hh"


#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"

#include <string>

#include "MicroException.hh"

#include "mysql/Mysql.hh"

using namespace std;

#define STRIPWIDTH 250
#define PADWIDTH 10000
#define MAXSTRIPENTRY 10000
#define X 0
#define Y 1
    
//---------------------------------------------
//
//---------------------------------------------
int main(int argc, char**argv)

{

  /*m2 prototype time cut*/
  int dt_min = 6;
  int dt_max = 9;
  
  /*telescope energy cut*/
  int adc_cut = 30;

  if ( argc != 2  && argc != 3) {
    FILE_LOG(logERROR)  << "usage:" << endl;
    FILE_LOG(logERROR)  << "geomVerif rootFile [stripZRotation(90)]" << endl;
    exit(1);
  }


  string file = argv[1];

  bool rotation =  false;
  if ( argc == 3)
  {
    if ( atoi(argv[2]) == 90 )
    {rotation = true; }
    else { cout << " stripZRotation should be 90 or NULL "  << endl; exit(-1);}
  }
    

  TApplication *theApp;

  theApp = new TRint("App", &argc, argv);

  bool noisy[9][96];
  for (int i=0;i<9;i++){ for (int j=0;j<96;j++){ noisy[i][j] = 0;}}

  /*HV pad*/
  noisy[6][77] = 1;
  noisy[7][77] = 1;
  noisy[8][77] = 1;

  /*floating pad*/
  noisy[6][9] = 1;
  noisy[8][33] = 1;


  bool super_hot[10];
  int nsuper_hot_tel = 0;
  int nhit[10] = {0};
  bool align = 0;
  bool clusterInStrip[2] = {false,false};

  int nchannel = 0;
  int hardid = 0;
  int chipid = 0;

  Float_t xpos = 0.;  
  Float_t ypos = 0.;  
  Float_t zpos = 0.;  

  int adc = 0;
  int chbid = 0;
  int ntree = 0;

  Float_t pad_x = 0.;
  Float_t pad_y = 0.;

  Float_t xStripNum = 0.;  
  Float_t yStripNum = 0.;  

  Int_t pad[2][MAXSTRIPENTRY];
  Int_t strip[2][MAXSTRIPENTRY];
  Int_t nEntries[2] = {0,0};

  Long64_t t,t0,t1,t2,t3,dt;

  TH2I * hxy_tel_pad_temp = new TH2I("hxy_tel_pad_temp","Hit position distribution in telescope pad chamber for super hot events (Nhit=1);y (pad);x (pad)",32,0,32,32,0,32);
  TH2I * hx_cor = new TH2I("hx_cor","Correlation Prototype - telescope;Xtarget (pad);Xm2 (pad)",100,0,1000000,100,0,1000000);
  TH2I * hy_cor = new TH2I("hy_cor","Correlation Prototype - telescope;Ytarget (pad);Ym2 (pad)",100,0,1000000,100,0,1000000);
  TH1I * hx_res = new TH1I("hx_res","Vertical residuals Prototype - telescope;Xm2 - Xtarget (pad)",101,-500000,500000);
  TH1I * hy_res = new TH1I("hy_res","Horizontal residuals Prototype - telescope;Ym2 - Ytarget (pad)",101,-500000,500000);
  
 
  TFile *f = new TFile(file.c_str());
  cout<<endl;
  cout<<"Reading file : "<<file<<endl;
  cout<<endl;

  if (f == NULL) { exit(-1); }
  TIter nextkey(f->GetListOfKeys());
  TKey *key;

  while (key = (TKey*)nextkey())
  {

    TTree *tree = (TTree*)key->ReadObj();                
    if ( tree == NULL ) { exit(-1); }
    MTRun* run = (MTRun*)tree->GetUserInfo()->FindObject("MTRun");

    int nEvent = tree->GetEntries();
    ntree++;
    cout<<"  Tree "<<ntree<<" contains "<<nEvent<<" events"<<endl;

    MTEvent *evt = new MTEvent();

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

    for (int i=0;i<nEvent-1 ;i++)
    {
      if (i%100==0){cout<<"Progress: "<<i/1./nEvent*100<<" \r"<<flush;}
      set<Float_t> stripX;
      set<Float_t> stripY;
      clusterInStrip[X] = false;
      clusterInStrip[Y] = false;

      tree->GetEntry(i);
      nchannel = evt->GetNchannel();

      nsuper_hot_tel = 0;

      for (int k=0;k<10;k++){nhit[k] = 0; super_hot[k] = 0;}

      // Pour chaque evenement:
      // 1 On cherche dans le telescope 1 hit par chambre
      // Avec iun filtre sutr les chambres brillante
      // et une coupuresur la valeur de l'adc

      /*Look in telescope*/
      for (int j=0;j<nchannel;j++)
      {
        MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
        chbid = channel->GetChamberId();
        if (chbid < 10)
        {
          if (!noisy[chbid-1][hardid])
          {
            xpos = channel->GetX();//*1e-4;
            ypos = channel->GetY();//*1e-4;
            adc = channel->GetAnalogValue();
            if (adc>=adc_cut)
            {              
              nhit[chbid-1]++;
              if (chbid == 7 || chbid == 8 || chbid == 9)
              {
                hxy_tel_pad_temp->Fill(ypos,xpos);
              }
              if (chbid == 8)
              {
                pad_x = xpos;
                pad_y = ypos;
              }
              if (  (  chbid == 1 && rotation == false) || ( chbid == 2 && rotation == true)) // Strip gassiplex HORIZONTAL

              {
                 stripX.insert(xpos); 
              }
              if (  (  chbid == 2 && rotation == false) || ( chbid == 1 && rotation == true)) // Strip gassiplex HORIZONTAL
              {
                 stripY.insert(ypos); 
              }
            } // end if (adc>=adc_cut)
          }  // end if (!noisy[chbid-1][hardid])
        } // end  if (chbid < 10)
      } // end for (int j=0;j<nchannel;j++)

      // Find cluster in strip chamber
      // eliminate noise by filter on particle cluster formed with 3 to 5 consecutive strips
      if ( stripX.size() >= 3 &&  stripX.size() <=5 )
      {
  
// Over X
        clusterInStrip[X] = true;
        Float_t lastStrip = -1.;
        xStripNum = 0.;

        for( set<Float_t>::iterator ii=stripX.begin(); ii!=stripX.end(); ++ii)
        {
          Float_t foo = (*ii ) / STRIPWIDTH;
          xStripNum = xStripNum + foo;
          if ( lastStrip != -1 && foo != lastStrip + 1. )
          {
            clusterInStrip[X] = false; 
          }
          lastStrip = foo;
          
        }
        xStripNum = xStripNum / stripX.size();
      }

      if ( stripY.size() >= 3 &&  stripY.size() <=5 )
      {

// Over Y
        clusterInStrip[Y] = true;
        Float_t lastStrip = -1.;
        yStripNum = 0.;

        for( set<Float_t>::iterator ii=stripY.begin(); ii!=stripY.end(); ++ii)
        {
          Float_t foo = (*ii ) / STRIPWIDTH;
          yStripNum = yStripNum + foo;
          if ( lastStrip != -1 && foo != lastStrip + 1. )
          {
            clusterInStrip[Y] = false;
          }
          lastStrip = foo;

        }
        yStripNum = yStripNum / stripY.size();
 

      }
 

      //Find tracks in pad telescope
      //Only tracks which have only one Pad hit per chamber  
      for (int k=6;k<9;k++)
      {
        if (nhit[k] != 0)
        {
          if (nhit[k] == 1)
          {
           super_hot[k] = 1;
          }
        }

        if (super_hot[k])
        {
          nsuper_hot_tel++;
        }
      }


      //If track, look in m2 prototype
      if (nsuper_hot_tel == 3)
      {
        //GetRMS(1) == 0 sigifie que tous les points ont la mewme valeur sur l'axe 1
        align = (hxy_tel_pad_temp->GetRMS(1) == 0 && hxy_tel_pad_temp->GetRMS(2) == 0); // GetRMS(1) RMS sur l'axe 1 ( x ) et GetRMS(2)  RMS sur l'axe 2 ( y )

        if (align)
        {
          if ( clusterInStrip[X] && nEntries[X] < MAXSTRIPENTRY )
          {
            strip[X][nEntries[X]] = (int)xStripNum;
            pad[X][nEntries[X]] = (Int_t)(pad_x/PADWIDTH);
            nEntries[X]++;
          }

          if ( clusterInStrip[Y] && nEntries[Y] < MAXSTRIPENTRY )
          {
            strip[Y][nEntries[Y]] = (int)yStripNum;
            pad[Y][nEntries[Y]] = (Int_t)(pad_y/PADWIDTH);
            nEntries[Y]++;
          }
          //Extrapolate track to m2 prototype

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

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

            if (chbid == 10)
            {

              t2 = channel->GetBcIdDif();
              t3 = channel->GetBcIdHit();
              dt = t2 - t3;

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

              if (dt>=dt_min && dt<=dt_max)
              {
                hx_cor->Fill(pad_x,xpos);
                hy_cor->Fill(pad_y,ypos);

                hx_res->Fill(xpos - pad_x);
                hy_res->Fill(ypos - pad_y);

              } // end if (dt>=dt_min && dt<=dt_max)
            } // end if (chbid == 10)

          } // end (int j=0;j<nchannel;j++)

        } // end if align
      } // end if (nsuper_hot_tel == 3)

      if ( i != nEvent-2 )
      {
         hxy_tel_pad_temp->Reset(); 
      }
    } // end for (int i=0;i<nEvent-1;i++)
  }  // end while TKey


  TCanvas c1;
  c1.Divide(2,3);
  c1.cd(1);
  cout << " n entries X "  << nEntries[X] << endl;
  TGraph *gx = new TGraph(nEntries[X],strip[X],pad[X]);
  gx->SetTitle ( "Correlation pad - strip X" ); 
  gx->Draw("A*");
  gx->Fit("pol1");
 
  cout << " n entries Y "  << nEntries[Y] << endl;
  TGraph *gy = new TGraph(nEntries[Y],strip[Y],pad[Y]);
  gy->SetTitle ( "Correlation pad - strip Y" ); 
  c1.cd(2);
  gy->Draw("A*");
  gy->Fit("pol1");

  c1.cd(3);
  hx_cor->Draw();
  hx_cor->ProfileX()->Fit("pol1");
  
  c1.cd(4);
  hy_cor->Draw();;
  hy_cor->ProfileX()->Fit("pol1");

  c1.cd(5);
  hx_res->Draw();

  c1.cd(6);
  hy_res->Draw();


  theApp->Run();
  delete theApp;

  return 0;
}

---