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

Go to the documentation of this file.

/* @version $Revision: 1376 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-11-18 16:11:02 +0100 (Fri, 18 Nov 2011) $ */
#include <map>
#include "TKey.h"
#include <TH2F.h>
#include <TROOT.h>
#include <TRint.h>
#include <TKey.h>
#include <TList.h>
#include <TTree.h>
#include <TFile.h>
#include <stdlib.h>
#include <TSystem.h>
#include <iostream>
#include <sstream>
#include "Log.hh"

#include "MicroException.hh"
#include <string>

#include "mysql/Mysql.hh"


#include "root/CaloRun.hh"
#include "root/CaloEvent.hh"
#include "root/CaloHit.hh"
#include "root/MTRun.hh"
#include "root/MTEvent.hh"
#include "root/MTChannel.hh"

using namespace std;


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

   if ( argc !=2  ) {
   FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
   exit(1);
  }

  Mysql mysql;



   int search = 1;
   TH2F * hxy = new TH2F("hxy","",100,0,100,100,0,100);
   TH1I * hnhit = new TH1I("hnhit","Number of hits",100,0,100);
   TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in region of interest",100,0,100);
   TH1F * hxres = new TH1F("hxres","Residual in x direction;xm2 - xtel",201,-100,100);
   TH1F * hyres = new TH1F("hyres","Residual in x direction;xm2 - xtel",201,-100,100);

   string rootName;
   rootName.assign(argv[1]);
   int nbHit = 0;
   int nbEvt = 0;
  
  TFile f(rootName.c_str());
  TIter nextkey(f.GetListOfKeys());
  TKey *key;

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

    TString keyname = key->GetName();
    if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
    else{

    TTree *tree = (TTree*)key->ReadObj();                
    cout << "---------------------------- NEW TTREE ---------------------------"<< endl;
    nbEvt =0;
    nbHit = 0;

    float xoffset = 0;
    float yoffset = 0;

    try {
      CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
      ////CaloRun *run=(CaloRun *)tree->GetUserInfo()->First();

      //xoffset = run->GetOffsetX();
      //yoffset = run->GetOffsetY();

      //cout << "RunId[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; 
      string beam;
      try {
//        ui32 runId = mysql.getRunIdFromDhcalId(715573);
        beam  = mysql.getRunEnergy(run->GetRunId());
      }
      catch (MicroException e) {}
  
      ////cout << "Beam[" << beam << "]" << endl;
      //cout << "Offset telescope - m2[" << xoffset <<  ", "  << yoffset  << "] " << endl ; 
    }
    catch ( ... ) {}

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

    CaloHit* hit =NULL;
    int nbEntries = tree->GetEntries();
    //cout << "**** Event nb entries " << dec << nbEntries <<  endl;
    bool micromegasHit = false;
    for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
    {
      tree->GetEntry(evtNum);
      UInt_t nbRPCHit=0;
      nbEvt++;
      //cout << "**** Event Information "<< dec << evtNum+1 << "/" << nbEntries<< " ****" <<  endl;
      //cout << "    event num in TTree:" << evtNum <<endl;
      //cout << "    event id:" << evt->GetEventId() <<endl;

      //cout << "    number of channel hit:" << evt->GetNHits() <<"for event[" << evt->GetEventId() << "]" <<   endl;
      int nbHits = evt->GetNHits();
      //cout << "nb hits:" << nbHits << endl;

      float xtel = 0;// evt->GetXTelescope();
      float ytel =  0;//evt->GetYTelescope();

      float xextrapol = xtel - xoffset;
      float yextrapol = ytel - yoffset;

      hnhit->Fill(nbHits);

      //number of hits in region of interest
      int nhit_roi = 0;
      float xpos = 0;
      float ypos = 0;
      float xresidual = 0;
      float yresidual = 0;

      for(int i=0 ; i < nbHits ; i++)
      {
        hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
              xpos = hit->GetX();
              ypos = hit->GetY();
              hxy->Fill(xpos,ypos);
        if ( hit->GetLayer() >= 50)
        {
          micromegasHit = true;
        }
        if ( hit->GetLayer() < 50 ) { nbRPCHit++ ;} 


              xresidual = xpos - xextrapol;
              yresidual = ypos - yextrapol;
      
              hxres->Fill(xresidual);
              hyres->Fill(yresidual);

        if ((fabs(xresidual)<=search) && (fabs(yresidual)<=search))
          {
            nhit_roi++;
          }
      } // end hit

      hnhit_roi->Fill(nhit_roi);

    if ( evt->GetNHits()> 130 &&  micromegasHit && nbRPCHit < 100 ) 
    {
       cout << "    number of channel hit:" << evt->GetNHits() <<"for event[" << evt->GetEventId() << "] RPCHIT[" <<  nbRPCHit << "]" <<  endl;
    } 
    } // end event
  } // end run ( TTree )
}
  TFile tf("temp.root","RECREATE");
  hxres->Write();
  hyres->Write();
  hxy->Write();
  hnhit->Write();
  hnhit_roi->Write();

}

---