Micromegas LAPP reconstruction

00001 /* @version $Revision: 1376 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2011-11-18 16:11:02 +0100 (Fri, 18 Nov 2011) $ */
00002 #include <map>
00003 #include "TKey.h"
00004 #include <TH2F.h>
00005 #include <TProfile.h>
00006 #include <TROOT.h>
00007 #include <TRint.h>
00008 #include <TKey.h>
00009 #include <TList.h>
00010 #include <TTree.h>
00011 #include <TFile.h>
00012 #include <TChain.h>
00013 #include "TChainElement.h"
00014 #include <TMath.h>
00015 #include <stdlib.h>
00016 #include <TSystem.h>
00017 #include <iostream>
00018 #include <sstream>
00019 #include "Log.hh"
00020 
00021 
00022 #include "MicroException.hh"
00023 #include <string>
00024 
00025 #include "mysql/Mysql.hh"
00026 
00027 
00028 #include "root/CaloRun.hh"
00029 #include "root/CaloEvent.hh"
00030 #include "root/CaloHit.hh"
00031 #include "root/MTRun.hh"
00032 #include "root/MTEvent.hh"
00033 #include "root/MTChannel.hh"
00034 
00035 #include "TGraphErrors.h"
00036 #include "TGraphAsymmErrors.h"
00037 #include <TF1.h>
00038 
00039 using namespace std;
00040 const int maxEvents = 3000000; // How many events to read
00041 const int maxTrees = 100;  // How many trees to read
00042 const int maxHits = 200; // Max nb of considered hits
00043 
00044 
00045 int main(int argc, char**argv){
00046 
00047   if ( argc !=2  ) {
00048     FILE_LOG(logERROR)  << "usage: example rootFile " << endl;
00049     exit(1);
00050   }
00051 
00052   int search = 1;
00053 
00054   TH2F * hxy = new TH2F("hxy","",100,0,100,100,0,100);
00055   TH1I * hnhit = new TH1I("hnhit","Number of all hits",200,0,1000);
00056   TH1F * hDensity = new TH1F("hDensity","Hits Density ",200,0,50.);
00057   TH1I * hnhit_roi = new TH1I("hnhit_roi","Number of hits in tight DeltaT",130,0,130);
00058   TH1I * hnHitForw = new TH1I("hnhitForw","MIP Number of hits Forw",50,0,50);
00059   TProfile * tp_nhit = new TProfile("tp_nhit","",60,0,60,0,200,"");
00060   TH1F * hxdir = new TH1F("hxdir"," x Angle between hit and all the next ", 200, -3., 3.);
00061   TH1F * hydir = new TH1F("hydir"," y Angle between hit and all the next ", 200, -3., 3.);
00062   TProfile * hLayerMult = new TProfile("hLayerMult","Average Multiplicity around track per layer Y",50,0,50, 0.,1000,"");
00063   TH1F * hChi2X = new TH1F("hChi2X"," Chi2 X ", 20, 0., 20.);
00064   TH1F * hChi2Y = new TH1F("hChi2Y"," Chi2 Y ", 20, 0., 20.);
00065   TH1F * hResidY = new TH1F("hResidY"," Residual Y ", 60, -30., 30.);
00066   TH1F * hResidX = new TH1F("hResidX"," Residual X ", 60, -30., 30.);
00067   TProfile * pResidX = new TProfile("pResidX","Average Residual per layer X",50,0,50,-1000.,1000,"");
00068   TProfile * pResidY = new TProfile("pResidY","Average Residual per layer Y",50,0,50,-1000.,1000,"");
00069   TH1F * hPenteY = new TH1F("hPentY"," Slope Y ", 200, -100.,100.);
00070   TH1F * hPenteX = new TH1F("hPentX"," Slope X ", 200, -100, 100.);
00071   TH1F * hOffsetX = new TH1F("hOffsetX"," Offset X ", 200, 0, 1000.);
00072   TH1F * hOffsetY = new TH1F("hOffsetY"," Offset Y ", 200, 0, 1000.);
00073 
00074   TH2F * hxyLayer0 = new TH2F("hxyLayer0","X verus Y fitted track. Layer 0",100,0,1000,100,0,1000);
00075   TH1I * hz = new TH1I("hz","",100,0,100);
00076   TH1I * hGlNbHtsFound  = new TH1I("hGlNbHtsFound"," Nb of Hits found in Chamber", 50, 0, 50);
00077   TH1I * hGlNbHtsExtrap  = new TH1I("hGlNbHtsExtrap"," Nb of Hits extrapolated in Chamber", 50, 0, 50);
00078   TH1I * hNbHitsFit = new TH1I("hNbHitsFit","Nb of Hits Used in Fit", 100, 0, 200 );
00079 
00080   // Book Layer Hit maps and efficiencies
00081   TH2I * hxyMap[50];
00082   TH2I * hTrackxyMap[50];
00083   TH2F * hEffMap[50];
00084   TH2I * hRegxyMap[50] ;
00085   TH2I * hxyEvtMap[50] ;
00086 
00087   TString name,title;
00088   int NbinsEffx=12;
00089   int NbinsEffy=12;
00090   int hbin=96;
00091   for (int i=0;i<50;i++)
00092     {
00093       name = Form("hxyMap_%i",i);
00094       title = Form("Hit distribution in Layer %i",i);
00095       hxyMap[i]= new TH2I(name,title,96,0,hbin,96,0,hbin);
00096       name = Form("hTrackxyMap_%i",i);
00097       hTrackxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
00098       name = Form("hEffMap_%i",i);            
00099       hEffMap[i]= new TH2F(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
00100       name = Form("hRegxyMap_%i",i);
00101       hRegxyMap[i]= new TH2I(name,title,NbinsEffx,0,hbin,NbinsEffy,0,hbin);
00102     }
00103 
00104   string rootName;
00105   rootName.assign(argv[1]);
00106   int nbHit = 0;
00107   int nbEvt = 0;
00108 
00109   int NbHitLayer[100];
00110   int IsPenetratingMIP;
00111   //  float Hits[2][maxHits];
00112   Double_t Zlayer[maxHits]; // Before it was Float
00113   Double_t Xlayer[maxHits];
00114   Double_t Ylayer[maxHits];
00115   Double_t ErPos[maxHits];
00116   TString OutFilename = "Muons"+rootName+".root";
00117   TFile tf(OutFilename,"RECREATE");
00118   int NbTrees = 0;
00119   int RunId=1;
00120 
00121   TChain chain("T");
00122   for (int i =0; i<5 ; i++){
00123     //    if(i==2)continue;
00124     TString SubFile = Form("_I%i",i);
00125     TString fname = Form("/lapp_data/LC/data/TB2012/SPS_H2_may_2012/Production/CaloHits/CaloHit_DHCAL_"+rootName+SubFile+"_0.slcio.root");
00126     chain.Add(fname);
00127     cout << fname << endl;
00128   }
00129 
00130   TObjArray *fileElements=chain.GetListOfFiles();
00131   TIter next(fileElements);
00132   TChainElement *chEl=0;
00133   int NchainFiles=0;
00134   while (( chEl=(TChainElement*)next() )) 
00135     {
00136       TFile f(chEl->GetTitle());
00137       //  TFile f(rootName.c_str());
00138       TIter nextkey(f.GetListOfKeys());
00139       TKey *key;
00140       NchainFiles++;
00141       cout << " Reading a new file in the chain " << NchainFiles << endl;
00142       while (key = (TKey*)nextkey()) 
00143         {
00144 
00145           TString keyname = key->GetName();
00146           if ((keyname == "hResidualX") || (keyname == "hResidualY")){}
00147           else
00148             {
00149 
00150               TTree *tree = (TTree*)key->ReadObj();                
00151               cout << "---------------------------- NEW TTREE -----------------"<< endl;
00152               nbHit = 0;
00153               NbTrees++;
00154               if( NbTrees > maxTrees ) continue;
00155               CaloRun *run=NULL;
00156 
00157               try {
00158                 CaloRun *run=dynamic_cast<CaloRun *>(tree->GetUserInfo()->First());
00159                 if( run != NULL ){cout << "Run[" << run->GetRunId() <<  "], temperature[" << run->GetTemperature() <<  "] pression[" << run->GetPressure() << "] " << endl ; }
00160 
00161               }
00162               catch ( ... ) {}
00163               CaloEvent *evt =  new CaloEvent();
00164               TBranch *branch= tree->GetBranch("CaloEvent");
00165               branch->SetAddress(&evt);
00166 
00167               CaloHit* hit =NULL;
00168               int nbEntries = tree->GetEntries();
00169               cout << " Nb of events to read " << nbEntries << endl;
00170               TString name = Form("deltaTRun_%i",RunId);
00171               TString title = Form("delta distribution Run %i",RunId);
00172               TH1I * deltaTRun = new TH1I(name,title,20,-10,10);
00173     
00174               for ( int evtNum = 0; evtNum < nbEntries ; evtNum++)
00175                 {
00176                   tree->GetEntry(evtNum);
00177                   nbEvt++;
00178                   //              if(nbEvt%2000==0){cout << nbEvt << endl;}
00179                   if(nbEvt>maxEvents) continue;
00180                   /*
00181                     cout << "**** Event Information "<< dec << evtNum+1 << "/" << nbEntries<< " ****" <<  endl;
00182                     cout << "    event num in TTree:" << evtNum <<endl;
00183                     cout << "    event id:" << evt->fEventId <<endl;
00184 
00185                     cout << "    number of channel hit:" << evt->GetNHits() <<endl;
00186                   */
00187                   int nbHits = evt->GetNHits();
00188                   //      cout << "nb hits:" << nbHits << endl;
00189 
00190                   hnhit->Fill(nbHits);
00191 
00192                   if (nbHits > 120) continue;
00193 
00194                   //number of hits in region of interest
00195                   int nhit_roi = 0;
00196                   int nhit_fit = 0;
00197                   float xpos = 0;
00198                   float ypos = 0;
00199                   float zpos = 0;
00200                   unsigned int ncol = 0;
00201                   unsigned int nrow = 0;
00202                   int zSlot = 0;
00203                   int deltat=0;
00204                   int maxDeltat = evt->GetDeltaTmax();
00205 
00206                   for (int layer=0 ; layer <100 ; layer++){
00207                     NbHitLayer[layer]=0;}
00208                   for (int i=0 ; i<maxHits ;i++){
00209                     Zlayer[i]=0;
00210                     Xlayer[i]=0;
00211                     Ylayer[i]=0;
00212                     double theta = (0.0136/100.)*sqrt((i+1)*2/1.757)*(1+0.038*log((i+1)*2/1.757));
00213                     ErPos[i] = sqrt(pow(10./sqrt(12.),2)+pow(2*28.*theta,2));
00214                   }
00215 
00216                   for (int i=0; i<50; i++) // Fill Maps for this event
00217                     {
00218                       name = Form("hxyEvtMap_%i",i);
00219                       title = Form("Hit distribution in Layer %i",i);
00220                       hxyEvtMap[i]= new TH2I(name,title,96,0,96,96,0,96);
00221                     }
00222 
00223                   for(int i=0 ; i < nbHits ; i++)
00224                     {
00225                       hit = (CaloHit*)evt->GetHits()->UncheckedAt(i);
00226                       deltat=hit->GetDeltaT();
00227                       deltaTRun->Fill(maxDeltat-deltat);
00228 
00229                       xpos =  hit->GetX()/1000 ; // unit = mm
00230                       ypos =  hit->GetY()/1000 ;
00231                       zSlot = hit->GetLayer() ;
00232                       zpos =  hit->GetZ()/1000;
00233                       ncol = hit->GetColInChamber();
00234                       nrow = hit->GetRowInChamber();
00235 
00236                       if(zSlot <48 && abs(deltat-maxDeltat)>2) continue; // This is for RPCs
00237                       if(zSlot >=48 && abs(deltat-maxDeltat)>5) continue; // For Micromegas
00238 
00239                       //Eliminate Bad zones
00240                       //  if(zSlot==33 && xpos>64 ) continue;
00241                       // if(zSlot==49 && (xpos>=32 &&xpos<40) && ( ypos>=48 && ypos <56) ) continue;
00242 
00243                       hz->Fill(zSlot);
00244 
00245                       NbHitLayer[zSlot]++; 
00246 
00247                       hxy->Fill(ypos,xpos);
00248 
00249                       nhit_roi++; // Nb of good hits
00250                       hxyEvtMap[zSlot]->Fill(ncol,nrow,1);
00251                       hxyMap[zSlot]->Fill(ncol,nrow,1);
00252 
00253                       // Prepare for fitting a muon
00254                       // Choose layers to make the fit
00255                       if(nhit_fit==0){
00256                         Zlayer[nhit_fit] = zpos; // Convert coordinates to mm
00257                         Xlayer[nhit_fit] = xpos;
00258                         Ylayer[nhit_fit] = ypos;
00259                         nhit_fit++;
00260                       }
00261                       else{
00262                         if (nhit_fit < maxHits){
00263                           if(NbHitLayer[zSlot] >1 ) continue ; // keep only 1st hit in layers
00264                           // Select hits that have almost the same slope as the previous one            
00265                           Double_t slopy = TMath::ATan2( (ypos-Ylayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
00266                           Double_t slopx = TMath::ATan2( (xpos-Xlayer[nhit_fit-1]) , (zpos-Zlayer[nhit_fit-1]) );
00267                           hxdir->Fill(slopx); hydir->Fill(slopy);
00268                           //    hSlopyY->Fill(slopy);  hSlopx->Fill(slopx);
00269                           if ( TMath::Abs(slopx)> 1. || TMath::Abs(slopy)>1. ) continue;
00270                           Zlayer[nhit_fit] = zpos; 
00271                           Xlayer[nhit_fit] = xpos;
00272                           Ylayer[nhit_fit] = ypos;
00273                           nhit_fit++;
00274                         }
00275                       }
00276                     } // end hit
00277 
00278                   hnhit_roi->Fill(nhit_roi);
00279 
00280                   // Compute Hit Density to find MIP events
00281                   int NbLayers=0;
00282                   for (int layer=0 ; layer < 50 ; layer++){
00283                     if(NbHitLayer[layer]>0)
00284                       {
00285                         NbLayers++;
00286                         tp_nhit->Fill(layer,NbHitLayer[layer]);   
00287                       }
00288                   }
00289 
00290                   float HitDensity = 1. * nhit_roi/NbLayers;
00291 
00292                   hDensity->Fill(HitDensity);
00293                  
00294                   //Look for penetrating MIPS
00295 
00296                   int NbLayersForw=0;
00297                   int NbLayersBack=0;
00298 
00299                   IsPenetratingMIP = 0;
00300                   if(1 < HitDensity < 3){
00301 
00302                     for (int layer=0 ; layer < 10 ; layer++){
00303                       if(0 < NbHitLayer[layer] && NbHitLayer[layer] <= 3){NbLayersForw++;}
00304                     }
00305                     for (int layer=38 ; layer < 48 ; layer++){
00306                       if(0 < NbHitLayer[layer]  && NbHitLayer[layer] <= 3){NbLayersBack++;}
00307                     }
00308 
00309                     hnHitForw->Fill(NbLayersForw);
00310                     hnHitForw->Fill(NbLayersBack+20);
00311 
00312                     //Define penetrating MIP
00313                     if(NbLayersForw>6 && NbLayersBack > 6) { IsPenetratingMIP = 1;}
00314                   } // end Penetrating MIPS
00315 
00316 
00317                   // Fit a straight line for penetrating MIPs 
00318                   if(IsPenetratingMIP==1 &&  (nhit_fit > 25) ){
00319                     hNbHitsFit->Fill(nhit_fit); // How many hits are used for the fit     
00320                     float zmax = 1372.;
00321                     TGraphErrors *grx = new TGraphErrors(nhit_fit,Zlayer,Xlayer,0,ErPos);
00322                     TF1 *f1 = new TF1("f1","pol1",0,zmax);
00323                     grx->Fit("f1","Q");
00324                     double betaX = f1->GetParameter(0);
00325                     double alfaX = f1->GetParameter(1);
00326                     double chi2X = f1->GetChisquare();
00327                     delete f1; delete grx;
00328                     TGraphErrors *gry = new TGraphErrors(nhit_fit,Zlayer,Ylayer,0,ErPos);
00329                     TF1 *f2 = new TF1("f2","pol1",0,zmax);
00330                     gry->Fit("f2","Q");
00331                     double betaY = f2->GetParameter(0);
00332                     double alfaY = f2->GetParameter(1);
00333                     double chi2Y = f2->GetChisquare();
00334                     // if(nbEvt<100) gry->Write();
00335                     delete f2; delete gry;
00336                     hPenteX->Fill(alfaX*1000.);
00337                     hPenteY->Fill(alfaY*1000.);
00338                     hOffsetX->Fill(betaX);
00339                     hOffsetY->Fill(betaY);
00340                     hChi2X->Fill(chi2X/(nhit_fit-2));
00341                     hChi2Y->Fill(chi2Y/(nhit_fit-2));
00342                     hxyLayer0->Fill(betaX,betaY); // Make X Y Map at the first layer
00343 
00344                     // Compute residuals to fitted track
00345                     for (int j=1; j<nhit_fit ; j++){
00346                       float ResX = Xlayer[j]-( alfaX*Zlayer[j]+betaX );
00347                       float ResY = Ylayer[j]-( alfaY*Zlayer[j]+betaY );
00348                       hResidY->Fill(ResY);
00349                       hResidX->Fill(ResX);
00350                       pResidX->Fill(Zlayer[j]/28,ResX);
00351                       pResidY->Fill(Zlayer[j]/28,ResY);
00352                     }
00353 
00354                     // Go to compute Efficiency with good tracks
00355                     if( chi2X/(nhit_fit-2) < 6. && chi2Y/(nhit_fit-2) < 6 ) {
00356                     for(int i=0; i<50 ; i++){
00357                       int NpadsLook = 3;
00358                       float Xextrapol = alfaX*i*28. + betaX ; // X extrapolate to chamber i in mm
00359                       float Yextrapol = alfaY*i*28. + betaY ; // Y extrapolate to chamber i in mm
00360 
00361                       // Find Central pad from computed x,y position
00362                       Int_t xpadCentr = -1000;
00363                       Int_t ypadCentr = -1000;
00364                       if( i>47 && ( fabs(Yextrapol-322.5) < 20 || fabs(Yextrapol-650) < 20 ) ) continue;
00365                       if( i>47 && ( fabs(Xextrapol-485.0) < 20  ) ) continue;
00366                       if(i<48){
00367                         xpadCentr = (int)(Xextrapol/10.4);
00368                         ypadCentr = (int)(Yextrapol/10.4);
00369                       }
00370                       else{
00371                         if(Xextrapol<=480) {xpadCentr = (int)(Xextrapol/10.);}
00372                         else {xpadCentr = (int)( (Xextrapol-5)/10 );}
00373                         if(Yextrapol<=320.) {ypadCentr = (int)( (Yextrapol/10.) );}
00374                         if(Yextrapol>325 && Yextrapol<=645.) {ypadCentr = (int)( (Yextrapol-5.)/10.);}
00375                         if(Yextrapol>=646) {ypadCentr = (int)((Yextrapol-10.)/10.);}
00376                       }
00377                       // check that we are inside the chamber 
00378                       if(xpadCentr < 0  || xpadCentr > 95 ) continue;
00379                       if(ypadCentr < 0  || ypadCentr > 95 ) continue;
00380                       // define fiducial volume
00381                       if(xpadCentr <= NpadsLook || xpadCentr>95-NpadsLook) continue ;
00382                       if(ypadCentr <= NpadsLook || ypadCentr>95-NpadsLook) continue ;
00383 
00384                       hTrackxyMap[i]->Fill(ypadCentr,xpadCentr,1);  //Fill extrapolated point
00385                       hGlNbHtsExtrap->Fill(i);
00386                       //Check now if there is a hit around this point in a region of NpadsLook
00387                       int xbinL = max(xpadCentr - NpadsLook,0); // For the PS data we must increase the search zone in x
00388                       int xbinH = min(xpadCentr + NpadsLook,95);
00389                       int ybinL = max(ypadCentr - NpadsLook,0);
00390                       int ybinH = min(ypadCentr + NpadsLook,95);
00391                       int NhitsRegSearch = hxyEvtMap[i]->Integral(ybinL,ybinH,xbinL,xbinH);
00392                       if( NhitsRegSearch >0 ) {
00393                         hRegxyMap[i]->Fill(ypadCentr,xpadCentr,1);  //Fill extrapolated point
00394                         hGlNbHtsFound->Fill(i) ;
00395                         hLayerMult->Fill(i, NhitsRegSearch) ; }
00396                     
00397                     } // end of efficiency coputing end layer loop
00398                     } // end of Chi2 selection
00399                   } // end of penetrating MIP
00400  
00401                   if (nbEvt%10000 == 0) {
00402                     cout << "Processing event: " << nbEvt
00403                          << endl;  
00404                   }  
00405                   for (int i=0; i<50 ; i++){
00406                     delete hxyEvtMap[i];
00407                   } 
00408                   tf.cd(0);
00409                 }  // end event 
00410               cout << " Finished with events in this TREE " << endl;
00411               deltaTRun->Write(); 
00412               delete deltaTRun;
00413             } // else 
00414         } // end of tree or run
00415       cout << " Ending  file in chain " << endl;
00416     } // end of chain
00417   cout << " Going to finsh the work after" << nbEvt << " events "<< endl;
00418   tf.cd();
00419   hPenteY->Write();
00420   hPenteX->Write();
00421   hOffsetX->Write();
00422   hOffsetY->Write();
00423   hxdir->Write();
00424   hydir->Write();
00425   tp_nhit->Write();
00426   hLayerMult->Write();
00427   hxy->Write();
00428   hxyMap[33]->Write();
00429   hxyMap[31]->Write();
00430   hxyMap[48]->Write();
00431   hxyMap[49]->Write();
00432   hz->Write();
00433   hnhit->Write();
00434   hnhit_roi->Write();
00435   hDensity->Write();
00436   hnHitForw->Write();
00437   hNbHitsFit->Write();
00438   hChi2X->Write();
00439   hChi2Y->Write();
00440   hResidY->Write();
00441   hResidX->Write();
00442   pResidX->Write();
00443   pResidY->Write();
00444   hxyLayer0->Write();
00445   for (int i=0 ; i<50 ; i++){
00446     hRegxyMap[i]->Write();
00447     hTrackxyMap[i]->Write();
00448     hEffMap[i]->Divide(hRegxyMap[i], hTrackxyMap[i], 100., 1.,"b"); // Binomial errors
00449     hEffMap[i]->Write();
00450   }
00451   TH1F * hEffRPC = new TH1F("hEffRPC"," RPC Efficiencies for every chip ", 100, 0., 100.);
00452   TH1F * hEffmMegas = new TH1F("hEffmMegas"," mMegas Efficiencies for every chip ", 100, 0., 100.);
00453   for (int i = 0; i< 50; i++){
00454     for (int j = 1; j<13 ; j++){
00455       for (int k=1; k<13;  k++) {
00456         float eff = hEffMap[i]->GetBinContent(j,k);
00457         if(i<48) {  hEffRPC->Fill(eff);   }
00458         if(i>=48){  hEffmMegas->Fill(eff); }
00459       }// end columns
00460     } //end rows
00461   }// end Chamber
00462 
00463   hEffRPC->Write();
00464   hEffmMegas->Write();
00465 
00466   // Compute Global efficiency
00467   TH1F * hGlEffChamber = new TH1F("hGlEffChamber"," Chamber Global Efficiency ", 50, 0, 50 );
00468   hGlEffChamber->Divide(hGlNbHtsFound, hGlNbHtsExtrap,  100. , 1., "b"); // Binomial errors
00469   hGlNbHtsFound->Write();
00470   hGlNbHtsExtrap->Write();
00471   hGlEffChamber->Write();
00472   for (int j=1; j<=50 ; j++){
00473     double nbTrFound= hGlNbHtsFound->GetBinContent(j);
00474     double nbTrExtr=hGlNbHtsExtrap->GetBinContent(j);
00475     double ef=hGlEffChamber->GetBinContent(j);
00476     double erEff=hGlEffChamber->GetBinError(j);
00477     double w = nbTrFound/nbTrExtr;
00478     double binErr=w*(1-w)/nbTrExtr;
00479     cout << "chamber " << j << " Nb tracks found and extrapol " << nbTrFound << "  " << nbTrExtr << " Efficiency " << ef << " +- " << binErr << endl;    
00480   }
00481   TGraphAsymmErrors * gr = new TGraphAsymmErrors(hGlNbHtsFound, hGlNbHtsExtrap,"");
00482   gr->Write();   
00483 }
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/getChamberEfficiency.cpp
