/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/TB2012/Testbeam_april_2012/Prototype_m3/Beam/makeHistosM3.cpp

Go to the documentation of this file.

#include "tools/Log.hh"
#include "tools/Toolbox.hh"

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

#include "TFile.h"
#include "TTree.h"
#include "TKey.h"
#include "TBranch.h"
#include "TRandom.h"
#include "TProfile.h"
#include "TH2I.h"
#include "TF1.h"


#include <sstream>
#include <sys/types.h>
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>



using namespace std;





int main(int argc, char**argv){

  if ( argc != 3  )
    {

      cout<<endl;
      cout<<"**********************************************************************"<<endl;
      cout<<"Make histograms of hit position and time distribution for 3 thresholds"<<endl;
      cout<<"**********************************************************************"<<endl;
      cout<<endl;
      cout<<"      Two arguments: PATH and ROOTFILE"<<endl;
      cout<<"      e.g. ->  ./bin/root/makePlotsM3 PATH ROOTFILE"<<endl;
      cout<<endl;
      cout<<"**********************************************************************"<<endl;
      cout<<endl;
    }

  else
    {

      TString directory = argv[1];
      TString filename = argv[2];

      TString infile = directory + filename;

      //convert boardid into slab number
      int slab[6] = {2,0,0,1,1,2};

      TString time,name,title;
      int xpos,ypos,digit,nchannel,chbid,hardid,channelid,boardid,chipid;
      Long64_t t,t0,t1,t2,t3,dt,tabs;

      TH2I * hxy[50][3];
      TH1I * hc_nhit[50][3];
      TH1I * hdt[50][3];

      TH1I * hchb = new TH1I("hcbh",";chamber id",60,0,60);
      TH1I * hdt_all = new TH1I("hdt_all","Time to readout in all chambers;time (bcid)",100000,0,100000);

      for (int i=0;i<50;i++)
        {
          for (int k=0;k<3;k++)
            {
              name = Form("hc_nhit_%i_dac%i",i+1,k);
              title = Form("Hit per channel in chb %i - thr %i;channel",i+1,k);
              hc_nhit[i][k] = new TH1I(name,title,9216,0,9216);

              name = Form("hxy_%i_dac%i",i+1,k);
              title = Form("Hit distribution in chb %i - thr %i;channel",i+1,k);
              hxy[i][k] = new TH2I(name,title,96,0,96,96,0,96);

              name = Form("hdt_%i_dac_%i",i+1,k);
              title = Form("Time to readout in chb %i - thr %i;time (bcid)",i+1,k);
              hdt[i][k] = new TH1I(name,title,1e4,0,1e6);
            }
        }






      TFile * tfile_input = new TFile(infile);
      cout<<endl;
      cout<<"Reading file : "<<infile<<endl;
      cout<<endl;

      TIter nextkey(tfile_input->GetListOfKeys());
      TKey *key;

      int ntree = 0;



      while (key = (TKey*)nextkey())
        {
          TTree *tree = (TTree*)key->ReadObj();                
          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);

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

              //cout<<"Event "<<i<<" contains "<<nchannel<<" hits"<<endl;

              /*Look in prototype*/
              for (int j=0;j<nchannel;j++)
                {
                  MTChannel* channel = (MTChannel*)evt->GetChannels()->UncheckedAt(j);
                  chbid = channel->GetChamberId();

                  hchb->Fill(chbid);

                  xpos = channel->GetX()*1e-4 + 48;
                  ypos = channel->GetY()*1e-4 + 48;

                  digit = channel->GetDigitalValue();
                  hardid = channel->GetHardId();
                  chipid = channel->GetChipId();

                  boardid = channel->GetBoardId()%6;
                  boardid = slab[boardid];

                  channelid = (chipid-1)*64 + hardid + 3072*boardid;

                  t2 = channel->GetBcId_Dif();
                  t3 = channel->GetBcId_Hit();
                  dt = t2 - t3;

                  if (digit==3)
                    {
                      for (int k=0;k<3;k++)
                        {
                          hc_nhit[chbid-1][k]->Fill(channelid);
                          hxy[chbid-1][k]->Fill(xpos,ypos);
                          hdt[chbid-1][k]->Fill(dt);
                        }
                    }
                  else
                    {
                      if (digit==2)
                        {
                          for (int k=0;k<2;k++)
                            {
                              hc_nhit[chbid-1][k]->Fill(channelid);
                              hxy[chbid-1][k]->Fill(xpos,ypos);
                              hdt[chbid-1][k]->Fill(dt);
                            }
                        }
                      else
                        {
                          if (digit==1)
                            {
                              hc_nhit[chbid-1][0]->Fill(channelid);
                              hxy[chbid-1][0]->Fill(xpos,ypos);
                              hdt[chbid-1][0]->Fill(dt);
                            }
                        }
                    }
                  hdt_all->Fill(dt);
                }
            }
        }



      TString outfile = infile;

      for (int k=0;k<5;k++){outfile.Chop();}
      
      outfile += "_histos.root";

      TFile * tfile_output = new TFile(outfile,"RECREATE");

      cout<<endl;
      cout<<"Writing file : "<<outfile<<endl;
      cout<<endl;


      for (int i=0;i<50;i++)
        {
          for (int k=0;k<3;k++)
            {
              hc_nhit[i][k]->Write();
              hxy[i][k]->Write();
              hdt[i][k]->Write();
            }
        }
      




    }
}

---