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

Go to the documentation of this file.

//////////////////////////////////////////////
//////////////////////////////////////////////
// Author : D. Girard, 07/2012
/////////////////////////////////////////////
/////////////////////////////////////////////

//includes and constants
#include <map>
#include "TKey.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 <TCanvas.h>
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include <TH2F.h>

using namespace std;

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

  if ( argc !=2  ) {
    FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
    exit(1);
  }
  
  //User's parameters
  string rootName;
  rootName.assign(argv[1]);

  // example of data file
  // /lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/DHCAL_714405_I*_0.slcio.root 

  char nameofplot[1000];
  sprintf(nameofplot,"./ShowerProfile_714405_DetLayer48.gif");
  
  //Plotting window
  TCanvas* c1 = new TCanvas("c1","Shower Profile", 200, 0, 1200, 1200);
  
  //Histogramms
  TH2F* h02DProfx = new TH2F("02DProfx", "Basic Energy 2D Profile;Layer (z);X", 50, 0., 50.,100,0.,100.); 
  TH2F* h02DProfy = new TH2F("02DProfy", "Basic Energy 2D Profile;Layer (z);Y", 50, 0., 50.,100,0.,100.); 
  TH1F* h0EProf = new TH1F("0EProf", "Basic Energy Profile;Layer (z); Hits",50, 0., 50.);
  TH2F* h2DProfx = new TH2F("2DProfx", "Normal Energy 2D Profile;Layer (z);X", 50, 0., 50.,100,0.,100.);  
  TH2F* h2DProfy = new TH2F("2DProfy", "Normal Energy 2D Profile;Layer (z);Y", 50, 0., 50.,100,0.,100.);  
  TH1F* hEProf = new TH1F("EProf", "Normal Energy Profile;Layer (z); Hits",50, 0., 50.);
  TH2F* hKs2DProfx = new TH2F("Ks2DProfx", "Micromegas Energy 2D Profile (detector in layer 48);Reconsructed layer (z);X", 50, 0., 50.,100,0.,100.);  
  TH2F* hKs2DProfy = new TH2F("Ks2DProfy", "Micromegas Energy 2D Profile (detector in layer 48);Reconsructed layer (z);Y", 50, 0., 50.,100,0.,100.);  
  TH1F* hKsEProf = new TH1F("KsEProf", "Micromegas Energy Profile (detector in layer 48);Reconstructed layer (z); Hits",50, 0., 50.);
 
  //Variables

  int nbEvt = 0;
  int nbEvt_tot_used = 0;
  int nbEvt_ShowerStart[50];
  for(int layer=0;layer<50;layer++) nbEvt_ShowerStart[layer]=0;
  int nbEvt_Micromegas = 0;

  //first data reading structure
  int NbTrees = 0;
  TChain chain("T");
  for (int i =0; i<5 ; i++){
    TString SubFile = Form("_I%i",i);
    TString fname = Form("/lapp_data/LC/Detecteurs/MicroMegas/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+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;
      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;
              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 ( ... ) {}
              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;

              // A first loop over events to get the number of evts with a shower starting at each layer
              for (int evtNum = 0; evtNum < nbEntries; evtNum++){

                tree->GetEntry(evtNum);
        
                //if(evt->IsSquare(400,int asu,int abs_time, int chamber=48)) continue;
                //if(evt->IsBorder(80,int asu,int abs_time, int chamber=48)) continue;
                //if(evt->Ischip(40,int chip,int abs_time, int chamber=48)) continue;

                int nHits = evt->GetNHits();
                double HitX[nHits];
                double HitY[nHits];
                int zSlot[nHits];
                for(int i=0 ; i < nHits ; i++)
                  {
                      hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
                      HitX[i] =  hit->GetX() - 5000 ;
                      HitY[i] =  hit->GetY() - 5000 ;
                      zSlot[i] = hit->GetLayer();
                  }
                
                //Count the hits of the main shower (near the center) at each layer
                int MainShowerHitperLayer[50];
                for (int layer = 0; layer < 50; layer++) MainShowerHitperLayer[layer] = 0;
                for (int hitNum = 0; hitNum < nHits; hitNum++) if(HitX[hitNum] > 400000 && HitX[hitNum] < 600000 && HitY[hitNum] > 400000 && HitY[hitNum] < 600000) MainShowerHitperLayer[zSlot[hitNum]]++;
                
                //Find shower start of this event
                int ShowerStart = -1;
                for (int layer = 0 ; layer < 50 ; layer++){
                  if (MainShowerHitperLayer[layer] >= 2 && MainShowerHitperLayer[layer+1] >= 4 ){
                    ShowerStart=layer;
                    break;
                  }
                }
                
                //Remember the number of evts with a shower starting at each layer
                for(int layer=0;layer<50;layer++) if(ShowerStart == layer) nbEvt_ShowerStart[layer]++; 
 
                //Remember the number of evts with a shower starting in the calo
                if(ShowerStart>0) nbEvt_tot_used++;

                //Remember the number of evts with a shower starting in the "micromegas" interval
                if(ShowerStart<=48 && ShowerStart>0) nbEvt_Micromegas++;
                
              }//end first loop over events
            }//end else
        }//end key
    }//end chain

  cout<<endl;
  
  //Main data reading structure
  NbTrees = 0;
  TObjArray *fileElement=chain.GetListOfFiles();
  TIter nex(fileElement);
  TChainElement *chElt=0;
  NchainFiles=0;
  while (( chElt=(TChainElement*)nex() )) 
    {
      TFile f(chElt->GetTitle());
      TIter nextkey(f.GetListOfKeys());
      TKey *key;
      NchainFiles++;
      while (key = (TKey*)nextkey()) 
        {
          TString keyname = key->GetName();
          if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
          else
            {
              TTree *tree = (TTree*)key->ReadObj();                
              NbTrees++;
              CaloRun *run=NULL;

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

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

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

              // Main loop over events
              int ev=0;
              for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
                {
                  
                  nbEvt++;
                  tree->GetEntry(evtNum);
                  
                  //if(evt->IsSquare(800,int asu,int abs_time, int chamber=48)) continue;
                  //if(evt->IsBorder(80,int asu,int abs_time, int chamber=48)) continue;
                  //if(evt->Ischip(40,int chip,int abs_time, int chamber=48)) continue;
 
                  int nHits = evt->GetNHits();
                  double HitX[nHits];
                  double HitY[nHits];
                  int zSlot[nHits];
                  //double HitEng[nHits];
                  for(int i=0 ; i < nHits ; i++){
                    hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
                    HitX[i] =  hit->GetX() - 5000 ;
                    HitY[i] =  hit->GetY() - 5000 ;
                    zSlot[i] = hit->GetLayer();
                    //HitEng[i] = hit->GetThreshold();
                    //cout<<"i "<<i<<"; X "<<HitX[i]<<"; Y "<<HitY[i]<<"; z "<<zSlot[i]<<"; Z "<<HitX[i]<<"; E "<<HitEng[i]<<endl;
                  }
                  
                  //Find the number of hits at each layer
                  int MainShowerHitperLayer[50];
                  int TotalHitperLayer[50];
                  for (int layer = 0; layer < 50; layer++) {
                    MainShowerHitperLayer[layer] = 0;
                    TotalHitperLayer[layer] = 0;
                  }
                  
                  for (int hitNum = 0; hitNum < nHits; hitNum++) {
                    //Basic
                    TotalHitperLayer[zSlot[hitNum]]++;
                    
                    //Count the hits of the main shower (near the center) at each layer
                    if(HitX[hitNum] > 400000 && HitX[hitNum] < 600000 && HitY[hitNum] > 400000 && HitY[hitNum] < 600000) {
                      MainShowerHitperLayer[zSlot[hitNum]]++;
                    }
                  }
                                  
                  //Find shower start of this event
                  int ShowerStart = -1;
                  for (int layer = 0 ; layer < 48 ; layer++){
                    if (MainShowerHitperLayer[layer] >= 2 && MainShowerHitperLayer[layer+1] >= 4 ){
                      ShowerStart=layer;
                      break;
                    }             

                  }
                  //cout<<"sh "<<ShowerStart<<endl;
                      
                  // Loop over hits
                  for (int hit = 0; hit < nHits; hit++) {    
                    int thelayer = zSlot[hit];
                    //float Edep = HitEng[hit];
                    double x = HitX[hit]/10000.;
                    double y = HitY[hit]/10000.;

                    //if(evtNum<10 && hit<10) cout<<"thelayer "<<thelayer<<"; x "<<x<<"; y "<<y<<"; ShowerStart "<<ShowerStart<<"; nbEntries "<<nbEntries<<"; nbEvt_ShowerStart[1] "<<nbEvt_ShowerStart[1]<<"; nbEvt_ShowerStart[ShowerStart]"<<nbEvt_ShowerStart[ShowerStart]<<endl; //<<"; Edep "<<Edep<<endl;
      
                    //Basic energy profile:
                    //All events cumulated
                    if(ShowerStart>0){
                      h02DProfx->Fill(thelayer, x, 100000./nbEvt_tot_used); // I normalize al histos to 100 000 evts
                      h02DProfy->Fill(thelayer, y, 100000./nbEvt_tot_used);
                      h0EProf->Fill(thelayer, 100000./nbEvt_tot_used);
                    }
                    
                    //Normal method energy profile:
                    //Events whith ShowerStart=1;
                    //We move our detector at every layer; 
                    //->we plot in bin 48 the energy deposited in layer 48 for ShoweStart=1,
                    //->we plot in bin 48-1 the energy deposited in layer 48-1 for ShowerStart=1,
                    //=>we plot in bin i the energy deposited in layer i, for all events with ShowerStart=1,
                    //to normalize we divide by the number of events that we have used.
                    if (ShowerStart == 1){  
                      h2DProfx->Fill(thelayer,x, 100000./nbEvt_ShowerStart[1]);
                      h2DProfy->Fill(thelayer,y, 100000./nbEvt_ShowerStart[1]);
                      hEProf->Fill(thelayer,100000./nbEvt_ShowerStart[1]);
                    } 
                    
                    //MicroMegas method energy profile:
                    //We are fixed in the 48th layer; 
                    //->we plot in bin 48 the energy deposited in layer 48 for ShowerStart=1, to normalize we divide by the number of events used: nbEvt_ShowerStart[layer1]
                    //->we plot in bin 48-1 the energy deposited in layer 48 for ShowerStart=2, to normalize we divide by the number of events used: nbEvt_ShowerStart[layer2]
                    //=>we plot in bin 48-ShowerStart+1 the energy deposited in layer 48, to normalize we divide by the number of events used: nbEvt_ShowerStart[ShowerStart]
                    if(thelayer==48 && ShowerStart>0 && ShowerStart<=48){
                      hKs2DProfx->Fill(48-ShowerStart+1,x, 100000./nbEvt_ShowerStart[ShowerStart]);
                      hKs2DProfy->Fill(48-ShowerStart+1,y, 100000./nbEvt_ShowerStart[ShowerStart]);
                      hKsEProf->Fill(48-ShowerStart+1,100000./nbEvt_ShowerStart[ShowerStart]);
                    }
                  } // End of the loop over hits    
                  
                  cout << "Processing event number " 
                       << ev << "\t\r" 
                       << flush; 
                  ev++;   
                  
                } // End of the loop over events
              cout << " Finished with events in this TREE " << endl;
            } // End else
        } // End Key
    } //End Chain
  
  //1 Plot for all events
  c1->Divide(3,3);  
  
  //Basic histo, all events
  c1->cd(1);          
  char Xtitle[100];
  sprintf(Xtitle,"All ");
  char buf[10];
  sprintf(buf,"%d",nbEvt_tot_used);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  h02DProfx->GetXaxis()->SetTitle(Xtitle);
  h02DProfx->Draw("Zcol");
  
  //Normal histo, all events with ShowerStart=1
  c1->cd(2);          
  //char Xtitle[100];
  sprintf(Xtitle,"All events with showers start at layer 1; Sum of ");
  //char buf[10];
  sprintf(buf,"%d",nbEvt_ShowerStart[1]);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  h2DProfx->GetXaxis()->SetTitle(Xtitle);
  h2DProfx->Draw("Zcol");
  
  //Micromegas histo, all events with ShowerStart<48
  c1->cd(3);          
  sprintf(Xtitle,"All events with ShowerStart<48; Sum of ");
  sprintf(buf,"%d",nbEvt_Micromegas);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  hKs2DProfx->GetXaxis()->SetTitle(Xtitle);
  hKs2DProfx->Draw("Zcol");
  
  //Basic histo, all events
  c1->cd(4);          
  sprintf(Xtitle,"All ");
  sprintf(buf,"%d",nbEvt_tot_used);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  h02DProfy->GetXaxis()->SetTitle(Xtitle);
  h02DProfy->Draw("Zcol");
  
  //Normal histo, all events with ShowerStart=1
  c1->cd(5);          
  sprintf(Xtitle,"All events with showers start at layer 1; Sum of ");
  sprintf(buf,"%d",nbEvt_ShowerStart[1]);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  h2DProfy->GetXaxis()->SetTitle(Xtitle);
  h2DProfy->Draw("Zcol");
  
  //Micromegas histo, all events with ShowerStart<48
  c1->cd(6);          
  sprintf(Xtitle,"All events with ShowerStart<48; Sum of ");
  sprintf(buf,"%d",nbEvt_Micromegas);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  hKs2DProfy->GetXaxis()->SetTitle(Xtitle);
  hKs2DProfy->Draw("Zcol");
  
  //Basic histo, all events
  c1->cd(7);          
  h0EProf->SetLineColor(1);
  sprintf(Xtitle,"All ");
  sprintf(buf,"%d",nbEvt_tot_used);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  h0EProf->GetXaxis()->SetTitle(Xtitle);
  h0EProf->Draw();
  
  //Normal histo, all events with ShowerStart=1
  c1->cd(8);          
  hEProf->SetLineColor(1);
  sprintf(Xtitle,"All events with showers start at layer 1; Sum of ");
  sprintf(buf,"%d",nbEvt_ShowerStart[1]);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  hEProf->GetXaxis()->SetTitle(Xtitle);
  hEProf->Draw();
  
  //Micromegas histo, all events with ShowerStart<48
  c1->cd(9);          
  hKsEProf->SetLineColor(1);
  sprintf(Xtitle,"All events with ShowerStart<48; Sum of ");
  sprintf(buf,"%d",nbEvt_Micromegas);
  strcat(Xtitle, buf);
  strcat(Xtitle," events.        Layer (z)");
  hKsEProf->GetXaxis()->SetTitle(Xtitle);
  hKsEProf->Draw();
  
  //save the plot
  c1->SaveAs(nameofplot);
      
  cout << " Going to finsh the work after" << nbEvt << " events "<< endl;
} 

---