SiStripCommissioningAnalysis/test/CalibrationScanAnalysis.C

0001 #include <TFile.h>
0002 #include <TKey.h>
0003 #include <TClass.h>
0004 #include <TH1F.h>
0005 #include <TH2F.h>
0006 #include <TGraph.h>
0007 #include <TObjString.h>
0008 #include <TCanvas.h>
0009 #include <string>
0010 #include <map>
0011 #include <utility>
0012 #include <vector>
0013 #include <list>
0014 #include <algorithm>
0015 #include <iostream>
0016 #include <fstream>
0017 #include <sstream>
0018 #include <iomanip>
0019 #include <cmath>
0020 #include <cstring>
0021
0022 typedef std::pair<int, int> Parameters;
0023 typedef std::map<Parameters,TFile*> FileList;
0024 typedef std::map<Parameters,std::vector<TH1*> > SummaryV;
0025
0026 #define DATAPATH "/DQMData/Collate/SiStrip/"
0027 #define HISTOPATH "/DQMData/Collate/SiStrip/ControlView/"
0028
0029 //#define DEBUG_ON
0030
0031 class CalibrationScanAnalysis
0032 {
0033
0034   public:
0035     CalibrationScanAnalysis(bool tuneISHA = true, bool tuneVFS = true);
0036     virtual ~CalibrationScanAnalysis();
0037     void tuneISHA(bool tune) { tuneISHA_ = tune; }
0038     void tuneVFS(bool tune)  { tuneVFS_  = tune; }
0039     void addFile(const std::string&);
0040     void analyze();
0041     void sanitizeResult(unsigned int cut = 2, bool doItForISHA = true, bool doItForVFS = true);
0042     void print(Option_t* option = "") const;
0043     void draw(Option_t* option = "") const;
0044     void save(const char* fileName="-");
0045
0046   protected:
0047     void addFile(TFile* );
0048     void getSummaries(FileList::const_iterator);
0049     void sortByGeometry();
0050     void loadPresentValues();
0051     float getX(const TGraph*, const float&) const;
0052     bool checkInput() const;
0053     TH1F* fixHisto(std::vector<std::string>&,TH1*) const;
0054
0055   private:
0056     bool tuneISHA_, tuneVFS_;
0057     FileList files_;
0058     SummaryV summaries_;
0059     std::map<std::string, Parameters> result_;
0060     std::map<std::string, int> geometries_;
0061     std::map<std::string, Parameters> presentValues_;
0062 };
0063
0064 CalibrationScanAnalysis::CalibrationScanAnalysis(bool tuneISHA, bool tuneVFS):
0065     tuneISHA_(tuneISHA),tuneVFS_(tuneVFS) {
0066 }
0067
0068 CalibrationScanAnalysis::~CalibrationScanAnalysis() {
0069   // close and delete all files
0070   for(FileList::iterator file = files_.begin();file!=files_.end();++file) {
0071     // this will automatically delete histograms in summaries_
0072     file->second->Close();
0073     delete file->second;
0074   }
0075 }
0076
0077 void CalibrationScanAnalysis::addFile(const std::string& filename) {
0078   TFile* test = new TFile(filename.c_str());
0079   bool noFile = test->IsZombie();
0080   test->Close();
0081   delete test;
0082   if(!noFile) {
0083     TFile* newFile = new TFile(filename.c_str(),"UPDATE");
0084     addFile(newFile);
0085   }
0086 }
0087
0088 void CalibrationScanAnalysis::addFile(TFile* newFile) {
0089   int isha,vfs;
0090   TList* keyList = newFile->GetDirectory(DATAPATH)->GetListOfKeys();
0091   TIter next(keyList);
0092   TNamed* ishaObj = NULL;
0093   TNamed* vfsObj  = NULL;
0094   TNamed* obj = NULL;
0095   while ((obj = (TNamed*)(next()))) {
0096     if(strncmp(obj->GetName(),"<isha>",6)==0) ishaObj = (TNamed*)obj;
0097     if(strncmp(obj->GetName(),"<vfs>",5)==0)  vfsObj  = (TNamed*)obj;
0098   }
0099   if(!ishaObj || !vfsObj) {
0100      std::cerr << "Error: Unexpected file structure. ISHA/VFS values not found." << std::endl;
0101      newFile->Close();
0102      delete newFile;
0103      return;
0104   }
0105   isha = atoi(ishaObj->GetName()+8);
0106   vfs  = atoi(vfsObj->GetName()+7 );
0107   std::cout << "Loaded File for ISHA/VFS = " << isha << "/" << vfs << std::endl;
0108   files_[std::make_pair(isha,vfs)] = newFile;
0109 }
0110
0111 void CalibrationScanAnalysis::getSummaries(FileList::const_iterator file) {
0112   std::cout << "." << std::flush;
0113   std::vector<TH1*> result;
0114   TFile* input = file->second;
0115   TDirectory* directory = input->GetDirectory(HISTOPATH);
0116   TList* histograms = directory->GetListOfKeys();
0117   TIter next(histograms);
0118   TKey* key = NULL;
0119   while ((key = (TKey*)next())) {
0120     if(TClass(key->GetClassName()).InheritsFrom("TH1")) {
0121       TH1* h = (TH1*)key->ReadObj();
0122       result.push_back(h);
0123     }
0124   }
0125   summaries_[file->first] = result;
0126 }
0127
0128 float CalibrationScanAnalysis::getX(const TGraph* graph, const float& y) const {
0129    Double_t* arrayX = graph->GetX();
0130    Double_t* arrayY = graph->GetY();
0131    //first, look for an intersection
0132    for(int i=0;i<graph->GetN()-1;++i) {
0133      if((arrayY[i]-y)*(arrayY[i+1]-y)<0) {
0134        return (arrayX[i]+((arrayX[i+1]-arrayX[i])/(arrayY[i+1]-arrayY[i])*(y-arrayY[i])));
0135      }
0136    }
0137    // if none, look for a plateau
0138    float finalDelta = fabs(arrayY[graph->GetN()-1]-y);
0139    // allow for a 50% increase of the difference
0140    float delta = finalDelta*0.5;
0141    int lastpoint = graph->GetN()-1;
0142    for(int i=lastpoint-1;i>=0;--i) {
0143      if(fabs(arrayY[lastpoint]-arrayY[i])>delta)
0144        return arrayX[i+1];
0145    }
0146    // in last ressort, return the central value.
0147    return arrayX[lastpoint]-arrayX[0];
0148 }
0149
0150 void CalibrationScanAnalysis::analyze() {
0151
0152 #ifdef DEBUG_ON
0153   TFile* debugFile = new TFile("debug.root","RECREATE");
0154 #endif
0155
0156   // load data from files
0157   std::cout << "Loading data from files..." << std::endl;
0158   for(FileList::const_iterator it=files_.begin();it!=files_.end();++it) {
0159     getSummaries(it);
0160   }
0161   std::cout << endl;
0162   sortByGeometry();
0163   loadPresentValues();
0164
0165   // sanity check
0166   if(!checkInput()) return;
0167
0168   // check if both ISHA and VFS have to be tuned
0169   std::cout << "Preparing analysis..." << std::endl;
0170   int minISHA = 1000;
0171   int maxISHA = 0;
0172   int minVFS  = 1000;
0173   int maxVFS  = 0;
0174   for(FileList::const_iterator file=files_.begin();file!=files_.end();++file){
0175     int isha = file->first.first;
0176     int vfs  = file->first.second;
0177     minISHA = minISHA<isha ? minISHA : isha;
0178     maxISHA = maxISHA>isha ? maxISHA : isha;
0179     minVFS  = minVFS <vfs  ? minVFS  : vfs ;
0180     maxVFS  = maxVFS >vfs  ? maxVFS  : vfs ;
0181   }
0182   tuneISHA_ &= (minISHA!=maxISHA);
0183   tuneVFS_  &= (minVFS !=maxVFS );
0184   if(!tuneISHA_) std::cout << "ISHA tune disabled" << std::endl;
0185   if(!tuneVFS_ ) std::cout << "VFS  tune disabled" << std::endl;
0186   // two cases are possible:
0187   // ISHA tune: look at the rise time
0188   // VFS  tune: look at the tail
0189
0190   // number of APVs
0191   unsigned int nAPVs = (*(summaries_.begin()->second.begin()))->GetNbinsX();
0192
0193   // loop over the inputs to find individual values of ISHA ans VFS
0194   std::list<unsigned int> ishaValues;
0195   std::list<unsigned int> vfsValues;
0196   for(SummaryV::const_iterator summary = summaries_.begin(); summary!=summaries_.end(); ++summary) {
0197      ishaValues.push_back(summary->first.first);
0198      vfsValues.push_back(summary->first.second);
0199   }
0200   ishaValues.sort();
0201   vfsValues.sort();
0202   ishaValues.unique();
0203   vfsValues.unique();
0204
0205   // loop over apvs (bins)
0206   std::cout << "Running analysis..." << std::endl;
0207   for(unsigned int apv=1;apv<=nAPVs;++apv) {
0208      TGraph* g1 = new TGraph();
0209      TGraph* g2 = new TGraph();
0210      int ii=0;
0211      cout << "\r" << setw(5) << setfill('0') << apv << flush;
0212
0213      // loop over the VFS values
0214      for(std::list<unsigned int>::const_iterator vfs = vfsValues.begin(); vfs!=vfsValues.end(); ++vfs,++ii) {
0215        float tail = 0.;
0216        unsigned int npts = 0;
0217        for(SummaryV::const_iterator summary = summaries_.begin(); summary!=summaries_.end(); ++summary){
0218          if((unsigned int)summary->first.second==(*vfs)) {
0219            // determine which histogram are the rise time and the tail
0220            const std::vector<TH1*>& observables = summary->second;
0221            int tail_index = 0;
0222            int rise_index = 0;
0223            for( std::vector<TH1*>::const_iterator histo = observables.begin();histo<observables.end();++histo) {
0224               std::string name = (*histo)->GetName();
0225               if(name.find("CalibrationTail")!=std::string::npos) tail_index = histo-observables.begin();
0226               if(name.find("CalibrationRiseTime")!=std::string::npos) rise_index = histo-observables.begin();
0227            }
0228        //for vfs, we take the mean tail over the ISHA values at that point
0229        tail += observables[tail_index]->GetBinContent(apv);
0230        ++npts;
0231      }
0232        }
0233        // fill the graph
0234        g2->SetPoint(ii,(*vfs), tail/npts);
0235      }
0236 #ifdef DEBUG_ON
0237      std::string name2 = Form("graph%s%s",summaries_.begin()->second[0]->GetXaxis()->GetBinLabel(apv),"CalibrationTail");
0238      std::replace( name2.begin(), name2.end(), '.', '_' );
0239      g2->Write(name2.c_str());
0240 #endif
0241      // analyse the graphs
0242      float best_vfs  = tuneVFS_  ? getX(g2,50) :
0243                                    presentValues_[summaries_.begin()->second[0]->GetXaxis()->GetBinLabel(apv)].second;
0244      // now that VFS is optimized, take the ISHA values for the closest VFS point
0245      // for ISHA, we consider the rise time for VFS values close to the optimal
0246
0247      // find the closest point in the VFS scan
0248      float dist = 1000.;
0249      std::list<unsigned int>::const_iterator vfsPoint = vfsValues.begin();
0250      for(std::list<unsigned int>::const_iterator vfs = vfsValues.begin(); vfs!=vfsValues.end(); ++vfs) {
0251        if(dist>fabs((*vfs)-best_vfs)) {
0252          dist = fabs((*vfs)-best_vfs);
0253      vfsPoint = vfs;
0254        }
0255      }
0256      // loop over the ISHA values
0257      ii=0;
0258      for(std::list<unsigned int>::const_iterator isha = ishaValues.begin(); isha!=ishaValues.end(); ++isha,++ii) {
0259        for(SummaryV::const_iterator summary = summaries_.begin(); summary!=summaries_.end(); ++summary){
0260          if(((unsigned int)summary->first.second==(*vfsPoint))&&((unsigned int)summary->first.first==(*isha))) {
0261            // determine which histogram are the rise time and the tail
0262            const std::vector<TH1*>& observables = summary->second;
0263            int tail_index = 0;
0264            int rise_index = 0;
0265            for( std::vector<TH1*>::const_iterator histo = observables.begin();histo<observables.end();++histo) {
0266               std::string name = (*histo)->GetName();
0267               if(name.find("CalibrationTail")!=std::string::npos) tail_index = histo-observables.begin();
0268               if(name.find("CalibrationRiseTime")!=std::string::npos) rise_index = histo-observables.begin();
0269            }
0270        // fill the graph
0271        g1->SetPoint(ii,summary->first.first,observables[rise_index]->GetBinContent(apv));
0272 #ifdef DEBUG_ON
0273            std::string name1 = Form("graph%s%s",summaries_.begin()->second[0]->GetXaxis()->GetBinLabel(apv),"CalibrationRiseTime");
0274            std::replace( name1.begin(), name1.end(), '.', '_' );
0275            g1->Write(name1.c_str());
0276 #endif
0277      }
0278        }
0279      }
0280      // analyse the graphs
0281      float best_isha = tuneISHA_ ? getX(g1,53.5 ) :
0282                                    presentValues_[summaries_.begin()->second[0]->GetXaxis()->GetBinLabel(apv)].first;
0283
0284      // save the result
0285      result_[summaries_.begin()->second[0]->GetXaxis()->GetBinLabel(apv)] = std::make_pair((int)round(best_isha),(int)round(best_vfs));
0286
0287      // cleaning
0288      delete g1;
0289      delete g2;
0290   }
0291   std::cout << std::endl;
0292
0293 #ifdef DEBUG_ON
0294   debugFile->Write();
0295   debugFile->Close();
0296   delete debugFile;
0297 #endif
0298
0299 }
0300
0301 bool CalibrationScanAnalysis::checkInput() const {
0302
0303   // check that we have data
0304   std::cout << "Checking data integrity." << std::endl;
0305   std::cout << "Step 1/5" << std::endl;
0306   if(!summaries_.size()) {
0307     std::cerr << "Error: No summary histogram found." << std::endl
0308               << " Did you load any file ? " << std::endl;
0309     return 0;
0310   }
0311
0312   if(summaries_.size()<2) {
0313     std::cerr << "Error: Only one  summary histogram found." << std::endl
0314               << " Analysis does not make sense with only one measurement" << std::endl;
0315     return 0;
0316   }
0317
0318   // check that we have the same entries in each record,
0319   // check that the binning is the same in all histograms
0320   std::cout << "Step 2/5" << std::endl;
0321   int nbinsAll = -1;
0322   std::vector<std::string> namesAll;
0323   for(SummaryV::const_iterator summary = summaries_.begin(); summary!=summaries_.end(); ++summary) {
0324     const std::vector<TH1*>& observables = summary->second;
0325     for( std::vector<TH1*>::const_iterator histo = observables.begin();histo<observables.end();++histo) {
0326        std::string name = (*histo)->GetName();
0327        if(summary == summaries_.begin()) {
0328           namesAll.push_back(name);
0329        } else {
0330           if(find(namesAll.begin(),namesAll.end(),name)==namesAll.end()) {
0331             std::cerr << "Error: Found an histogram that is not common to all inputs: "
0332                       << name << std::endl;
0333             return 0;
0334           }
0335        }
0336        int nbins = (*histo)->GetNbinsX();
0337        if(nbinsAll<0) nbinsAll = nbins;
0338        if(nbins != nbinsAll) {
0339          std::cerr << "Error: The number of bins is not the same in all inputs." << std::endl;
0340 // non fatal
0341 //         return 0;
0342        }
0343     }
0344   }
0345
0346   // check that we have at least 2 histograms with measurements
0347   std::cout << "Step 3/5" << std::endl;
0348   if(namesAll.size()<2) {
0349     std::cerr << "Error: The number of available measurements is smaller than 2." << std::endl;
0350     return 0;
0351   }
0352
0353   // check that the bin labels are all the same
0354   std::cout << "Step 4/5" << std::endl;
0355   std::vector<std::string> labelsAll;
0356   for(SummaryV::const_iterator summary = summaries_.begin(); summary!=summaries_.end(); ++summary) {
0357     const std::vector<TH1*>& observables = summary->second;
0358     for( std::vector<TH1*>::const_iterator histo = observables.begin();histo<observables.end();++histo) {
0359        for(int i = 1;i <= (*histo)->GetNbinsX(); ++i) {
0360          std::string label = (*histo)->GetXaxis()->GetBinLabel(i);
0361          if(summary == summaries_.begin() && histo == observables.begin()) {
0362            labelsAll.push_back(label);
0363          } else {
0364            if(labelsAll[i-1] != label) {
0365 *((TH1F*)(*histo)) = TH1F(*(fixHisto(labelsAll,*histo)));
0366 /*
0367              std::cerr << "Error: Incoherency in bin labels. Bin " << i
0368                        << " of " << (*histo)->GetName() << " is " << label
0369                        << " and not " << labelsAll[i] << "." << std::endl;
0370              return 0;
0371 */
0372            }
0373          }
0374        }
0375     }
0376   }
0377
0378   // check that all APVs have an associated geometry
0379   std::cout << "Step 5/5" << std::endl;
0380    for(std::vector<std::string>::const_iterator apvLabel = labelsAll.begin();
0381        apvLabel != labelsAll.end(); ++apvLabel) {
0382      if(geometries_.find(*apvLabel)==geometries_.end()) {
0383        std::cerr << "Error: Geometry unknown for APV " << *apvLabel << std::endl;
0384        // made this a non-fatal error
0385  //      return 0;
0386        std::string label = *apvLabel;
0387        ((CalibrationScanAnalysis*)this)->geometries_[label] = 0;
0388      }
0389    }
0390
0391   return 1;
0392
0393 }
0394
0395 void CalibrationScanAnalysis::sortByGeometry() {
0396
0397   //categorize APVs per module geometry
0398   std::cout << "Reading cabling from debug.log" << std::endl;
0399   ifstream debuglog("debug.log");
0400   char buffer[1024];
0401   while(debuglog.getline(buffer,1024)) {
0402     if(strncmp(buffer," FED:cr/sl/id/fe/ch/chan",23)==0) {
0403
0404       // Decode input
0405       int fecCrate,fecSlot,fecRing,ccuAddr,ccuChan,channel1,channel2,detid,tmp;
0406       sscanf(strstr(buffer,"FEC:cr/sl/ring/ccu/mod"), "FEC:cr/sl/ring/ccu/mod=%d/%d/%d/%d/%d", &fecCrate,&fecSlot,&fecRing,&ccuAddr, &ccuChan);
0407       sscanf(strstr(buffer,"apvs"), "apvs=%d/%d", &channel1,&channel2);
0408       sscanf(strstr(buffer,"dcu/detid"), "dcu/detid=%x/%x", &tmp,&detid);
0409
0410       // Construct bin label
0411       std::stringstream bin1;
0412       bin1        << std::setw(1) << std::setfill('0') << fecCrate;
0413       bin1 << "." << std::setw(2) << std::setfill('0') << fecSlot;
0414       bin1 << "." << std::setw(1) << std::setfill('0') << fecRing;
0415       bin1 << "." << std::setw(3) << std::setfill('0') << ccuAddr;
0416       bin1 << "." << std::setw(2) << std::setfill('0') << ccuChan;
0417       bin1  << "." << channel1;
0418       std::stringstream bin2;
0419       bin2        << std::setw(1) << std::setfill('0') << fecCrate;
0420       bin2 << "." << std::setw(2) << std::setfill('0') << fecSlot;
0421       bin2 << "." << std::setw(1) << std::setfill('0') << fecRing;
0422       bin2 << "." << std::setw(3) << std::setfill('0') << ccuAddr;
0423       bin2 << "." << std::setw(2) << std::setfill('0') << ccuChan;
0424       bin2 << "." << channel2;
0425
0426       // Decode the detid -> sensor geometry
0427       int subdet = (detid>>25)&0x7;
0428       int ring = 0;
0429       if(subdet == 6) ring = (detid>>5)&0x7;
0430       if(subdet == 4) ring = (detid>>9)&0x3;
0431       int geom = ring + ((subdet==6) ? 3 : 0);
0432       if(subdet==3) geom +=10 ;
0433       if(subdet==5) geom +=15 ;
0434
0435       // Save
0436       geometries_[bin1.str()] = geom;
0437       geometries_[bin2.str()] = geom;
0438     }
0439   }
0440 }
0441
0442 void CalibrationScanAnalysis::loadPresentValues() {
0443
0444   //categorize APVs per module geometry
0445   std::cout << "Reading present ISHA/VFS values from debug.log" << std::endl;
0446   ifstream debuglog("debug.log");
0447   char buffer[1024];
0448   while(debuglog.getline(buffer,1024)) {
0449     if(strncmp(buffer,"Present values for ISHA/VFS",27)==0) {
0450       // Decode input
0451       int isha, vfs;
0452       char apv_addr[256];
0453       sscanf(strstr(buffer,"APV"),"APV %s : %d %d", apv_addr,&isha,&vfs);
0454       // Save
0455       std::string apv_address = apv_addr;
0456       presentValues_[apv_address] = std::make_pair(isha,vfs);
0457     }
0458   }
0459
0460 }
0461
0462 void CalibrationScanAnalysis::sanitizeResult(unsigned int cut, bool doItForISHA, bool doItForVFS) {
0463
0464   // create and fill the utility histograms (similar to the draw method)
0465   std::cout << "Applying sanity constraints on the results." << std::endl;
0466   std::map<int,TH2F*> histos;
0467   for(std::map<std::string, int>::iterator it = geometries_.begin(); it!= geometries_.end(); ++it) {
0468     if(histos.find(it->second)==histos.end()) {
0469       TH2F* histo = new TH2F(Form("modulesGeometry%d",it->second),
0470                              Form("Module Geometry %d",it->second),255,0,255,255,0,255);
0471       histos[it->second] = histo;
0472     }
0473   }
0474
0475   // first loop to compute mean and rms
0476   for(std::map<std::string, Parameters>::const_iterator apvValue = result_.begin();
0477       apvValue != result_.end(); ++apvValue) {
0478     histos[geometries_[apvValue->first]]->Fill(apvValue->second.first,apvValue->second.second);
0479   }
0480
0481   // second loop to cut at x RMS
0482   int lowVFS,highVFS,lowISHA,highISHA,geom;
0483   for(std::map<std::string, Parameters>::iterator apvValue = result_.begin();
0484       apvValue != result_.end(); ++apvValue) {
0485     geom = geometries_[apvValue->first];
0486     lowISHA  = (int)round(histos[geom]->GetMean(1) -
0487                           cut*histos[geom]->GetRMS(1));
0488     highISHA = (int)round(histos[geom]->GetMean(1) +
0489                           cut*histos[geom]->GetRMS(1));
0490     lowVFS   = (int)round(histos[geom]->GetMean(2) -
0491                           cut*histos[geom]->GetRMS(2));
0492     highVFS  = (int)round(histos[geom]->GetMean(2) +
0493                           cut*histos[geom]->GetRMS(2));
0494     if((apvValue->second.first<lowISHA || apvValue->second.first>highISHA) && doItForISHA) {
0495       apvValue->second.first = (int)round((lowISHA+highISHA)/2.);
0496     }
0497     if((apvValue->second.second<lowVFS || apvValue->second.second>highVFS) && doItForVFS) {
0498       apvValue->second.second = (int)round((lowVFS+highVFS)/2.);
0499     }
0500   }
0501
0502   // finaly delete the temporary histograms
0503   for(std::map<int,TH2F*>::iterator it = histos.begin(); it!=histos.end(); ++it) {
0504     delete it->second;
0505   }
0506
0507 }
0508
0509 void CalibrationScanAnalysis::print(Option_t*) const {
0510
0511   // input
0512   std::cout << "Analysis of ISHA/VFS performed using the following inputs: " << std::endl;
0513   std::cout << "ISHA \t VFS \t File" << std::endl;
0514   for(FileList::const_iterator file=files_.begin();file!=files_.end();++file){
0515     int isha = file->first.first;
0516     int vfs  = file->first.second;
0517     std::string filename = file->second->GetName();
0518     std::cout << isha << "\t " << vfs << "\t " << filename << std::endl;
0519   }
0520   std::cout << std::endl << std::endl;
0521
0522   // output
0523   std::cout << "Resulting values: " << std::endl;
0524   std::cout << "APV \t \t ISHA \t VFS" << std::endl;
0525   for(std::map<std::string, Parameters>::const_iterator result = result_.begin();
0526       result != result_.end(); ++result) {
0527     std::cout << result->first << "\t "
0528               << result->second.first << "\t "
0529               << result->second.second << std::endl;
0530   }
0531
0532 }
0533
0534 void CalibrationScanAnalysis::draw(Option_t*) const {
0535
0536   std::cout << "Drawing results..." << std::endl;
0537
0538   // first create the histograms
0539   std::cout << "   - first create the 2D histograms" << std::endl;
0540   new TCanvas;
0541   std::map<int,TH2F*> histos;
0542   for(std::map<std::string, int>::const_iterator it = geometries_.begin(); it!= geometries_.end(); ++it) {
0543     if(histos.find(it->second)==histos.end()) {
0544       TH2F* histo = new TH2F(Form("modulesGeometry%d",it->second),
0545                              Form("Module Geometry %d",it->second),255,-1.25,0.6625,255,0,255);
0546       histo->SetDirectory(0);
0547       histo->GetXaxis()->SetTitle("VFS");
0548       histo->GetYaxis()->SetTitle("ISHA");
0549       histo->SetMarkerStyle(7);
0550       histo->SetMarkerColor(2+it->second);
0551       histos[it->second] = histo;
0552     }
0553   }
0554
0555   // loop over apvs
0556   std::cout << "   - loop over apvs" << std::endl;
0557   for(std::map<std::string, Parameters>::const_iterator apvValue = result_.begin();
0558       apvValue != result_.end(); ++apvValue) {
0559     histos[geometries_.find(apvValue->first)->second]->Fill(-1.25+apvValue->second.second*0.0075,apvValue->second.first);
0560   }
0561
0562   // draw the histograms
0563   std::cout << "   - draw the histograms" << std::endl;
0564   for(std::map<int,TH2F*>::iterator h = histos.begin(); h != histos.end(); ++h) {
0565     h->second->Draw(h == histos.begin() ? "" : "same");
0566   }
0567
0568   // draw the histogram with the mean per geometry
0569   std::cout << "   - draw the histogram with the mean per geometry" << std::endl;
0570   TH2F* histo = new TH2F("Geometries","Geometries",255,-1.25,0.6625,255,0,255);
0571   histo->SetDirectory(0);
0572   histo->GetXaxis()->SetTitle("VFS");
0573   histo->GetYaxis()->SetTitle("ISHA");
0574   histo->SetMarkerStyle(20);
0575   histo->SetMarkerColor(2);
0576   for(std::map<int,TH2F*>::iterator h = histos.begin(); h != histos.end(); ++h) {
0577     histo->Fill(h->second->GetMean(1),h->second->GetMean(2));
0578   }
0579   histo->Draw("same");
0580
0581 //////////////////////////////////////////////////////////////
0582
0583
0584   // first create the histograms
0585   std::cout << "   - create the 1D histograms" << std::endl;
0586   new TCanvas;
0587   std::map<int,TH1F*> histosVFS;
0588   std::map<int,TH1F*> histosISHA;
0589   for(std::map<std::string, int>::const_iterator it = geometries_.begin(); it!= geometries_.end(); ++it) {
0590     if(histosVFS.find(it->second)==histosVFS.end()) {
0591       TH1F* histoVFS = new TH1F(Form("VFSmodulesGeometry%d",it->second),
0592                                 Form("VFS for Module Geometry %d",it->second),255,0,255);
0593       histoVFS->SetDirectory(0);
0594       histosVFS[it->second] = histoVFS;
0595       TH1F* histoISHA = new TH1F(Form("ISHAmodulesGeometry%d",it->second),
0596                                  Form("ISHA for Module Geometry %d",it->second),255,0,255);
0597       histoISHA->SetDirectory(0);
0598       histosISHA[it->second] = histoISHA;
0599     }
0600   }
0601
0602   // loop over apvs
0603   std::cout << "   - loop over apvs" << std::endl;
0604   for(std::map<std::string, Parameters>::const_iterator apvValue = result_.begin();
0605       apvValue != result_.end(); ++apvValue) {
0606     histosISHA[geometries_.find(apvValue->first)->second]->Fill(apvValue->second.first);
0607     histosVFS[geometries_.find(apvValue->first)->second]->Fill(apvValue->second.second);
0608   }
0609
0610   // draw the histograms
0611   std::cout << "   - draw the histograms" << std::endl;
0612   for(std::map<int,TH1F*>::iterator h = histosISHA.begin(); h != histosISHA.end(); ++h) {
0613     h->second->Draw(h == histosISHA.begin() ? "" : "same");
0614   }
0615
0616   new TCanvas;
0617
0618   for(std::map<int,TH1F*>::iterator h = histosVFS.begin(); h != histosVFS.end(); ++h) {
0619     h->second->Draw(h == histosVFS.begin() ? "" : "same");
0620   }
0621
0622 }
0623
0624 void CalibrationScanAnalysis::save(const char* fileName) {
0625
0626   std::cout << "Saving results..." << std::endl;
0627 /*
0628   // save in the input files
0629   for(FileList::const_iterator it=files_.begin();it!=files_.end();++it) {
0630     TFile* input = it->second;
0631     TDirectory* directory = input->GetDirectory(HISTOPATH);
0632     directory->cd();
0633     TList* histograms = directory->GetListOfKeys();
0634     TIter next(histograms);
0635     TKey* key = NULL;
0636     while ((key = (TKey*)next())) {
0637       if(TClass(key->GetClassName()).InheritsFrom("TH1")) {
0638         TH1* h = (TH1*)key->ReadObj();
0639         std::string name = h->GetName();
0640         if(name.find("CalibrationTail")!=std::string::npos) {
0641           TH1F* ishaOutput = new TH1F("isha","isha",h->GetNbinsX(),0,h->GetNbinsX());
0642           TH1F* vfsOutput = new TH1F("vfs","vfs",h->GetNbinsX(),0,h->GetNbinsX());
0643           for(int i=1;i<=h->GetNbinsX();++i) {
0644             ishaOutput->SetBinContent(i,result_[h->GetXaxis()->GetBinLabel(i)].first);
0645             vfsOutput ->SetBinContent(i,result_[h->GetXaxis()->GetBinLabel(i)].second);
0646             ishaOutput->GetXaxis()->SetBinLabel(i,h->GetXaxis()->GetBinLabel(i));
0647             vfsOutput ->GetXaxis()->SetBinLabel(i,h->GetXaxis()->GetBinLabel(i));
0648           }
0649           break;
0650         }
0651       }
0652     }
0653     input->Write();
0654   }
0655 */
0656   // save in a file for TkConfigurationDb
0657   std::ostream * output;
0658   TString filen = fileName;
0659   if (filen == "")
0660     return;
0661   if (filen == "-")
0662     output = &std::cout;
0663   else
0664     output = new std::ofstream(fileName);
0665   FileList::const_iterator it=files_.begin();
0666   TFile* input = it->second;
0667   TDirectory* directory = input->GetDirectory(HISTOPATH);
0668   directory->cd();
0669   TList* histograms = directory->GetListOfKeys();
0670   TIter next(histograms);
0671   TKey* key = NULL;
0672   while ((key = (TKey*)next())) {
0673     if(TClass(key->GetClassName()).InheritsFrom("TH1")) {
0674       TH1* h = (TH1*)key->ReadObj();
0675       std::string name = h->GetName();
0676       if(name.find("CalibrationTail")!=std::string::npos) {
0677         for(int i=1;i<=h->GetNbinsX();++i) {
0678       std::string address = h->GetXaxis()->GetBinLabel(i);
0679       address = address.substr(address.find('.')+1);
0680       std::replace(address.begin(),address.end(),'.',' ');
0681       *output << address << " "
0682               << result_[h->GetXaxis()->GetBinLabel(i)].first << " "
0683                   << result_[h->GetXaxis()->GetBinLabel(i)].second << std::endl;
0684         }
0685         break;
0686       }
0687     }
0688   }
0689 }
0690
0691 TH1F* CalibrationScanAnalysis::fixHisto(std::vector<std::string>& names,TH1* histo) const
0692 {
0693    // prepare an histogram to replace input
0694    TH1F* newHisto = new TH1F(histo->GetName(),histo->GetTitle(),names.size(),histo->GetXaxis()->GetXmin(),histo->GetXaxis()->GetXmax());
0695    std::cout << "fixing histo " << histo->GetName() << " at " << histo << " by " << newHisto << std::endl;
0696    for(std::vector<std::string>::iterator name=names.begin();name!=names.end();++name) {
0697      newHisto->GetXaxis()->SetBinLabel(name-names.begin()+1,name->c_str());
0698      int pos = histo->GetXaxis()->FindBin(name->c_str());
0699      if(pos!=-1) {
0700        newHisto->SetBinContent(pos,histo->GetBinContent(pos));
0701      }
0702    }
0703    return newHisto;
0704 }