/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/getChamberEfficiency.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 "TGraphAsymmErrors.h"
#include <TF1.h>
#include "tools/Arguments.hh"

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


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

  Arguments arg(argc,argv," getNbOfHits run1 [run2] max [run20] -p particleType");

  if ( arg.getNbOfArgs() < 1) { 
    arg.usage();
    exit(1);
  }
  // Analyse the command line to find what to do
  
  int NbOfRuns = arg.getNbOfArgs();

  string rootName[20];
  for (int nr=0 ; nr < NbOfRuns ; nr ++ ){rootName[nr].assign(arg.getArgument(nr+1)); }

  string  RootFileType = "Eff";  // With Efficiency and Multiplicity Corrections
  if ( arg.getOption("-r").size() != 0 )
    {
      RootFileType  = arg.getOption("-r");
    }
  string  ParticleType = "Muons";
  if ( arg.getOption("-p").size() != 0 )
    {
      ParticleType  = arg.getOption("-p");
    }
 
  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);
  TProfile * tp_nhit = new TProfile("tp_nhit","",60,0,60,0,200,"");
  TH1F * hxdir = new TH1F("hxdir"," x Angle between hit and all the next ", 200, -3., 3.);
  TH1F * hydir = new TH1F("hydir"," y Angle between hit and all the next ", 200, -3., 3.);
  TProfile * hLayerMult = new TProfile("hLayerMult","Average Multiplicity around track per layer Y",50,0,50, 0.,1000,"");
  TH1F * hChi2X = new TH1F("hChi2X"," Chi2 X ", 20, 0., 20.);
  TH1F * hChi2Y = new TH1F("hChi2Y"," Chi2 Y ", 20, 0., 20.);
  TH1F * hResidY = new TH1F("hResidY"," Residual Y ", 60, -30., 30.);
  TH1F * hResidX = new TH1F("hResidX"," Residual X ", 60, -30., 30.);
  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, -100.,100.);
  TH1F * hPenteX = new TH1F("hPentX"," Slope X ", 200, -100, 100.);
  TH1F * hOffsetX = new TH1F("hOffsetX"," Offset X ", 200, 0, 1000.);
  TH1F * hOffsetY = new TH1F("hOffsetY"," Offset Y ", 200, 0, 1000.);

  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);
  TH1I * hGlNbHtsFound  = new TH1I("hGlNbHtsFound"," Nb of Hits found in Chamber", 50, 0, 50);
  TH1I * hGlNbHtsExtrap  = new TH1I("hGlNbHtsExtrap"," Nb of Hits extrapolated in Chamber", 50, 0, 50);
  TH1I * hNbHitsFit = new TH1I("hNbHitsFit","Nb of Hits Used in Fit", 100, 0, 200 );

  // Book Layer Hit maps and efficiencies
  TH2I * hxyMap[50];
  TH2I * hTrackxyMap[50];
  TH2F * hEffMap[50];
  TH2I * hRegxyMap[50] ;
  TH2I * hxyEvtMap[50] ;

  TString name,title;
  int NbinsEffx=12;
  int NbinsEffy=12;
  int hbin=96;
  for (int i=0;i<50;i++)
    {
      name = Form("hxyMap_%i",i);
      title = Form("Hit distribution in Layer %i",i);
      hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin);
      name = Form("hTrackxyMap_%i",i);
      hTrackxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
      name = Form("hEffMap_%i",i);            
      hEffMap[i]= new TH2F(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
      name = Form("hRegxyMap_%i",i);
      hRegxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
    }

 
  int nbHit = 0;
  int nbEvt;

  int NbHitLayer[100];
  int IsPenetratingMIP;
  //  float Hits[2][maxHits];
  Double_t Zlayer[maxHits]; // Before it was Float
  Double_t Xlayer[maxHits];
  Double_t Ylayer[maxHits];
  Double_t ErPos[maxHits];
  TString OutFilename = "/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/"+ParticleType+rootName[0]+RootFileType+".root";  TFile tf(OutFilename,"RECREATE");
  cout << " Writting file : " << OutFilename << endl;

  int NbTrees = 0;
  int RunId=1;

  TChain chain("T");
  cout << " Reading Files : " << endl;
  for ( int nr=0; nr< NbOfRuns; nr++){
  for (int i =0; i<5 ; i++){
    //    if(i==2)continue;
    TString SubFile = Form("_I%i",i);
    TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName[nr]+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());
      //  TFile f(rootName.c_str());
      TIter nextkey(f.GetListOfKeys());
      TKey *key;
      NchainFiles++;
      cout << " Reading a new file in the chain " << NchainFiles << endl;
      while (key = (TKey*)nextkey()) 
        {

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

              TTree *tree = (TTree*)key->ReadObj();                
              cout << "---------------------------- NEW TTREE -----------------"<< endl;
              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 ( ... ) {}
              CaloEvent *evt =  new CaloEvent();
              TBranch *branch= tree->GetBranch("CaloEvent");
              branch->SetAddress(&evt);

              CaloHit* hit =NULL;
              int nbEntries = tree->GetEntries();
              cout << " Nb of events to read " << nbEntries << endl;
              TString name = Form("deltaTRun_%i",RunId);
              TString title = Form("delta distribution Run %i",RunId);
              TH1I * deltaTRun = new TH1I(name,title,20,-10,10);

              nbEvt = 0;

              for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
                {
                  tree->GetEntry(evtNum);
                  nbEvt++;
                  //              if(nbEvt%2000==0){cout << nbEvt << endl;}
                  if(nbEvt>maxEvents) break;
                  /*
                    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;
                  unsigned int ncol = 0;
                  unsigned int nrow = 0;
                  int zSlot = 0;
                  int deltat=0;
                  int maxDeltat = evt->GetDeltaTmax();

                  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;
                    double theta = (0.0136/100.)*sqrt((i+1)*2/1.757)*(1+0.038*log((i+1)*2/1.757));
                    ErPos[i] = sqrt(pow(10./sqrt(12.),2)+pow(2*28.*theta,2));
                  }

                  for (int i=0; i<50; i++) // Fill Maps for this event
                    {
                      name = Form("hxyEvtMap_%i",i);
                      title = Form("Hit distribution in Layer %i",i);
                      hxyEvtMap[i]= new TH2I(name,title,96,0,96,96,0,96);
                    }

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

                      xpos =  hit->GetX()/1000 ; // unit = mm
                      ypos =  hit->GetY()/1000 ;
                      zSlot = hit->GetLayer() ;
                      zpos =  hit->GetZ()/1000;
                      ncol = hit->GetColInChamber();
                      nrow = hit->GetRowInChamber();

                      if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs
                      if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas

                      //Eliminate Bad zones
                      //  if(zSlot==33 && xpos>64 ) continue;
                      // if(zSlot==49 && (xpos>=32 &&xpos<40) && ( ypos>=48 && ypos <56) ) continue;

                      hz->Fill(zSlot);

                      NbHitLayer[zSlot]++; 

                      hxy->Fill(ypos/10.,xpos/10.);

                      nhit_roi++; // Nb of good hits
                      hxyEvtMap[zSlot]->Fill(ncol,nrow,1);
                      hxyMap[zSlot]->Fill(ncol,nrow,1);

                      // Prepare for fitting a muon
                      // Choose layers to make the fit
                      if(nhit_fit==0){
                        Zlayer[nhit_fit] = zpos; // Convert coordinates to mm
                        Xlayer[nhit_fit] = xpos;
                        Ylayer[nhit_fit] = ypos;
                        nhit_fit++;
                      }
                      else{
                        if (nhit_fit < maxHits){
                          if(NbHitLayer[zSlot] >1 ) continue ; // keep only 1st hit in layers
                          // Select hits that have almost the same slope as the previous one            
                          Double_t slopy = TMath::ATan2( (ypos-Ylayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
                          Double_t slopx = TMath::ATan2( (xpos-Xlayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
                          hxdir->Fill(slopx); hydir->Fill(slopy);
                          //    hSlopyY->Fill(slopy);  hSlopx->Fill(slopx);
                          //  cout << slopx << " " << slopy << endl;
                          if ( TMath::Abs(slopx)> 5. || TMath::Abs(slopy)>5. ) continue;
                          Zlayer[nhit_fit] = zpos; 
                          Xlayer[nhit_fit] = xpos;
                          Ylayer[nhit_fit] = ypos;
                          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,NbHitLayer[layer]);   
                      }
                  }
                  float HitDensity = 1. * nhit_roi/NbLayers;

                  hDensity->Fill(HitDensity);
                 
                  //Look for penetrating MIPS

                  int NbLayersForw=0;
                  int NbLayersBack=0;

                  IsPenetratingMIP = 0;
                  if( HitDensity > 1 && HitDensity < 3){

                    for (int layer=0 ; layer < 10 ; layer++){
                      if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;}
                    }
                    for (int layer=38 ; layer < 48 ; layer++){
                      if(0 < NbHitLayer[layer]  && NbHitLayer[layer] <= 3){NbLayersBack++;}
                    }

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

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

                  // Fit a straight line for penetrating MIPs 
                  if(IsPenetratingMIP==1 ){ hNbHitsFit->Fill(nhit_fit);} // How many hits are used for the fit    
                  if(IsPenetratingMIP==1 &&  (nhit_fit > 25) ){
                    float zmax = 1372.;
                    TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos);
                    TF1 *f1 = new TF1("f1","pol1",0,zmax);
                    grx->Fit("f1","Q");
                    double betaX = f1->GetParameter(0);
                    double alfaX = f1->GetParameter(1);
                    double chi2X = f1->GetChisquare();
                    delete f1; delete grx;
                    TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos);
                    TF1 *f2 = new TF1("f2","pol1",0,zmax);
                    gry->Fit("f2","Q");
                    double betaY = f2->GetParameter(0);
                    double alfaY = f2->GetParameter(1);
                    double chi2Y = f2->GetChisquare();
                    // if(nbEvt<100) gry->Write();
                    delete f2; delete gry;
                    hPenteX->Fill(alfaX*1000.);
                    hPenteY->Fill(alfaY*1000.);
                    hOffsetX->Fill(betaX);
                    hOffsetY->Fill(betaY);
                    hChi2X->Fill(chi2X/(nhit_fit-2));
                    hChi2Y->Fill(chi2Y/(nhit_fit-2));
                    hxyLayer0->Fill(betaX,betaY); // Make X Y Map at the first layer

                    // 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);
                    }

                    // Go to compute Efficiency with good tracks
                    if( chi2X/(nhit_fit-2) < 6. && chi2Y/(nhit_fit-2) < 6 ) {
                    for(int i=0; i<50 ; i++){
                      int NpadsLook = 3;
                      float Xextrapol = alfaX*i*28. + betaX ; // X extrapolate to chamber i in mm
                      float Yextrapol = alfaY*i*28. + betaY ; // Y extrapolate to chamber i in mm

                      // Find Central pad from computed x,y position
                      Int_t xpadCentr = -1000;
                      Int_t ypadCentr = -1000;
                      if( i>47 && ( fabs(Yextrapol-322.5) < 20 || fabs(Yextrapol-650) < 20 ) ) continue;
                      if( i>47 && ( fabs(Xextrapol-485.0) < 20  ) ) continue;
                      if(i<48){
                        xpadCentr = (int)(Xextrapol/10.4);
                        ypadCentr = (int)(Yextrapol/10.4);
                      }
                      else{
                        if(Xextrapol<=480) {xpadCentr = (int)(Xextrapol/10.);}
                        else {xpadCentr = (int)( (Xextrapol-5)/10 );}
                        if(Yextrapol<=320.) {ypadCentr = (int)( (Yextrapol/10.) );}
                        if(Yextrapol>325 && Yextrapol<=645.) {ypadCentr = (int)( (Yextrapol-5.)/10.);}
                        if(Yextrapol>=646) {ypadCentr = (int)((Yextrapol-10.)/10.);}
                      }
                      // check that we are inside the chamber 
                      if(xpadCentr < 0  || xpadCentr > 95 ) continue;
                      if(ypadCentr < 0  || ypadCentr > 95 ) continue;
                      // define fiducial volume
                      if(xpadCentr <= NpadsLook || xpadCentr>95-NpadsLook) continue ;
                      if(ypadCentr <= NpadsLook || ypadCentr>95-NpadsLook) continue ;

                      hTrackxyMap[i]->Fill(ypadCentr,xpadCentr,1);  //Fill extrapolated point
                      hGlNbHtsExtrap->Fill(i);
                      //Check now if there is a hit around this point in a region of NpadsLook
                      int xbinL = max(xpadCentr - NpadsLook,0); // For the PS data we must increase the search zone in x
                      int xbinH = min(xpadCentr + NpadsLook,95);
                      int ybinL = max(ypadCentr - NpadsLook,0);
                      int ybinH = min(ypadCentr + NpadsLook,95);
                      int NhitsRegSearch = hxyEvtMap[i]->Integral(ybinL,ybinH,xbinL,xbinH);
                      if( NhitsRegSearch >0 ) {
                        hRegxyMap[i]->Fill(ypadCentr,xpadCentr,1);  //Fill extrapolated point
                        hGlNbHtsFound->Fill(i) ;
                        hLayerMult->Fill(i, NhitsRegSearch) ; }
                    
                    } // end of efficiency coputing end layer loop
                    } // end of Chi2 selection
                  } // end of penetrating MIP
 
                  if (nbEvt%10000 == 0) {
                    cout << "Processing event: " << nbEvt
                         << endl;  
                  }  
                  for (int i=0; i<50 ; i++){
                    delete hxyEvtMap[i];
                  } 
                  tf.cd(0);
                }  // end event 
              cout << " Finished with events in this TREE " << endl;
              deltaTRun->Write(); 
              delete deltaTRun;
            } // else 
        } // end of tree or run
      cout << " Ending  file in chain " << endl;
    } // end of chain
  cout << " Going to finsh the work after" << nbEvt << " events "<< endl;
  tf.cd();
  hPenteY->Write();
  hPenteX->Write();
  hOffsetX->Write();
  hOffsetY->Write();
  hxdir->Write();
  hydir->Write();
  tp_nhit->Write();
  hLayerMult->Write();
  hxy->Write();
  for ( int nch=0; nch<=49; nch++){hxyMap[nch]->Write();}
//   hxyMap[33]->Write();
//   hxyMap[31]->Write();
//   hxyMap[48]->Write();
//   hxyMap[49]->Write();
  hz->Write();
  hnhit->Write();
  hnhit_roi->Write();
  hDensity->Write();
  hnHitForw->Write();
  hNbHitsFit->Write();
  hChi2X->Write();
  hChi2Y->Write();
  hResidY->Write();
  hResidX->Write();
  pResidX->Write();
  pResidY->Write();
  hxyLayer0->Write();
  for (int i=0 ; i<50 ; i++){
    hRegxyMap[i]->Write();
    hTrackxyMap[i]->Write();
    hEffMap[i]->Divide(hRegxyMap[i], hTrackxyMap[i], 100., 1.,"b"); // Binomial errors
    hEffMap[i]->Write();
  }
  TH1F * hEffRPC = new TH1F("hEffRPC"," RPC Efficiencies for every chip ", 100, 0., 100.);
  TH1F * hEffmMegas = new TH1F("hEffmMegas"," mMegas Efficiencies for every chip ", 100, 0., 100.);
  for (int i = 0; i< 50; i++){
    for (int j = 1; j<13 ; j++){
      for (int k=1; k<13;  k++) {
        float eff = hEffMap[i]->GetBinContent(j,k);
        if(i<48) {  hEffRPC->Fill(eff);   }
        if(i>=48){  hEffmMegas->Fill(eff); }
      }// end columns
    } //end rows
  }// end Chamber

  hEffRPC->Write();
  hEffmMegas->Write();

  // Compute Global efficiency
  TH1F * hGlEffChamber = new TH1F("hGlEffChamber"," Chamber Global Efficiency ", 50, 0, 50 );
  hGlEffChamber->Divide(hGlNbHtsFound, hGlNbHtsExtrap,  100. , 1., "b"); // Binomial errors
  hGlNbHtsFound->Write();
  hGlNbHtsExtrap->Write();
  hGlEffChamber->Write();
  for (int j=1; j<=50 ; j++){
    double nbTrFound= hGlNbHtsFound->GetBinContent(j);
    double nbTrExtr=hGlNbHtsExtrap->GetBinContent(j);
    double ef=hGlEffChamber->GetBinContent(j);
    double erEff=hGlEffChamber->GetBinError(j);
    double w = nbTrFound/nbTrExtr;
    double binErr=w*(1-w)/nbTrExtr;
    cout << "chamber " << j << " Nb tracks found and extrapol " << nbTrFound << "  " << nbTrExtr << " Efficiency " << ef << " +- " << binErr << endl;    
  }
  TGraphAsymmErrors * gr = new TGraphAsymmErrors(hGlNbHtsFound, hGlNbHtsExtrap,"");
  gr->Write();   
} 

---