/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/residual.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 <TChain.h>
#include <TChainElement.h>
#include <TMath.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;


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

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


  int mip = 0;
  int penetrating_mip = 0;
  float ngoodchi2 = 0;
  float ntest48 = 0;
  float neff48 = 0;

  TString name,title;

  TF1 * fzx = new TF1("fzx","pol1",0,60);
  TF1 * fzy = new TF1("fzy","pol1",0,60);

  TH1I * hx_evt = new TH1I("hx_evt","Horizontal distribution;x (pad number)",1000000,0,1000000);
  TH1I * hy_evt = new TH1I("hy_evt","Vertical distribution;y (pad number)",1000000,0,1000000);

  TH1I * hx_all = new TH1I("hx_all","Horizontal distribution;x (pad number)",1000000,0,1000000);
  TH1I * hy_all = new TH1I("hy_all","Vertical distribution;y (pad number)",1000000,0,1000000);

  TH1F * heff = new TH1F("heff","Efficiency per layer",60,0,60);
  TH1F * hmult = new TH1F("hmult","Multiplicity per layer",60,0,60);

  TH1I * hnhitcalo = new TH1I("hnhitcalo","Number of hits - All calo",1000,0,1000);
  TH1I * hdtcalo = new TH1I("hdtcalo","Time to event time - All calo",41,-20,20);
  TH1F * hhit_dens = new TH1F("hhit_dens","Number of hit per layer",100,0,50);

  TH1F * hzx_slope = new TH1F("hzx_slope","horizontal slope",101,-3,3);
  TH1F * hzy_slope = new TH1F("hzy_slope","vertical slope",101,-3,3);

  TH1F * hzx_const = new TH1F("hzx_const","horizontal const",100,0,100);
  TH1F * hzy_const = new TH1F("hzy_const","vertical const",100,0,100);
  
  TH1F * hzx_chi2 = new TH1F("hzx_chi2","horizontal chi2/ndf",1000,0,100);
  TH1F * hzy_chi2 = new TH1F("hzy_chi2","vertical chi2/ndf",1000,0,100);

  TH2I * hxy[50];
  TH1I * hdt[50];
  TH1I * hnhit[50];
  TH1F * hresx[50];
  TH1F * hresy[50];
  int nhit[50];
  float neff[50];

  for (int i=0;i<50;i++)
    {
      nhit[i] = 0;

      name = Form("hdt_%i",i);
      title = Form("Time to event time - chb %i",i);
      hdt[i]= new TH1I(name,title,41,-20,20);

      name = Form("hresx_%i",i);
      title = Form("Horizontal residuals - chb %i",i);
      hresx[i]= new TH1F(name,title,101,-500000,500000);

      name = Form("hresy_%i",i);
      title = Form("Vertical residuals - chb %i",i);
      hresy[i]= new TH1F(name,title,101,-500000,500000);

      name = Form("hnhit_%i",i);
      title = Form("Number of CaloHit - chb %i",i);
      hnhit[i]= new TH1I(name,title,100,0,100);

      name = Form("hxy_%i",i);
      title = Form("CaloHit position - chb %i",i);
      hxy[i] = new TH2I(name,title,100,0.,1000000.,100,0.,1000000.);
                                   
      
    }






  string rootName;
  rootName.assign(argv[1]);

  TChain chain("T");
  for (int i=0 ; i<6 ; i++)
    {
      if ( i != 5 ) continue;
      TString SubFile = Form("_I%i",i);
      TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root");
      //TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/Calo_HCAL_714408"+SubFile+"_0.slcio.root");
      chain.Add(fname);
      cout << fname << endl;
    }

  TObjArray *fileElements=chain.GetListOfFiles();
  TIter next(fileElements);
  TChainElement *chEl=0;
  int NchainFiles=0;

  while (( chEl=(TChainElement*)next() ))
    {
      TFile f(chEl->GetTitle());
      TIter nextkey(f.GetListOfKeys());
      TKey *key;

      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;
            int nbEvt =0;
            int nbHit = 0;
            NbTrees++;
            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;
            CaloEvent *evt =  new CaloEvent();
            TBranch *branch= tree->GetBranch("CaloEvent");
            branch->SetAddress(&evt);

            CaloHit* hit =NULL;
            int nbEntries = tree->GetEntries();

         // nbEntries = 30;
          for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
          {
                TGraph * tgzx = new TGraph();
                TGraph * tgzy = new TGraph();

                tree->GetEntry(evtNum);
                nbEvt++;

                int nbHits = evt->GetNHits();
                int maxDeltat = evt->GetDeltaTmax();

                cout<<"evt "<<evtNum<<" -> nhit = "<<nbHits<<" \r"<<flush;

                hnhitcalo->Fill(nbHits);

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

                    int deltat = hit->GetDeltaT();
                    int xpos =  hit->GetX();// * 1e-4;
                    int ypos =  hit->GetY();// * 1e-4;
                    int zSlot = hit->GetLayer() ;
                    
                    hdtcalo->Fill(maxDeltat-deltat);
                    hdt[zSlot]->Fill(maxDeltat-deltat);
                    
                    if (zSlot >= 48)
                      {
                        if  (fabs(maxDeltat-deltat) <= 2)
                          {
                            nhit[zSlot]++;
                            hx_evt->Fill(xpos);
                            hy_evt->Fill(ypos);
                          }
                      }
                    else
                      {
                        if  (fabs(maxDeltat-deltat) <= 1)
                          {
                            nhit[zSlot]++;
                            hx_evt->Fill(xpos);
                            hy_evt->Fill(ypos);
                          }
                      }
                  }

                int bin_max_x = hx_evt->GetMaximumBin();
                int bin_max_y = hy_evt->GetMaximumBin();


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

                    int deltat = hit->GetDeltaT();
                    int xpos =  hit->GetX();// * 1e-4;
                    int ypos =  hit->GetY();// * 1e-4;
                    int zSlot = hit->GetLayer() ;
                    bool close = 0;
                    int close_cut = 5;
                    
                    if (zSlot >= 48)
                      {
                        //inverse DIF
                        //if (xpos<=32){xpos+=64;}
                        //else{if (xpos>64){xpos-=64;}}
                        
                        if  (fabs(maxDeltat-deltat) <= 2)
                          {
                            //hxy[zSlot]->Fill(xpos,ypos);
                            close =  (fabs(bin_max_x-1 - xpos) <= close_cut) && (fabs(bin_max_y-1 - ypos) <= close_cut);
                            
                            if (nhit[zSlot] <= 4 && close)
                              {
                                tgzx->SetPoint(tgzx->GetN(),zSlot+1,xpos);
                                tgzy->SetPoint(tgzy->GetN(),zSlot+1,ypos);
                              }                   
                          }
                      }
                    else
                      {
                        if  (fabs(maxDeltat-deltat) <= 1)
                          {
                            //hxy[zSlot]->Fill(xpos,ypos);
                            close =  (fabs(bin_max_x-1 - xpos) <= close_cut) && (fabs(bin_max_y-1 - ypos) <= close_cut);
                            
                            if (nhit[zSlot] <= 4 && close)
                              {
                                tgzx->SetPoint(tgzx->GetN(),zSlot+1,xpos);
                                tgzy->SetPoint(tgzy->GetN(),zSlot+1,ypos);
                              }
                          }
                      }
                  }
                
                int nhot_layer = 0;
                int nhit_tot = 0;
                int nfirst_ten = 0;
                int nlast_ten = 0;
                
                for (int chb=0;chb<50;chb++)
                  {
                    hnhit[chb]->Fill(nhit[chb]);
                    
                    if (nhit[chb]!=0)
                      {
                        if (chb<10){nfirst_ten++;}
                        if (chb>=35 && chb<45){nlast_ten++;}
                        
                        nhot_layer++;
                        nhit_tot += nhit[chb];
                      }
                  }
                

                float hit_density = 1. * nhit_tot / nhot_layer;
                hhit_dens->Fill(hit_density);

                
                if (hit_density <= 4)
                  {
                    mip++;
                    
                    if (nfirst_ten>=7 && nlast_ten>=7)
                      {
                        penetrating_mip++;

                        fzx->SetParameters(bin_max_x,0);
                        fzy->SetParameters(bin_max_y,0);
                        
                        tgzx->Fit(fzx,"q");
                        tgzy->Fit(fzy,"q");
                        
                        float chi2zx = fzx->GetChisquare()/fzx->GetNDF();
                        float chi2zy = fzy->GetChisquare()/fzy->GetNDF();
                        
                        hzx_chi2->Fill(chi2zx);
                        hzy_chi2->Fill(chi2zy);
                        
                        hzx_slope->Fill(fzx->GetParameter(1));
                        hzy_slope->Fill(fzy->GetParameter(1));
                        
                        hzx_const->Fill(fzx->GetParameter(0));
                        hzy_const->Fill(fzy->GetParameter(0));
                        
                        if (chi2zx<=1 && chi2zy<=1)
                          {
                            ngoodchi2++;
                            
                            for(int i=0 ; i < nbHits ; i++)
                              {
                                hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
                                
                                int deltat = hit->GetDeltaT();
                                int xpos =  hit->GetX();// * 1e-4;
                                int ypos =  hit->GetY();// * 1e-4;
                                int zSlot = hit->GetLayer() ;
                                
                                if (zSlot >= 48)
                                  {
                                    //inverse DIF
                                    //if (xpos<=32){xpos+=64;}
                                    //else{if (xpos>64){xpos-=64;}}
                                    
                                    if  (fabs(maxDeltat-deltat) <= 2)
                                      {
                                        if (nhit[zSlot]<=4)
                                          {
                                            hxy[zSlot]->Fill(xpos,ypos);

                                            hresx[zSlot]->Fill(fzx->Eval(zSlot+1) - xpos);
                                            hresy[zSlot]->Fill(fzy->Eval(zSlot+1) - ypos);
                                          }
                                      }
                                  }
                                else
                                  {
                                    if  (fabs(maxDeltat-deltat) <= 1)
                                      {
                                        if (nhit[zSlot]<=4)
                                          {
              //cout << "zSlot[" << zSlot << "]"  << endl;
              //cout << "xpos[" << xpos << "]"  << endl;
              //cout << "ypos[" << ypos << "]"  << endl;
              //cout << "hxy[zSlot]" << hxy[zSlot] << "]"  << endl;
                                            hxy[zSlot]->Fill(xpos,ypos);
                                            hresx[zSlot]->Fill(fzx->Eval(zSlot+1) - xpos);
                                            hresy[zSlot]->Fill(fzy->Eval(zSlot+1) - ypos);                              
                                          }
                                      }
                                  }
                              }
                            
                            for (int chb=0;chb<50;chb++)
                              {
                                if (nhit[chb]!=0)
                                  {
                                    neff[chb]++;
                                    hmult->Fill(chb,nhit[chb]);
                                  }
                              }
                            
                            if (nhit[49]!=0)
                              {
                                ntest48++;
                                if (nhit[48]!=0)
                                  {
                                    neff48++;
                                  }
                              }
                            
                            
                          }
                        
                        name = Form("tgzx_%i",evtNum);
                        title = Form("Hit in zx evt %i",evtNum);
                        tgzx->SetName(name);
                        tgzx->SetTitle(title);
                        
                        name = Form("tgzy_%i",evtNum);
                        title = Form("Hit in zy evt %i",evtNum);
                        tgzy->SetName(name);
                        tgzy->SetTitle(title);
                        
                        tgzx->SetMarkerStyle(20);
                        tgzy->SetMarkerStyle(20);
                        
                        
                        
                        for (int b=0;b<1000000;b++)
                          {
                            int content = hx_evt->GetBinContent(b+1);
                            hx_all->Fill(b+1 - bin_max_x,content);
                            
                            content = hy_evt->GetBinContent(b+1);
                            hy_all->Fill(b+1 - bin_max_y,content);
                          }
                      }
                  }
                
                for (int chb=0;chb<50;chb++)
                  {
                    nhit[chb] = 0;
                  }
                
                hx_evt->Reset();
                hy_evt->Reset();
                
                tgzx->Delete();
                tgzy->Delete();
                
              }
          }
        }
    }

  cout<<endl;
  cout<<"MIP = "<<mip<<endl;
  cout<<"penetrating = "<<penetrating_mip<<endl;
  cout<<"good chi2 = "<<ngoodchi2<<endl;
  cout<<endl;

  for (int chb=0;chb<50;chb++)
    {
      float eff = neff[chb]/(1.*ngoodchi2);
      float eff_err = eff*(1-eff)/sqrt(1.*ngoodchi2);

      heff->SetBinContent(chb+1,eff);
      heff->SetBinError(chb+1,eff_err);

      hmult->SetBinContent(chb+1,hmult->GetBinContent(chb+1)/neff[chb]);

      cout<<"Chb "<<chb+1<<" eff = "<<eff<<" +/- "<<eff_err<<" (Nevt = "<<ngoodchi2<<")"<<endl;
    }

  float eff48 = neff48/ntest48;
  float eff48_err = eff48 * (1-eff48) / sqrt(ntest48);

  cout<<endl;
  cout<<"Chamber 49 (uM#2) = "<<neff48<<" / "<<ntest48<<" = "<<eff48<<" +/- "<<eff48_err<<endl;
  cout<<endl;




  TFile tf("Muons100GeV_deadzone.root","RECREATE");

  hx_all->Write();
  hy_all->Write();

  heff->Write();
  hmult->Write();


  hzy_slope->Write();

  hzx_const->Write();
  hzy_const->Write();

  hzx_chi2->Write();
  hzy_chi2->Write();

  hhit_dens->Write();
  hdtcalo->Write();
  hnhitcalo->Write();

  for (int i=0;i<50;i++)
    {
      hdt[i]->Write();
      hxy[i]->Write();
      hnhit[i]->Write();
      hresx[i]->Write();
      hresy[i]->Write();
    }
}

---