/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/yannis.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 <TProfile.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"

#include "TGraphErrors.h"
#include <TF1.h>

using namespace std;
const int maxEvents = 10000; // How many events to read
const int maxTrees = 2;  // How many trees to read
const int maxHits = 200; // Max nb of considered hits


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

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

   int search = 1;
   TH2F * hxy = new TH2F("hxy","",100,0,100,100,0,100);
   TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000);
   TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,50.);
   TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",130,0,130);
   TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50);
   TH2F * hxy33 = new TH2F("hxy33","Chambre 34 Map",100,0,100,100,0,100);
   TH2F * hxy41 = new TH2F("hxy41","Chambre 49 Map",100,0,100,100,0,100);
   TH2F * hxy49 = new TH2F("hxy49","Chambre 49 Map",100,0,100,100,0,100);
   TProfile * tp_nhit = new TProfile("tp_nhit","",60,0,60,0,200,"");
   TH1F * hChi2X = new TH1F("hChi2X"," Chi2 X ", 200, 0., 10000.);
   TH1F * hChi2Y = new TH1F("hChi2Y"," Chi2 Y ", 200, 0., 10000.);
   TH1F * hResidY = new TH1F("hResidY"," Residual Y ", 200, -100., 100.);
   TH1F * hResidX = new TH1F("hResidX"," Residual X ", 200, -100., 100.);
   TProfile * pResidX = new TProfile("pResidX","Average Residual per layer X",50,0,50,-1000.,1000,"");
   TProfile * pResidY = new TProfile("pResidY","Average Residual per layer Y",50,0,50,-1000.,1000,"");
   TH1F * hPenteY = new TH1F("hPentY"," Slope Y ", 200, -0.1, 0.1);
   TH1F * hPenteX = new TH1F("hPentX"," Slope X ", 200, -0.1, 0.1);
   TH2F * hxyLayer0 = new TH2F("hxyLayer0","X verus Y fitted track. Layer 0",100,0,1000,100,0,1000);

   TH1I * hz = new TH1I("hz","",100,0,100);

   string rootName;
   rootName.assign(argv[1]);
   int nbHit = 0;
   int nbEvt = 0;

   int NbHitLayer[100];
   bool IsPenetratingMIP;

   float Zlayer[maxHits];
   float Xlayer[maxHits];
   float Ylayer[maxHits];
   float ErPos[maxHits];
  
  TFile f(rootName.c_str());
  TIter nextkey(f.GetListOfKeys());
  TKey *key;

 TFile tf("debug_Pions8GeV.root","RECREATE");
 int NbTrees = 0;
  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;
    NbTrees++;
    if( NbTrees > maxTrees ) continue;
    CaloRun *run=NULL;

    try {
      CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
      if ( run !- NULL )
      {
        cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; 
      }
    }
    catch ( ... ) {}
    int RunId=1;
    TString name = Form("deltaTRun_%i",RunId);
    TString title = Form("delta distribution Run %i",RunId);
    TH1I * deltaTRun = new TH1I(name,title,20,-10,10);

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

    CaloHit* hit =NULL;
    int nbEntries = tree->GetEntries();
    
    for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
    {
      tree->GetEntry(evtNum);
      nbEvt++;
      if(nbEvt>maxEvents) continue;
      /*
      cout << "**** Event Information "<< dec << evtNum+1 << "/" << nbEntries<< " ****" <<  endl;
      cout << "    event num in TTree:" << evtNum <<endl;
      cout << "    event id:" << evt->fEventId <<endl;

      cout << "    number of channel hit:" << evt->GetNHits() <<endl;
      */
      int nbHits = evt->GetNHits();
      //      cout << "nb hits:" << nbHits << endl;

      hnhit->Fill(nbHits);


      if (nbHits > 120) continue;


      //number of hits in region of interest
      int nhit_roi = 0;
      int nhit_fit=0;
      float xpos = 0;
      float ypos = 0;
      float zpos=0;
      int deltat=0;

      for (int layer=0 ; layer <100 ; layer++){
        NbHitLayer[layer]=0;}
      for (int i=0 ; i<maxHits ;i++){
        Zlayer[i]=0;
        Xlayer[i]=0;
        Ylayer[i]=0;
        ErPos[i] = 30./sqrt(12.);
      }

      //Look for the Max Deltat
      //Find the 1st deltat and Book accordingly a histogram in order to find the maximum deltaT
        hit = (CaloHit*)evt->GetHits()->UncheckedAt(1);
        int deltat1 = hit->GetDeltaT();
        TH1I * hEvtDeltat = new TH1I("hEvtDeltat","Event DeltaT",20,deltat1-10,deltat1+10);

      for (int i=0 ; i <nbHits ; i++){
        hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
        deltat=hit->GetDeltaT();
        hEvtDeltat->Fill(deltat);
      }
      int maxDeltat = deltat1-10 + hEvtDeltat->GetMaximumBin()-1;
      //      cout << "Max DeltaT " << maxDeltat << endl;
      //     hEvtDeltat->Write();
            delete hEvtDeltat;

      for(int i=0 ; i < nbHits ; i++)
      {
        hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);

        deltat=hit->GetDeltaT();
        deltaTRun->Fill(maxDeltat-deltat);
        if(abs(deltat-maxDeltat)>2) continue;

        xpos = hit->GetX()/10000+48;
        ypos = hit->GetY()/10000+48;

        zpos = hit->GetZ();
        if (zpos == 1370000){zpos = 1372000;}
        zpos = zpos / 28000;

        hz->Fill(zpos);

        //cout  << zpos << " " << " Y " << ypos <<endl;

        NbHitLayer[(int)zpos]++; 

        hxy->Fill(xpos,ypos);
        if((int)zpos == 33 ){hxy33->Fill(xpos,ypos);}
        if((int)zpos == 41 ){hxy41->Fill(xpos,ypos);}
        if((int)zpos == 49 ){hxy49->Fill(xpos,ypos);}
        //      cout << "deltat "<< deltat << endl;       

        nhit_roi++;
        // Prepare for fitting a muon
        // Choose layers to make the fit
        if (nhit_fit < maxHits){
          if(NbHitLayer[(int)zpos] >1 ) continue ; // keep only 1st hit in layers
          Zlayer[nhit_fit] = zpos*28; // Convert coordinates to mm
          Xlayer[nhit_fit] = xpos*10.;
          Ylayer[nhit_fit] = ypos*10.;
          nhit_fit++;
        }
          
      } // end hit
      hnhit_roi->Fill(nhit_roi);

      // Compute Hit Density to find MIP events
      int NbLayers=0;
      for (int layer=0 ; layer < 50 ; layer++){
        if(NbHitLayer[layer]>0)
          {
            NbLayers++;
            tp_nhit->Fill(layer+1,NbHitLayer[layer]);     
          }
      }

      float HitDensity = 1. * nhit_roi/NbLayers;

      //          cout << "Event num " << nbEvt << " Nb of Hits " << nbHits << " nb Layers " << NbLayers <<  " Density " << HitDensity << endl; 

      hDensity->Fill(HitDensity);

      //Look for penetrating MIPS

      int NbLayersForw=0;
      int NbLayersBack=0;

      if(1 < HitDensity < 3){

        for (int layer=0 ; layer < 10 ; layer++){
          //cout<<layer<<" "<<NbHitLayer[layer]<<endl;
          if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;}
        }
        for (int layer=35 ; layer < 45 ; layer++){
          //cout<<layer<<" "<<NbHitLayer[layer]<<endl;
          if(0 < NbHitLayer[layer]  && NbHitLayer[layer] <= 3){NbLayersBack++;}
        }

        //cout<<NbLayersForw<<" "<<NbLayersBack<<endl;

        hnHitForw->Fill(NbLayersForw);
        hnHitForw->Fill(NbLayersBack+20);

        //Define penetrating MIP
        IsPenetratingMIP = 0;
        if(NbLayersForw>6 && NbLayersBack > 6){ IsPenetratingMIP = 1;}
      } // end Penetrating MIPS


      // Fit a straight line for penetrating MIPs 
       if(IsPenetratingMIP){
         //      for (int i=0;i<nhit_fit;i++){cout << i << " Z " << Zlayer[i] << " Y "<< Ylayer[i] << " X " << Xlayer[i] << endl;}
         float zmax = 50.* 28.;
        TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos);
        TF1 *f1 = new TF1("f1","pol1",0,zmax);
        grx->Fit("f1","Q");
        float betaX = f1->GetParameter(0);
        float alfaX = f1->GetParameter(1);
        double chi2X = f1->GetChisquare();
        TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos);
        TF1 *f2 = new TF1("f2","pol1",0,zmax);
        gry->Fit("f2","Q");
        float betaY = f2->GetParameter(0);
        float alfaY = f2->GetParameter(1);
        double chi2Y = f2->GetChisquare();
        hChi2X->Fill(chi2X);
        hChi2Y->Fill(chi2Y);
        hxyLayer0->Fill(betaX,betaY); // Make X Y Map at the first layer

        //              cout << " Pente X" << alfaX << "  offset " << betaX << " chi2 " << chi2X << endl;
        //              cout << " Pente Y" << alfaY << "  offset " << betaY << " chi2 " << chi2Y << endl;

        // Compute residuals to fitted track
        for (int j=1; j<nhit_fit ; j++){
          float ResX = Xlayer[j]-( alfaX*Zlayer[j]+betaX );
          float ResY = Ylayer[j]-( alfaY*Zlayer[j]+betaY );
          hResidY->Fill(ResY);
          hResidX->Fill(ResX);
          pResidX->Fill(Zlayer[j]/28,ResX);
          pResidY->Fill(Zlayer[j]/28,ResY);
        }
      }
 
      if (nbEvt%1000 == 0) {
      cout << "Processing event: " << nbEvt
           << endl;  
      }  
    } // end event

    deltaTRun->Write();

    

    } // end run ( TTree )
  
    
}
  tf.cd();
  tp_nhit->Write();
  hxy->Write();
  hxy33->Write();
  hxy41->Write();
  hxy49->Write();
  hz->Write();
  hnhit->Write();
  hnhit_roi->Write();
  hDensity->Write();
  hnHitForw->Write();
  hChi2X->Write();
  hChi2Y->Write();
  hResidY->Write();
  hResidX->Write();
  pResidX->Write();
  pResidY->Write();
  hxyLayer0->Write();
}

---