OfflineValidation/scripts/merge_TrackerOfflineValidation.C

0001 /*
0002
0003 NOTE: This macro is based on the hadd macro provided by root.
0004 The merging of a TTree is different from what hadd does to merge trees.
0005 The Tree is once read out for entries which are the same in all input
0006 variables - like global position variables for each module - and the values
0007 are copied into a new Ttree.
0008 The mean and RMS values for the different residuals are taken from the merged
0009 histograms and filled into the tree after the merging procedure.
0010
0011 File names can be given in two ways:
0012 1)
0013 Provide a comma separated list of filenames as argument to hadd.
0014 2)
0015 Use default empty argument. You will be asked for common file name 'prefix'
0016 and the number of files to be merged.
0017 File names will then be generated as
0018 <prefix>1.root, <prefix>2.root, <prefix>3.root, etc.
0019 until the number specified.
0020
0021 NOTE:
0022 Includes from CMSSW release as needed below need CMSSW environment as
0023 you can get by these lines (add them to your rootlogon.C):
0024
0025 #include "TSystem.h"
0026 if (gSystem->Getenv("CMSSW_RELEASE_BASE") != '\0') {
0027   printf("\nLoading CMSSW FWLite...\n");
0028   gSystem->Load("libFWCoreFWLite");
0029   FWLiteEnabler::enable();
0030 }
0031
0032 */
0033
0034 // This line works only if we have a CMSSW environment...
0035 #include "Alignment/OfflineValidation/interface/TkOffTreeVariables.h"
0036
0037 #include "TROOT.h"
0038 #include <string>
0039 #include <map>
0040 #include <iostream>
0041 #include "TChain.h"
0042 #include "TFile.h"
0043 #include "TH1.h"
0044 #include "TTree.h"
0045 #include "TKey.h"
0046 #include "TF1.h"
0047 #include "TString.h"
0048 #include "TObjString.h"
0049 #include "TMath.h"
0050 // global variables:
0051 bool copiedTree_ = false;
0052 std::map<unsigned int, TkOffTreeVariables> map_;
0053 std::map<unsigned int, TString> idPathMap_;
0054
0055 // methods
0056 void MergeRootfile( TDirectory *target, TList *sourcelist);
0057 void RewriteTree( TDirectory *target,TTree *tree);
0058 std::pair<float,float> FitResiduals(TH1 *h, float meantmp, float rmstmp);
0059 float getMedian(const TH1 *histo);
0060 //////////////////////////////////////////////////////////////////////////
0061 // master method
0062 //////////////////////////////////////////////////////////////////////////
0063 int hadd(const char *filesSeparatedByCommaOrEmpty = "", const char * outputFile = "") {
0064 //void merge_TrackerOfflineValidation(const char *filesSeparatedByCommaOrEmpty = "") {
0065
0066   TString fileNames(filesSeparatedByCommaOrEmpty);
0067   TObjArray *names = 0;
0068
0069   if (fileNames.Length() == 0) {
0070     // nothing specified, so read in
0071     TString inputFileName="";
0072     Int_t nEvt = 0;
0073     //get file name
0074     std::cout<<"Type in general file name (ending is added automatically)."<<std::endl;
0075     std::cout<<"i.e. Validation_CRUZET_data_ as name, code adds 1,2,3 + .root"<<std::endl;
0076     std::cin>>inputFileName;
0077
0078     std::cout << "Type number of files to merge." << std::endl;
0079     std::cin>>nEvt;
0080
0081     names = new TObjArray;
0082     names->SetOwner();
0083     for (Int_t i = 1; i != nEvt+1 ; ++i){
0084       TString fileName = inputFileName;
0085       names->Add(new TObjString((fileName += i) += ".root"));
0086     }
0087   } else {
0088     // decode file name from comma separated list
0089     names = fileNames.Tokenize(","); // is already owner
0090   }
0091
0092   int iFilesUsed=0, iFilesSkipped=0;
0093   const TString keyName = "TrackerOfflineValidationStandalone";
0094   TList *FileList = new TList();
0095   for (Int_t iFile = 0; iFile < names->GetEntriesFast(); ++iFile) {
0096     TFile *file = TFile::Open(names->At(iFile)->GetName());
0097     if (file) {
0098       // check if file is ok and contains data
0099       if (file->FindKey(keyName)) {
0100     FileList->Add(file);
0101     std::cout << names->At(iFile)->GetName() << std::endl;
0102     iFilesUsed++;
0103       }
0104       else {
0105     std::cout << names->At(iFile)->GetName()
0106           << " --- does not contain data, skipping file " << std::endl;
0107     iFilesSkipped++;
0108       }
0109     } else {
0110       cout << "File " << names->At(iFile)->GetName() << " does not exist!" << endl;
0111       delete names; names = 0;
0112       return 1;
0113     }
0114   }
0115   delete names;
0116
0117   TString outputFileString;
0118
0119   if (strlen(outputFile)!=0)
0120     outputFileString = TString(outputFile);
0121   else
0122     outputFileString = "merge_output.root";
0123   TFile *Target = TFile::Open( outputFileString, "RECREATE" );
0124   MergeRootfile( Target, FileList);
0125   std::cout << "Finished merging of histograms." << std::endl;
0126
0127   Target->cd( keyName );
0128   TTree *tree = new TTree("TkOffVal","TkOffVal");
0129   RewriteTree(Target,tree);
0130   tree->Write();
0131   std::cout << std::endl;
0132   std::cout << "Read " << iFilesUsed << " files, " << iFilesSkipped << " files skipped" << std::endl;
0133   if (iFilesSkipped>0)
0134     std::cout << " (maybe because there was no data in those files?)" << std::endl;
0135   std::cout << "Written Tree, now saving hists..." << std::endl;
0136
0137   //  Target->SaveSelf(kTRUE);
0138   std::cout << "Closing file " << Target->GetName() << std::endl;
0139
0140   // There is a bug in this file(?) and the script gets
0141   // stuck in the line "delete Target"
0142   // without the followin
0143   std::cout << endl;
0144   std::cout << "-----------------------------------------------------------------" << std::endl;
0145   std::cout << "---  Because of a bug in the code, ROOT gets stuck when  --------" << std::endl;
0146   std::cout << "---  closing the file.  In that case, please do the following ---" << std::endl;
0147   std::cout << "--- 1) hit CTRL-Z -----------------------------------------------" << std::endl;
0148   std::cout << "--- 2) check the ID number or the process root with ps---------- " << std::endl;
0149   std::cout << "--- 3) kill the root process ------------------------------------" << std::endl;
0150   std::cout << "--- 4) continue plotting with fg (only a few minutes left)------ " << std::endl;
0151   std::cout << "-- ------------------------------------------------------------- " << std::endl;
0152
0153
0154   //  gROOT->GetListOfFiles()->Remove(Target);
0155
0156   // tested, does not work
0157   //  Target->cd();
0158   Target->Close();
0159   //delete Target;
0160
0161   // The abort() command is ugly, but much quicker than the clean return()
0162   // Use of return() can take 90 minutes, while abort() takes 10 minutes
0163   // (merging 20 jobs with 1M events in total)
0164   // BUT abort creates an error signal and incorrect functioning
0165   // at higher level
0166   // abort();
0167
0168   std::cout << "Now returning from merge_TrackerOfflineValidation.C" << std::endl;
0169   return 0;
0170 }
0171
0172 /////////////////////////////////////////////////////////////////////////
0173 /////////////////////////////////////////////////////////////////////////
0174 /////////////////////////////////////////////////////////////////////////
0175 void MergeRootfile( TDirectory *target, TList *sourcelist) {
0176
0177   cout << target->GetPath() << endl;
0178   TString path( (char*)strstr( target->GetPath(), ".root:" ) );
0179   path.Remove( 0, 7 );
0180   //  TString tmp;
0181   TFile *first_source = (TFile*)sourcelist->First();
0182   first_source->cd( path );
0183   TDirectory *current_sourcedir = gDirectory;
0184   //gain time, do not add the objects in the list in memory
0185   Bool_t status = TH1::AddDirectoryStatus();
0186   TH1::AddDirectory(kFALSE);
0187
0188   // loop over all keys in this directory
0189
0190   TIter nextkey( current_sourcedir->GetListOfKeys() );
0191   TKey *key, *oldkey=0;
0192   while ( (key = (TKey*)nextkey())) {
0193
0194     //keep only the highest cycle number for each key
0195     if (oldkey && !strcmp(oldkey->GetName(),key->GetName())) continue;
0196
0197     // read object from first source file
0198     first_source->cd( path );
0199     TObject *obj = key->ReadObj();
0200
0201     if ( obj->IsA()->InheritsFrom( "TH1" ) ) {
0202       // descendant of TH1 -> merge it
0203
0204       // But first note the path for this moduleId!
0205       const TString histName(obj->GetName());
0206       if (histName.Contains("module_")){
0207     // GF: fragile! have to look after 'module_'!
0208     const TString tmp = histName(histName.Last('_')+1, 9);
0209     unsigned int moduleId = tmp.Atoi(); // check for valid ID???
0210     idPathMap_[moduleId] = path + '/'; // only for first hist...?
0211       }
0212
0213       TH1 *h1 = (TH1*)obj;
0214
0215       // loop over all source files and add the content of the
0216       // correspondant histogram to the one pointed to by "h1"
0217       TFile *nextsource = (TFile*)sourcelist->After( first_source );
0218       while ( nextsource ) {
0219
0220         // make sure we are at the correct directory level by cd'ing to path
0221         nextsource->cd( path );
0222         TKey *key2 = (TKey*)gDirectory->GetListOfKeys()->FindObject(h1->GetName());
0223         if (key2) {
0224       TH1 *h2 = (TH1*)key2->ReadObj();
0225       h1->Add( h2 );
0226       delete h2;
0227         }
0228
0229         nextsource = (TFile*)sourcelist->After( nextsource );
0230       }
0231     }
0232     else if ( obj->IsA()->InheritsFrom( "TTree" ) ) {
0233       if (!copiedTree_ && TString("TkOffVal") == obj->GetName()) {
0234     // get tree structure and 'const' entries for each module once
0235     // ('constant' means not effected by merging)
0236     TTree* tree =(TTree*)obj;
0237
0238     TkOffTreeVariables *treeMem = 0; // ROOT will initilise
0239     tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
0240
0241     std::cout << "Copy info from first TTree." << std::endl;
0242     Long64_t nentries = tree->GetEntriesFast();
0243     for (Long64_t jentry = 0; jentry < nentries; ++jentry) {
0244       if (jentry%1000==0)
0245         std::cout << "Copy entry " << jentry << " / " << nentries << std::endl;
0246       tree->GetEntry(jentry);
0247       // Erase everything not common:
0248       treeMem->clearMergeAffectedPart();
0249       map_[treeMem->moduleId] = *treeMem;
0250     }
0251     std::cout << "Done copy info from first TTree." << std::endl;
0252     copiedTree_ = true;
0253       }
0254     } else if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0255       // it's a subdirectory
0256       //      cout << "Found subdirectory " << obj->GetName() << endl;
0257
0258       // create a new subdir of same name and title in the target file
0259       target->cd();
0260       TDirectory *newdir = target->mkdir( obj->GetName(), obj->GetTitle() );
0261
0262       // newdir is now the starting point of another round of merging
0263       // newdir still knows its depth within the target file via
0264       // GetPath(), so we can still figure out where we are in the recursion
0265       MergeRootfile( newdir, sourcelist );
0266
0267     } else {
0268
0269       // object is of no type that we know or can handle
0270       cout << "Unknown object type, name: "
0271            << obj->GetName() << " title: " << obj->GetTitle() << endl;
0272     }
0273
0274     // now write the merged histogram (which is "in" obj) to the target file
0275     // note that this will just store obj in the current directory level,
0276     // which is not persistent until the complete directory itself is stored
0277     // by "target->Write()" below
0278     if ( obj ) {
0279       target->cd();
0280
0281       if(obj->IsA()->InheritsFrom( "TTree" ))cout<<"TTree will be merged separately"<<endl;
0282       else
0283     obj->Write( key->GetName() );
0284
0285     }
0286
0287   } // while ( ( TKey *key = (TKey*)nextkey() ) )
0288
0289   // save modifications to target file
0290   target->SaveSelf(kTRUE);
0291   TH1::AddDirectory(status);
0292
0293   return;
0294 }
0295
0296 //////////////////////////////////////////////////////////////////////////////////
0297 void RewriteTree( TDirectory *target, TTree *tree)
0298
0299 {
0300   TkOffTreeVariables *treeVar = new TkOffTreeVariables;
0301   tree->Branch("TkOffTreeVariables", &treeVar);
0302
0303   // As above: gain time, do not add the objects in the list in memory.
0304   // In addition avoids saving hists twice...
0305   const Bool_t status = TH1::AddDirectoryStatus();
0306   TH1::AddDirectory(kFALSE);
0307
0308   std::cout << "Create merged TTree:" << std::endl;
0309   TH1 *h = 0;
0310   unsigned int counter = 0;
0311   for (std::map<unsigned int,TkOffTreeVariables>::const_iterator it = map_.begin(),
0312      itEnd = map_.end(); it != itEnd; ++it, ++counter) {
0313     treeVar->clear();
0314     // first get 'constant' values from map ('constant' means not effected by merging)
0315     (*treeVar) = it->second; // (includes module id)
0316
0317     // now path name:
0318     const TString &path = idPathMap_[treeVar->moduleId];
0319
0320     // Read out mean and rms values from merged histograms,
0321     // fit for fitted values etc.
0322
0323     //////////////////////////////////////////////////////////////
0324     // xPrime
0325     //////////////////////////////////////////////////////////////
0326     target->GetObject(path + treeVar->histNameX, h);
0327     if (h) {
0328       treeVar->entries = static_cast<UInt_t>(h->GetEntries());   //entries are const for all histos
0329       treeVar->meanX = h->GetMean(); //get mean value from histogram
0330       treeVar->rmsX = h->GetRMS();   //get RMS value from histogram
0331
0332       int numberOfBins = h->GetNbinsX();
0333       treeVar->numberOfUnderflows = h->GetBinContent(0);
0334       treeVar->numberOfOverflows = h->GetBinContent(numberOfBins+1);
0335       treeVar->numberOfOutliers =  h->GetBinContent(0)+h->GetBinContent(numberOfBins+1);
0336
0337       const std::pair<float,float> meanSigma = FitResiduals(h, h->GetMean(), h->GetRMS());
0338       treeVar->fitMeanX = meanSigma.first;
0339       treeVar->fitSigmaX= meanSigma.second;
0340       treeVar->medianX = getMedian(h);
0341       delete h; h = 0;
0342     } else {
0343       std::cout << "Module " << treeVar->moduleId << " without hist X: "
0344             << path << treeVar->histNameX << std::endl;
0345     }
0346
0347     ////////////////////////////////////////////////////////////////
0348     // normalized xPrime
0349     ////////////////////////////////////////////////////////////////
0350     target->GetObject(path + treeVar->histNameNormX, h);
0351     if (h) {
0352       treeVar->meanNormX = h->GetMean(); //get mean value from histogram
0353       treeVar->rmsNormX = h->GetRMS();   //get RMS value from histogram
0354
0355       const std::pair<float,float> meanSigma = FitResiduals(h, h->GetMean(), h->GetRMS());
0356       treeVar->fitMeanNormX  = meanSigma.first; // get mean value from histogram
0357       treeVar->fitSigmaNormX = meanSigma.second; // get sigma value from histogram
0358
0359       double stats[20];
0360       h->GetStats(stats);
0361       //    treeVar->chi2PerDofX = stats[3]/(stats[0]-1); // GF: why -1?
0362       if (stats[0]) treeVar->chi2PerDofX = stats[3]/stats[0];
0363       delete h; h = 0;
0364     } else {
0365       std::cout << "Module " << treeVar->moduleId << " without hist normX: "
0366         << path << treeVar->histNameNormX << std::endl;
0367     }
0368
0369     ////////////////////////////////////////////////////////////////
0370     // local x if present
0371     ////////////////////////////////////////////////////////////////
0372     if (treeVar->histNameLocalX.size()) {
0373       target->GetObject(path + treeVar->histNameLocalX, h);
0374       if (h) {
0375     treeVar->meanLocalX = h->GetMean();      //get mean value from histogram
0376     treeVar->rmsLocalX = h->GetRMS();        //get RMS value from histogram
0377
0378     delete h; h = 0;
0379       } else {
0380     std::cout << "Module " << treeVar->moduleId << " without hist local X: "
0381           << path << treeVar->histNameLocalX << std::endl;
0382       }
0383     }
0384
0385     ////////////////////////////////////////////////////////////////
0386     // local normalized x if present
0387     ////////////////////////////////////////////////////////////////
0388     if (treeVar->histNameNormLocalX.size()) {
0389       target->GetObject(path + treeVar->histNameNormLocalX, h);
0390       if (h) {
0391     treeVar->meanNormLocalX = h->GetMean();//get mean value from histogram
0392     treeVar->rmsNormLocalX = h->GetRMS();  //get RMS value from histogram
0393     delete h; h = 0;
0394       } else {
0395     std::cout << "Module " << treeVar->moduleId << " without hist normLocal X"
0396           << path << treeVar->histNameNormLocalX << std::endl;
0397       }
0398     }
0399
0400     ////////////////////////////////////////////////////////////////
0401     // yPrime if existing (pixel and, if configured, for strip)
0402     ////////////////////////////////////////////////////////////////
0403     if (treeVar->histNameY.size()) {
0404       target->GetObject(path + treeVar->histNameY, h);
0405       if (h) {
0406     treeVar->meanY = h->GetMean();      //get mean value from histogram
0407     treeVar->rmsY = h->GetRMS();        //get RMS value from histogram
0408
0409     const std::pair<float,float> meanSigma = FitResiduals(h, h->GetMean(), h->GetRMS());
0410     treeVar->fitMeanY = meanSigma.first;
0411     treeVar->fitSigmaY= meanSigma.second;
0412     treeVar->medianY = getMedian(h);
0413     delete h; h = 0;
0414       } else {
0415     std::cout << "Module " << treeVar->moduleId << " without hist Y "
0416           << path << treeVar->histNameY << std::endl;
0417       }
0418     }
0419
0420     ////////////////////////////////////////////////////////////////
0421     // normalized yPrime
0422     ////////////////////////////////////////////////////////////////
0423     if (treeVar->histNameNormY.size()) {
0424       target->GetObject(path + treeVar->histNameNormY, h);
0425       if (h) {
0426     treeVar->meanNormY = h->GetMean(); //get mean value from histogram
0427     treeVar->rmsNormY = h->GetRMS();   //get RMS value from histogram
0428
0429     const std::pair<float,float> meanSigma = FitResiduals(h, h->GetMean(), h->GetRMS());
0430     treeVar->fitMeanNormY  = meanSigma.first; // get mean value from histogram
0431     treeVar->fitSigmaNormY = meanSigma.second; // get sigma value from histogram
0432
0433     double stats[20];
0434     h->GetStats(stats);
0435     if (stats[0]) treeVar->chi2PerDofY = stats[3]/stats[0];
0436     delete h; h = 0;
0437       } else {
0438     std::cout << "Module " << treeVar->moduleId << " without hist norm Y "
0439           << path << treeVar->histNameNormY << std::endl;
0440       }
0441     }
0442
0443     tree->Fill();
0444   } // end loop on modules
0445
0446   TH1::AddDirectory(status); // back to where we were before
0447 }
0448
0449 ///////////////////////////////////////////////////////////////////////////
0450 std::pair<float,float> FitResiduals(TH1 *h,float meantmp,float rmstmp)
0451 {
0452   std::pair<float,float> fitResult(9999., 9999.);
0453   if (!h || h->GetEntries() < 20) return fitResult;
0454
0455   // first fit: two RMS around mean
0456   TF1 func("tmp", "gaus", meantmp - 2.*rmstmp, meantmp + 2.*rmstmp);
0457   if (0 == h->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
0458     float mean  = func.GetParameter(1);
0459     float sigma = func.GetParameter(2);
0460
0461     // second fit: three sigma of first fit around mean of first fit
0462     func.SetRange(mean - 3.*sigma, mean + 3.*sigma);
0463     // I: integral gives more correct results if binning is too wide (slow)
0464     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
0465     if (0 == h->Fit(&func, "QNRL")) { // I")) {
0466       fitResult.first = func.GetParameter(1);
0467       fitResult.second = func.GetParameter(2);
0468     }
0469   }
0470
0471   return fitResult;
0472 }
0473
0474 float getMedian(const TH1 *histo)
0475 {
0476
0477   float median = 999;
0478   int nbins = histo->GetNbinsX();
0479
0480  //extract median from histogram
0481
0482    double *x = new double[nbins];
0483   double *y = new double[nbins];
0484   for (int j = 0; j < nbins; j++) {
0485     x[j] = histo->GetBinCenter(j+1);
0486     y[j] = histo->GetBinContent(j+1);
0487   }
0488   median = TMath::Median(nbins, x, y);
0489
0490
0491   delete[] x; x = 0;
0492   delete[] y; y = 0;
0493
0494   return median;
0495
0496 }