Micromegas LAPP reconstruction

00001 /* @version $Revision: 1566 $ * @modifiedby $Author: jacquem $ * @lastmodified $Date: 2012-03-08 14:00:10 +0100 (Thu, 08 Mar 2012) $ */
00002 /* 
00003  * File:   beamProfile.cpp
00004  * Author: jacquem
00005  *
00006  * Created on May 6, 2010, 10:03 AM
00007  */
00008 
00009 // System includes
00010 #include <iostream>
00011 #include <stdlib.h>
00012 #include <map>
00013 #include <string>
00014 
00015 // Root includes
00016 #include <TFile.h>
00017 #include <TTree.h>
00018 #include <TKey.h>
00019 #include <TH1I.h>
00020 #include <TH2I.h>
00021 #include <TText.h>
00022 #include <TCanvas.h>
00023 #include <TPaveLabel.h>
00024 
00025 // Framework inlcudes
00026 #include "Run.hh"
00027 #include "Detector.hh"
00028 #include "XMLTool.hh"
00029 #include "Toolbox.hh"
00030 #include "MicroException.hh"
00031 #include "Log.hh"
00032 
00033 #include "root/MTRun.hh"
00034 #include "root/MTEvent.hh"
00035 #include "root/MTChannel.hh"
00036 #include "root/MTChamber.hh"
00037 #include "root/MTDetector.hh"
00038 #include "tools/Arguments.hh"
00039 
00040 using namespace std;
00041 
00042 
00043 
00044 int main(int argc, char** argv) {
00045 
00046 
00047   FILELog::ReportingLevel() = FILELog::FromString(INFO);
00048 
00049   Arguments arg(argc,argv,"monitoring inputFile outputFile -l evtNum");
00050   if (arg.getNbOfArgs() != 2 ) { arg.usage() ; exit(-1);}
00051   string rootName = arg.getArgument(1);
00052   string outputFileName = arg.getArgument(2);
00053 
00054    
00055   vector<TH1I* > dtVec;
00056 
00057  map<unsigned int, TH1F* > tempAsu1Map; // key is DifId
00058  map<unsigned int, TH1F* > tempAsu2Map; // key is DifId
00059  map<unsigned int, TH1F* > tempDifMap; // key is DifId
00060 
00061 
00062   UInt_t eventLimit = 0;
00063   if ( arg.getOption("-l").size() != 0 )
00064   {
00065     eventLimit = atoi(arg.getOption("-l").c_str());
00066   }
00067 
00068   // XY histos
00069   std::map<ui16,TH2I*> xyKeyChamber; // key is chamberId
00070   bool construct = false;
00071 
00072   TFile f(rootName.c_str());
00073 
00074   UInt_t nBoarderEvt = 0;
00075   UInt_t nCRCError = 0;
00076   UInt_t nBadChipId = 0;
00077   UInt_t nUnknowError = 0;
00078   UInt_t nNotValidEvent = 0;
00079  
00080 
00081   TH2I* xy = NULL;
00082   TIter nextkey(f.GetListOfKeys());
00083   TKey *key;
00084   while (key = (TKey*)nextkey()) 
00085   {
00086     TTree *tree = (TTree*)key->ReadObj();
00087     FILE_LOG(logINFO  ) << "New TTree [" << tree->GetName() << endl;
00088     MTRun *run=NULL;
00089     try
00090     {
00091       run=dynamic_cast<MTRun *>(tree->GetUserInfo()->First());
00092       if (construct == false)
00093       {
00094         MTDetector* det = run->GetDetector();
00095         if ( det == NULL ) { exit(-1);}
00096         // loop over chamber config
00097         const std::map<UInt_t,MTChamber*>& chambers = det->GetChambers();
00098         for( std::map<UInt_t,MTChamber*>::const_iterator ii=chambers.begin(); ii!=chambers.end(); ++ii)
00099         {
00100           const MTChamber& chamber= *((*ii).second);
00101           
00102           string xyname("hxy_");
00103           Toolbox::addIntToString(xyname,chamber.GetId());
00104 
00105           string  xytitle("XY distribution of hits in chamber ");
00106           Toolbox::addIntToString(xytitle,chamber.GetId());
00107 
00108           
00109           // create histogram
00110           if ( chamber.GetId() >=48 ) // MICROMEGAS
00111           {                                                   // 1 m -> 1 000 000 um
00112             xy = new TH2I(xyname.c_str(),xytitle.c_str(),193,0,965000,194,0,970000);
00113           }
00114           else // RPC
00115           {
00116             Int_t max = 10400*96; // 10400 micro * 96 pad
00117             xy = new TH2I(xyname.c_str(),xytitle.c_str(),96,0,max,96,0,max);
00118           }
00119 
00120           xy->SetXTitle("y");
00121           xy->SetYTitle("x");
00122           //cout << "Create XY dist for " << xytitle << endl;
00123 
00124           // store it to the map
00125           xyKeyChamber.insert(make_pair(chamber.GetId(),xy));
00126         }
00127       construct = true;
00128       }
00129     }
00130     catch (MicroException e) {}
00131 
00132     MTEvent *evt =  new MTEvent();
00133     TBranch *branch= tree->GetBranch("MTEvent");
00134     branch->SetAddress(&evt);
00135     MTChannel* channel =NULL;
00136     int nbEntries = tree->GetEntries();
00137     if ( eventLimit == 0 ) { eventLimit = nbEntries; }
00138     unsigned int interval = (unsigned int) eventLimit / 10;
00139     if (interval == 0) interval = 1;
00140     for ( int evtNum = 0; evtNum <eventLimit ; evtNum++)
00141     {
00142       tree->GetEntry(evtNum);
00143  
00144 /*
00145     const map<unsigned int , float>& tempDif = evt->GetTemperatureDif();
00146     for ( map<unsigned int , float>::const_iterator it = tempDif.begin() ; it != tempDif.end(); it++)
00147     {
00148       if ( != 0) { asu1TempHist->Fill(evt->GetTemperatureAsu1()); } 
00149       cout << "    temperature dif: " << it->first << " dif:" << it->second <<endl;
00150     }
00151       if ( evt->GetTemperatureAsu2() != 0) { asu2TempHist->Fill(evt->GetTemperatureAsu2()); } 
00152       if ( evt->GetTemperatureDif() != 0) { difTempHist->Fill(evt->GetTemperatureDif()); } 
00153 */
00154       
00155       UInt_t nbPads=100;     // input   2/3 du nombre de pad format le cadre d'un asu 1 ASU -> ( 6 chips * 8 pads * 2 cotes + 4 chips * 8 pad * 2 cote ) / 3 * 2
00156       UInt_t asu;            // output
00157       UInt_t abs_time;       // output
00158       UInt_t chamber;       // output
00159 
00160       if ( evt->IsSquare(nbPads,asu,abs_time, chamber) == true ) { nBoarderEvt++; } 
00161       if ( evt->GetCrcIsCorrect() == false ) {nCRCError++; } 
00162       if ( (evt->GetValidFlag() && UNKNOW_NOT_VALID) ==  UNKNOW_NOT_VALID ) { nUnknowError++; }
00163       if ( (evt->GetValidFlag() && BAD_CHIP_ID) ==  BAD_CHIP_ID ) { nBadChipId++; }
00164       if ( ! evt->IsValid() ) { nNotValidEvent++; }
00165 
00166 
00167 
00168 
00169       map<UInt_t, UInt_t> dtMap; // key is BcId_Dif() - BcId_Hit. Value is number of entry
00170       UInt_t maxDt = 0;
00171       UInt_t minDt = 0xFFFFFFFF;
00172 
00173       int nbChannel = evt->GetNchannel();
00174       for(int i=0;i<nbChannel  ;i++)
00175       {
00176         channel = (MTChannel*)evt->GetChannels()->UncheckedAt(i);
00177         UInt_t t2 = channel->GetBcId_Dif();
00178         UInt_t t3 = channel->GetBcId_Hit();
00179         UInt_t dt = t2 - t3;
00180         if ( dt > maxDt ) { maxDt = dt ; }
00181         if ( dt < minDt ) { minDt = dt ; }
00182 
00183 
00184         if (dtMap.find(dt) == dtMap.end())
00185         { //if this dt doesn't exist yet
00186           dtMap.insert(make_pair(dt,1));
00187         }
00188         else
00189         { //if this abs time already exist
00190           dtMap[dt] = dtMap[dt]++;
00191         }
00192 
00193   
00194 
00195         TH2I *hist = xyKeyChamber[channel->GetChamberId()];
00196         if ( hist ) 
00197         {
00198           if ( channel->GetChamberId() >=48 ) // MICROMEGAS
00199           {
00200               hist->Fill(channel->GetY() + 2500,channel->GetX()+2500,channel->GetDigitalValue());
00201               hist->Fill(channel->GetY()+2500,channel->GetX()-2500,channel->GetDigitalValue());
00202               hist->Fill(channel->GetY()-2500,channel->GetX()+2500,channel->GetDigitalValue());
00203               hist->Fill(channel->GetY()-2500,channel->GetX()-2500,channel->GetDigitalValue());
00204           }
00205           else  //RPC
00206           { 
00207               hist->Fill(channel->GetY(),channel->GetX(),channel->GetDigitalValue());
00208           }
00209 
00210         }
00211       }
00212       
00213       if ( evtNum % interval == 0 )
00214       {
00215         //cout << "Min Dt[" << minDt <<  "], maxDt[" << maxDt << "]" << flush;
00216         string name = "dt_" ;
00217         name = name + tree->GetName();
00218         Toolbox::addIntToString(name,evtNum);
00219         string title = "BcId_Dif - BcId_Hit";
00220         TH1I *dtHist = new TH1I(name.c_str(),title.c_str(), maxDt-minDt,  0 ,maxDt);
00221         
00222         name = "dt_projection";
00223         name = name + tree->GetName();
00224         Toolbox::addIntToString(name,evtNum);
00225 
00226         TH1I *dtProjHist = new TH1I(name.c_str(),title.c_str(), nbChannel,  0 ,nbChannel);
00227         for ( map<UInt_t ,UInt_t>::const_iterator iterC = dtMap.begin(); iterC != dtMap.end(); iterC++)
00228         {
00229           Int_t val = (*iterC).second;
00230           dtHist->Fill((*iterC).first,val);
00231           dtProjHist->Fill(val);
00232         }
00233         dtVec.push_back(dtHist);
00234         dtVec.push_back(dtProjHist);
00235       }
00236 
00237       // Asu1, 2 and Dif temperature
00238       const map<unsigned int , float>& tempAsu1 = evt->GetTemperatureAsu1();
00239       for ( map<unsigned int , float>::const_iterator it = tempAsu1.begin() ; it != tempAsu1.end(); it++)
00240       {
00241         unsigned int difId = it->first ;
00242         float temp = it->second;
00243         TH1F* foo = NULL;
00244         if ( tempAsu1Map.find(difId) == tempAsu1Map.end())
00245         {
00246           string name; 
00247           std::ostringstream oss;
00248           oss << "Dif_" << difId;
00249           name = name + oss.str();
00250           name = name + "_Asu1_temp" ;
00251           
00252           foo = new TH1F(name.c_str(), "Asu1 temperature", 100,  0 ,100);
00253           tempAsu1Map.insert(make_pair(difId,foo));
00254         }
00255         else
00256         {
00257           foo = tempAsu1Map.find(difId)->second;
00258         }
00259         foo->Fill(temp);
00260       }
00261 
00262       const map<unsigned int , float>& tempDif = evt->GetTemperatureDif();
00263       for ( map<unsigned int , float>::const_iterator it = tempDif.begin() ; it != tempDif.end(); it++)
00264       {       
00265         unsigned int difId = it->first ;
00266         float temp = it->second;
00267         TH1F* foo = NULL;
00268         if ( tempDifMap.find(difId) == tempDifMap.end())
00269         { 
00270           string name;
00271           std::ostringstream oss;
00272           oss << "Dif_" << difId;
00273           name = name + oss.str();
00274           name = name + "_Dif_temp" ;
00275        
00276           foo = new TH1F(name.c_str(), "Dif temperature", 100,  0 ,100);
00277           tempDifMap.insert(make_pair(difId,foo));
00278         }
00279         else
00280         {
00281           foo = tempDifMap.find(difId)->second;
00282         }
00283         foo->Fill(temp);
00284       } 
00285 
00286       const map<unsigned int , float>& tempAsu2 = evt->GetTemperatureAsu2();
00287       for ( map<unsigned int , float>::const_iterator it = tempAsu2.begin() ; it != tempAsu2.end(); it++)
00288       {       
00289         unsigned int difId = it->first ;
00290         float temp = it->second;
00291         TH1F* foo = NULL;
00292         if ( tempAsu2Map.find(difId) == tempAsu2Map.end())
00293         { 
00294           string name;
00295           std::ostringstream oss;
00296           oss << "Dif_" << difId;
00297           name = name + oss.str();
00298           name = name + "_Asu2_temp" ;
00299        
00300           foo = new TH1F(name.c_str(), "Asu2 temperature", 100,  0 ,100);
00301           tempAsu2Map.insert(make_pair(difId,foo));
00302         }
00303         else
00304         {
00305           foo = tempAsu2Map.find(difId)->second;
00306         }
00307         foo->Fill(temp);
00308       } 
00309 
00310       
00311       FILE_LOG(logINFO  ) << "Write Event[" << evtNum << "] \t\t\t\t\t\t\t\t\t\t\t\r" << flush;
00312     }
00313  }
00314 
00315 
00316   TFile output(outputFileName.c_str(),"RECREATE");
00317 
00318 
00319 
00320 
00321    TCanvas *statistic = new TCanvas("statistic", "statistic",100,10,700,500);
00322    statistic->Range(0,0,20,24);
00323    statistic->SetFillColor(10);
00324    statistic->SetBorderSize(2);
00325 
00326    TPaveLabel *pl = new TPaveLabel(3,22,17,23.7, "Monitoring result","br");
00327    pl->SetFillColor(18);
00328    pl->Draw();
00329 
00330    TText t(0,0,"a");
00331    t.SetTextFont(62);
00332    t.SetTextSize(0.025);
00333    t.SetTextAlign(12);
00334 
00335    string msg = "Total of not valid Event: ";
00336    Toolbox::addIntToString(msg,   nNotValidEvent); msg = msg + ( " / "); Toolbox::addIntToString(msg,   eventLimit );
00337    float p = float(float(nNotValidEvent) / float( eventLimit) * 100 );
00338    msg = msg + ( "-> "); Toolbox::addFloatToString(msg,  p ); msg = msg + ( "% ");
00339    t.DrawText(2,20,msg.c_str());
00340 
00341    msg = "Number of Event with CRC incorrect: ";
00342    Toolbox::addIntToString(msg,   nCRCError); msg = msg + ( " / "); Toolbox::addIntToString(msg,   eventLimit );
00343    p = float(float(nCRCError) / float( eventLimit) * 100 );
00344    msg = msg + ( "-> "); Toolbox::addFloatToString(msg, p ); msg = msg + ( "% ");
00345    t.DrawText(2,18,msg.c_str());
00346    
00347    msg = "Number of not valid Event because of  bad chip id: ";
00348    Toolbox::addIntToString(msg,   nBadChipId); msg = msg + ( " / "); Toolbox::addIntToString(msg,   eventLimit );
00349    p = float(float(nBadChipId) / float( eventLimit) * 100 );
00350    msg = msg + ( "-> "); Toolbox::addFloatToString(msg,  p ); msg = msg + ( "% ");
00351    t.DrawText(2,16,msg.c_str());
00352 
00353    msg = "Number of not valid Event because of unknow error: ";
00354    Toolbox::addIntToString(msg,   nUnknowError); msg = msg + ( " / "); Toolbox::addIntToString(msg,   eventLimit );
00355    p = float(float(nUnknowError) / float( eventLimit) * 100 );
00356    msg = msg + ( "-> "); Toolbox::addFloatToString(msg,  p ); msg = msg + ( "% ");
00357    t.DrawText(2,14,msg.c_str());
00358 
00359    msg = "Number of boarder Event (with  a least 100 boarder pad): ";
00360    Toolbox::addIntToString(msg,   nBoarderEvt); msg = msg + ( " / "); Toolbox::addIntToString(msg,   eventLimit );
00361    p = float(float(nBoarderEvt) / float( eventLimit) * 100 );
00362    msg = msg + ( "-> "); Toolbox::addFloatToString(msg,  p ); msg = msg + ( "% ");
00363    t.DrawText(2,10,msg.c_str());
00364 
00365    statistic->cd();
00366    statistic->Write();
00367 
00368   for ( map<ui16,TH2I*>::const_iterator iterC = xyKeyChamber.begin(); iterC != xyKeyChamber.end(); iterC++)
00369   {
00370     TH2I *hist = (*iterC).second;
00371     if ( hist )  
00372     {  
00373       hist->Write();
00374     } 
00375   }
00376 
00377  
00378   TCanvas deltaT("deltatT");
00379   double racine = sqrt(dtVec.size());
00380   UInt_t intSqrt = (UInt_t)(racine+1);
00381 
00382   deltaT.Divide(intSqrt,intSqrt);
00383   
00384   UInt_t index = 1;
00385 
00386   for (vector<TH1I* >::const_iterator iterC = dtVec.begin(); iterC != dtVec.end(); iterC++)
00387   {
00388     TH1I *hist = *iterC;
00389     if ( hist )
00390     { 
00391       if ( index % 2 == 0) deltaT.GetPad(index)->SetLogx();
00392       deltaT.GetPad(index)->SetLogy();
00393       deltaT.cd(index);
00394       index++;
00395       hist->Draw();
00396     }
00397   }
00398   deltaT.Write();
00399 
00400 
00401   for ( map<unsigned int,TH1F*>::const_iterator iterC = tempAsu1Map.begin(); iterC != tempAsu1Map.end(); iterC++)
00402   {
00403     UInt_t difId = (*iterC).first;
00404     string title("temperature_dif_");
00405     Toolbox::addIntToString(title,difId);
00406     TCanvas temperature(title.c_str());
00407     temperature.Divide(3);
00408 
00409     TH1F *asu1 = (*iterC).second;
00410     TH1F *asu2 = tempAsu2Map.find(difId)->second;
00411     TH1F *dif = tempDifMap.find(difId)->second;
00412     
00413     if ( asu1 )  { temperature.cd(1); asu1->Draw(); } 
00414     if ( asu2 )  {  temperature.cd(2); asu2->Draw(); } 
00415     if ( dif )  {  temperature.cd(3); dif->Draw(); } 
00416     temperature.Write();
00417   }
00418 
00419 
00420   output.Close();
00421   f.Close(); 
00422    
00423 
00424   cout <<endl << endl << "Done" << endl;
00425   
00426 return 0;
00427 }
00428 
00429 
00430   
00431
/data3/calcul/jacquem/working_dir/Micromegas/micromegasFrameWork/src/analyse/root/monitoring.cpp
