Micromegas LAPP reconstruction

00001 /* @version $Revision: 1376 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-11-18 16:11:02 +0100 (Fri, 18 Nov 2011) $ */
00002 /* This program finds the Shower Start Yannis */
00003 #include <map>
00004 #include "TKey.h"
00005 #include <TH2F.h>
00006 #include <TProfile.h>
00007 #include <TROOT.h>
00008 #include <TRint.h>
00009 #include <TKey.h>
00010 #include <TList.h>
00011 #include <TTree.h>
00012 #include <TFile.h>
00013 #include <TChain.h>
00014 #include "TChainElement.h"
00015 #include <TMath.h>
00016 #include <stdlib.h>
00017 #include <TSystem.h>
00018 #include <iostream>
00019 #include <sstream>
00020 #include "Log.hh"
00021 
00022 
00023 #include "MicroException.hh"
00024 #include <string>
00025 
00026 #include "mysql/Mysql.hh"
00027 
00028 
00029 #include "root/CaloRun.hh"
00030 #include "root/CaloEvent.hh"
00031 #include "root/CaloHit.hh"
00032 #include "root/MTRun.hh"
00033 #include "root/MTEvent.hh"
00034 #include "root/MTChannel.hh"
00035 
00036 #include "TGraphErrors.h"
00037 #include "TGraphAsymmErrors.h"
00038 #include <TF1.h>
00039 
00040 using namespace std;
00041 const int maxEvents = 3000000; // How many events to read
00042 const int maxTrees = 100;  // How many trees to read
00043 const int maxHits = 200; // Max nb of considered hits
00044 
00045 
00046 int main(int argc, char**argv){
00047 
00048   if ( argc !=2  ) {
00049     FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
00050     exit(1);
00051   }
00052 
00053   int search = 1;
00054 
00055 
00056   TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,200.);
00057   TH2F * hxy = new TH2F("hxy","Hit Distribution",100,0,100.,100,0,100.);
00058   TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000);
00059   TH1I * hShowerStart = new TH1I("hShowerStart","Start of the Shower",50, 0, 50);
00060   TH1I * hShowerStartDa = new TH1I("hShowerStartDa","Start of the Shower",50, 0, 50);
00061   TH1I * hz = new TH1I("hz","Nb of Hits per Layer",50,0,50);
00062   TProfile * tp_nhit = new TProfile("tp_nhit","Average Nb of Hits / Layer",50,0,50,0,200,"");
00063   TProfile * hLongprof = new TProfile("hLongprof","Long Profile reconstructed from 4 layers",50,0,50,0,150,"");
00064   TProfile * hLongprof1 = new TProfile("hLongprof1","Long Profile from 1 layer",50,0,50,0,150,"");
00065   TProfile * hLongprofL1 = new TProfile("hLongprofL1","Long Profile when Shower Starts Layer 1",50,0,50,0,150,"");
00066   TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",200,0,2000.);
00067 
00068   TString name,title,name1,title1 ;
00069   int hbin=96;
00070   // Book Layer Hit maps
00071   TH2I * hxyMap[50];
00072   TH2I * hxyMapTot[50];
00073   for (int i=0;i<50;i++)
00074     {
00075       name = Form("hxyMap_%i",i);
00076       title = Form("Hit distribution in Layer %i",i);
00077       hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin);
00078       name1 = Form("hxyMapTot_%i",i);
00079       title1 = Form("All events Hit distribution in Layer %i",i);
00080       hxyMapTot[i]= new TH2I(name1,title1,96,0,hbin,96,0,hbin);
00081     }
00082   TH2I * hxyTot = new TH2I("hxyTot","All Layers of 1 Event",96,0,hbin,96,0,hbin);
00083 
00084   int NbHitLayer[100];
00085   int NbHitLayer77[100];
00086   float RG1G3 = -100.;
00087 
00088   string rootName;
00089   rootName.assign(argv[1]);
00090   int nbHit = 0;
00091   int nbEvt = 0;
00092   TString RootFiles = "/lapp_data/LC/Detecteurs/Sid/karyotak/DataTB2012/rootfiles/" ;
00093   TString OutFilename = RootFiles+"ShowerStart"+rootName+".root";
00094   TFile tf(OutFilename,"RECREATE");
00095   int NbTrees = 0;
00096   int RunId=1;
00097 
00098   TChain chain("T");
00099   for (int i =0; i<5 ; i++){
00100     TString SubFile = Form("_I%i",i);
00101     TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root");
00102     chain.Add(fname);
00103     cout << fname << endl;
00104   }
00105 
00106   TObjArray *fileElements=chain.GetListOfFiles();
00107   TIter next(fileElements);
00108   TChainElement *chEl=0;
00109   int NchainFiles=0;
00110   while (( chEl=(TChainElement*)next() )) 
00111     {
00112       TFile f(chEl->GetTitle());
00113       //  TFile f(rootName.c_str());
00114       TIter nextkey(f.GetListOfKeys());
00115       TKey *key;
00116       NchainFiles++;
00117       cout << " Reading a new file in the chain " << NchainFiles << endl;
00118       while (key = (TKey*)nextkey()) 
00119         {
00120 
00121           TString keyname = key->GetName();
00122           if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
00123           else
00124             {
00125 
00126               TTree *tree = (TTree*)key->ReadObj();                
00127               cout << "---------------------------- NEW TTREE -----------------"<< endl;
00128               nbHit = 0;
00129               NbTrees++;
00130               if( NbTrees > maxTrees ) continue;
00131               CaloRun *run=NULL;
00132 
00133               try {
00134                 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
00135                 if( run != NULL ){cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; }
00136 
00137               }
00138               catch ( ... ) {}
00139               CaloEvent *evt =  new CaloEvent();
00140               TBranch *branch= tree->GetBranch("CaloEvent");
00141               branch->SetAddress(&evt);
00142 
00143               CaloHit* hit =NULL;
00144               int nbEntries = tree->GetEntries();
00145               cout << " Nb of events to read " << nbEntries << endl;
00146               TString name = Form("deltaTRun_%i",RunId);
00147               TString title = Form("delta distribution Run %i",RunId);
00148               TH1I * deltaTRun = new TH1I(name,title,20,-10,10);
00149     
00150               for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
00151                 {
00152                   tree->GetEntry(evtNum);
00153                   nbEvt++;
00154                   for (int i=0; i<50; i++){hxyMap[i]->Reset();}
00155                   if(nbEvt>maxEvents) break;
00156 
00157                   int nbHits = evt->GetNHits();
00158                   hnhit->Fill(nbHits);
00159 
00160                   //number of hits in region of interest
00161                   float xpos = 0;
00162                   float ypos = 0;
00163                   float zpos = 0;
00164                   int ncol = 0;
00165                   int nrow = 0;
00166                   int zSlot = 0;
00167                   int deltat=0;
00168                   int nhit_roi = 0;
00169                   float Rad1lay=10000.;
00170                   int maxDeltat = evt->GetDeltaTmax();
00171                   int XOrg=1000; int YOrg=1000;
00172                   for (int layer=0 ; layer <100 ; layer++){NbHitLayer[layer]=0;NbHitLayer77[layer]=0;} 
00173 
00174                   for(int i=0 ; i < nbHits ; i++)
00175                     {
00176                       hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
00177                       deltat=hit->GetDeltaT();
00178                       deltaTRun->Fill(maxDeltat-deltat);
00179 
00180                       xpos =  hit->GetX()/10000. ; // unit = mm
00181                       ypos =  hit->GetY()/10000. ;
00182                       zSlot = hit->GetLayer() ;
00183                       zpos =  hit->GetZ()/10000.;
00184                       ncol = hit->GetColInChamber();
00185                       nrow = hit->GetRowInChamber();
00186 
00187                       if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs
00188                       if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas
00189 
00190                       nhit_roi++; // Nb of good hits
00191 
00192                       hz->Fill(zSlot);
00193 
00194                       hxy->Fill(ypos,xpos);
00195                       hxyMap[zSlot]->Fill(ncol,nrow,1);
00196                       hxyMapTot[zSlot]->Fill(ncol,nrow,1);
00197                       hxyTot->Fill(ncol,nrow);
00198 
00199                     } // end hit
00200                   hnhit_roi->Fill(nhit_roi);
00201 
00202                   if( nhit_roi<150 ) continue; // Reject muons already
00203 
00204                   int lG1=4 ; int lG2=15 ; 
00205                   int NbHitsG1 = 0;       int NbHitsG2 = 0;       int NbHitsG3 = 0;  
00206                   YOrg=hxyTot->ProjectionX()->GetMaximumBin()-1; // We get the central pad where beam goes for this event
00207                   XOrg=hxyTot->ProjectionY()->GetMaximumBin()-1;
00208                   //      cout << " YOrg " << YOrg << " XOrg " << XOrg << endl;
00209 
00210                   for (int kl=0 ; kl<=49 ; kl++){
00211                     NbHitLayer[kl] = NbHitLayer[kl] +  (int)hxyMap[kl]->GetEntries();
00212                     NbHitLayer77[kl] = hxyMap[kl]->Integral(YOrg-10,XOrg-10,YOrg+10, XOrg+10);
00213                     if(kl <=lG1 ) { NbHitsG1 = NbHitsG1 + (int)hxyMap[kl]->GetEntries();}
00214                     if(kl > lG1 && kl< lG2 ) { NbHitsG2 = NbHitsG2+ (int)hxyMap[kl]->GetEntries();}
00215                     if(kl > lG2 ) { NbHitsG3 = NbHitsG3 + (int)hxyMap[kl]->GetEntries();}
00216                   }
00217                   RG1G3 = (float)NbHitsG1 / NbHitsG3;
00218                 
00219                   // Compute Hit Density to find MIP events
00220                   int NbLayers=0;
00221                   for (int layer=0 ; layer < 50 ; layer++){
00222                     if(NbHitLayer[layer]>0)
00223                       {
00224                         NbLayers++;
00225                       }
00226                   }
00227 
00228                   float HitDensity = 1. * nhit_roi/NbLayers;
00229                   hDensity->Fill(HitDensity);
00230 
00231                   //Select Pion events
00232                   if( nhit_roi<150 ||  HitDensity<4  || RG1G3 > 1.0 ) continue;
00233 
00234                   // Find where the shower starts
00235 
00236                   NbLayers=0;
00237                   int ShowerStart=-10;
00238                   int ShowerStartDa=-10;
00239 
00240                   for (int layer=1 ; layer < 48 ; layer++){
00241                     if(NbHitLayer[layer]>0)
00242                       {
00243                         NbLayers++;
00244                         tp_nhit->Fill(layer,NbHitLayer[layer]); 
00245                         if(ShowerStartDa<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer]>=2 && NbHitLayer77[layer+1]>=4) {ShowerStartDa=layer; }
00246                         if(ShowerStart<0 && NbHitLayer77[layer-1]<=3 && NbHitLayer77[layer+1]-NbHitLayer77[layer]>2 && NbHitLayer77[layer+2]-NbHitLayer77[layer+1]>2) {
00247                           //if(ShowerStart<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer+1]>NbHitLayer77[layer] && NbHitLayer77[layer+2]>NbHitLayer77[layer+1]) {
00248                           //                      if(ShowerStart<0 && NbHitLayer77[layer-1]<=2 && NbHitLayer77[layer]>4) {
00249                           ShowerStart = layer; // cout << ShowerStart << endl ;
00250 
00251                         }
00252                       }
00253                   }
00254                   hShowerStart->Fill(ShowerStart);
00255                   hShowerStartDa->Fill(ShowerStartDa);
00256 
00257                   int uMgaslayer=45;
00258                   if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);}
00259                   if(ShowerStart<=uMgaslayer){hLongprof1->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);}
00260                   uMgaslayer=9;
00261                   if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);}
00262                   uMgaslayer=20;
00263                   if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);}
00264                   uMgaslayer=32;
00265                   if(ShowerStart<=uMgaslayer){hLongprof->Fill(uMgaslayer-ShowerStart,NbHitLayer[uMgaslayer]);}
00266 
00267                   if(ShowerStart==1){
00268                     for (int layer=0 ; layer < 48 ; layer++){
00269                       hLongprofL1->Fill(layer,NbHitLayer[layer]);}
00270                   }
00271          
00272                   if (nbEvt%10000 == 0) {
00273                     cout << "Processing event: " << nbEvt
00274                          << endl;  
00275                   }  
00276                   tf.cd(0);
00277                 }  // end event 
00278               cout << " Finished with events in this TREE " << endl;
00279               deltaTRun->Write(); 
00280               delete deltaTRun;
00281             } // else 
00282         } // end of tree or run
00283       cout << " Ending  file in chain " << endl;
00284     } // end of chain
00285   cout << " Going to finsh the work after" << nbEvt << " events "<< endl;
00286   tf.cd();
00287   hDensity->Write();
00288   hxy->Write(); 
00289   hz->Write();
00290   hnhit->Write();
00291   tp_nhit->Write();
00292   hnhit_roi->Write();
00293   hShowerStart->Write();
00294   hShowerStartDa->Write();
00295   hLongprof->Write();
00296   hLongprof1->Write();
00297   hLongprofL1->Write();
00298   for (int i=0 ; i<50 ; i++){
00299     hxyMapTot[i]->Write();
00300   }
00301 
00302 }
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/getShowerStart.cpp
