tschum/FWlite_Analysis/PlotTool.cc

#include "PlotTool.h"
//------------------------------------------------------------------
//Constructur: Set some global Variables to define canvas output
PlotTool::PlotTool() {

        this->SetClass("TChain", 100);

        samePad_trees = true;
        samePad_vars = false;
        samePad_cuts = false;
        sameCanv_trees = false;
        sameCanv_vars = false;
        sameCanv_cuts = false;

        showLegend = false;
        logY = true;
        addTrackJets = true;
        addHitDetInfo = true;
        diTrackMass = true;
        addTrigger = true;
        addTrackJets = false;
        addEventInfo = false;
        addTower = false;
        verbose = true;

        globalCuts="";

        recreateTree = false; //force recreation of friend tree

}

//------------------------------------------------------------------
//Fill object PlotTool with Chains constructed from files from given source 

int PlotTool::init(string fileName, string dirPath, string treeName,
                   string fileLabel) {
  this->New(this->GetEntries() );
  int currChain = this->GetEntries() - 1;

  if (currChain < 0)
    return currChain;

  TStopwatch timer;
  timer.Start();


  fileNames.clear();

  //check if alternative label is different from default(searchString)
  if (fileLabel=="") {
    fileLabel=fileName;
  }
  ((TChain*) this->At(currChain))->SetName(fileLabel.c_str());

  TList *files = new TList();
  TSystemFile* sysFile = 0;
  if(fileName.find(".") != string::npos) {
    ifstream f(fileName.c_str());
    if( ! f.is_open() ) return -1;
    string line;

    while (!f.eof()) {
      getline(f,line);
      sysFile = new TSystemFile(line.c_str(),dirPath.c_str());
      files->Add(sysFile);
    }


  } else {
    TSystemDirectory dir("sourceDir", dirPath.c_str());
    files = dir.GetListOfFiles();
  }
  if (files->GetEntries()>0) {
    TIter next(files);
    TSystemFile *file;
    TString fname;
    string filePath;

    if(verbose && fileName.find(".") == string::npos )  cout<<"Open"<<dirPath.c_str()<<" Search for .root files that contain: "
                                                            <<fileName.c_str()<<endl;
    if(verbose && fileName.find(".") != string::npos )  cout<<"Open"<<dirPath.c_str()<<" Search lines with .root in: "
                                                            <<fileName.c_str()<<endl;


    while ((file=(TSystemFile*)next())) {
      fname = file->GetName();
      if (!fname.EndsWith(".root"))
        continue;
      if (!fname.Contains(fileName.c_str()) && fileName.find(".") == string::npos )
        continue;

      filePath = dirPath;
      filePath += fname.Data();
      //fwlite::ChainEvent to lop over events, jets, etc
      fileNames.push_back(filePath.c_str());

      if (((TChain*) this->At(currChain))->AddFile(filePath.c_str(), -1,
                                                   treeName.c_str()) )
        if(verbose)     cout<<"Chained "<<((TChain*) this->At(currChain))->GetNtrees()<<" file(s) with "<<((TChain*) this->At(currChain))->GetEntries()<<" events."<<endl;
        else
          return -1;

    }
  } else
    return -1;

  for (int i=0; i<((TChain*) this->At(currChain))->GetListOfBranches()->GetEntries(); ++i) {

    string s(((TChain*) this->At(currChain))->GetListOfBranches()->At(i)->GetName());

    string branch_alias = s;
    string branch_name = s;
    if (s.find(".", s.size()-1) != string::npos)
      branch_name += "obj";

    size_t a = s.find("_");
    if (a != string::npos) {
      size_t b = s.find("_", a+1);
      if (b != string::npos) {
        size_t c = s.find("_", b+1);
        if (c != string::npos) {
          string _prod     =s.substr(0,a);
          string _label    =s.substr(a+1, b-a-1);
          string _instance =s.substr(b+1, c-b-1);
          branch_alias = _label;
          if(_instance.length() > 0 )  {
            branch_alias += "_";
            branch_alias += _instance;
            ((TChain*) this->At(currChain))->SetAlias(_instance.c_str(),
                                                      branch_name.c_str());
          }
          string branch_alias_full = _prod + "_" + branch_alias;
          ((TChain*) this->At(currChain))->SetAlias(branch_alias_full.c_str(),
                                                    branch_name.c_str());
        }
      }
    }

    ((TChain*) this->At(currChain))->SetAlias(branch_alias.c_str(),
                                              branch_name.c_str());

  }

  if( addEventInfo || addTrackJets || addTower || addTrackJets || addHitDetInfo || diTrackMass || addTrigger) {

    if(verbose)   cout<<"add friend tree with additional variables..."<<endl;

    //make file for tree friends (adding additional, computed branches)
    string friendFileName("/scratch/hh/current/cms/user/");
    friendFileName += gSystem->GetUserInfo()->fUser;
    friendFileName += "/temp/";
    friendFileName += fileLabel;
    friendFileName += ".root";
    string fileOpt = "update";
    if( recreateTree )    fileOpt = "recreate";

    TFile *f = new TFile(friendFileName.c_str(),fileOpt.c_str());

    if( f->IsZombie() ) return -1;

    string friendTreeName("friendTree");
    //char number = currChain;
    //friendTreeName += number;
    //TTree *friendTree = new TTree("friendTree","friendTree");

    if(! f->FindKey(friendTreeName.c_str())) {
      if(verbose)   cout<<"calculating additional variables..."<<endl;


      TTree *friendTree = new TTree(friendTreeName.c_str(),friendTreeName.c_str());
      fwlite::ChainEvent ev(fileNames);
      int _event, _run, _lumi;

      //tower data
      const int kMAX = 10000;
      int NobjTowCal;
      towet  = new float [ kMAX ];
      toweta = new float [ kMAX ];
      towphi = new float [ kMAX ];
      towen  = new float [ kMAX ];
      towem  = new float [ kMAX ];
      towhd  = new float [ kMAX ];
      towoe  = new float [ kMAX ];
      towid_phi = new int[ kMAX ];
      towid_eta = new int[ kMAX ];
      towid     = new int[ kMAX ];
      
      int nJetsKT;
      TrackJetKT = new float [100];
      int nTracks;
      TrackDetE = new float [3000];
      TrackDetEECAL = new float [3000];
      float DiTrackMass;
      bool Trigger;
      bool Trigger1;
      
      
      if( addEventInfo ) {
 

        friendTree->Branch("event", &_event, "event/I");
        friendTree->Branch("run", &_run, "run/I");
        friendTree->Branch("lumi", &_lumi, "lumi/I");

      }

      if( addTrackJets ) {


        friendTree->Branch("nJetsKT", &nJetsKT, "nJetsKT/I");
        friendTree->Branch("TrackJetKT", TrackJetKT, "TrackJetKT[nJetsKT]/F");
      }

      if( addTower) {

        // CaloTower branches
        friendTree->Branch( "NobjTowCal",&NobjTowCal,"NobjTowCal/I"            );
        friendTree->Branch( "TowId",     towid,      "TowId[NobjTowCal]/I"     );
        friendTree->Branch( "TowId_phi", towid_phi,  "TowId_phi[NobjTowCal]/I" );
        friendTree->Branch( "TowId_eta", towid_eta,  "TowId_eta[NobjTowCal]/I" );
        friendTree->Branch( "TowEt",     towet,      "TowEt[NobjTowCal]/F"     );
        friendTree->Branch( "TowEta",    toweta,     "TowEta[NobjTowCal]/F"    );
        friendTree->Branch( "TowPhi",    towphi,     "TowPhi[NobjTowCal]/F"    );
        friendTree->Branch( "TowE",      towen,      "TowE[NobjTowCal]/F"      );
        friendTree->Branch( "TowEm",     towem,      "TowEm[NobjTowCal]/F"     );
        friendTree->Branch( "TowHad",    towhd,      "TowHad[NobjTowCal]/F"    );
        friendTree ->Branch( "TowOE",     towoe,      "TowOE[NobjTowCal]/F"     );


      }

      if( addTrackJets){            
        friendTree->Branch("nJetsKT", &nJetsKT, "nJetsKT/I");
        friendTree->Branch("TrackJetKT", TrackJetKT, "TrackJetKT[nJetsKT]/F");
      }
      
      if( addHitDetInfo){    
        friendTree->Branch("nTracks", &nTracks, "nTracks/I");
        friendTree->Branch("TrackDetE", TrackDetE, "TrackDetE[nTracks]/F");
        friendTree->Branch("TrackDetEECAL", TrackDetEECAL, "TrackDetEECAL[nTracks]/F");
      }
      if(diTrackMass)  friendTree->Branch("DiTrackMass", &DiTrackMass, "DiTrackMass/F");
      if(addTrigger)  {
        friendTree->Branch("Trigger", &Trigger, "Trigger/O");
        friendTree->Branch("Trigger1", &Trigger1, "Trigger1/O");
      }

      int tenth = ev.size() / 10;
      int eventNum =0;
      for (ev.toBegin(); !ev.atEnd(); ++ev, ++eventNum) {


        if( addEventInfo ) {


          if(eventNum==0) {
            //    fwlite::Handle<edm::TriggerResults> hltHandle;
            //    hltHandle.getByLabel(ev, "TriggerResults");

            //    fwlite::Handle<edm::TriggerResults> hTriggerResults;

            //    hTriggerResults.getByLabel(ev, "TriggerResults","","TEST");
            //    fwlite::TriggerNames const&  triggerNames = ev.triggerNames(*hTriggerResults);


            //    //      std::vector<std::string> const& names = triggerNames.triggerNames();
            //    for (unsigned i = 0; i < triggerNames.size(); ++i) {
            //    std::cout << i << "  " << triggerNames.triggerName(i) << std::endl;
            //    }
          }

          _event = ev.id().event();
          _run   = ev.id().run();
          _lumi = ev.luminosityBlock();

        }

        if( addTrackJets)
          {
            fwlite::Handle<reco::CaloJetCollection> jets;
            jets.getByLabel(ev, "antikt5CaloJets");
            
            fwlite::Handle<reco::TrackCollection> tracks;
            tracks.getByLabel(ev, "generalTracks");
            
            if (!jets.isValid())
              continue;
            if (!tracks.isValid())
              continue;
            double trackSum;
            nJetsKT = 0;
            for (reco::CaloJetCollection::const_iterator jet = jets->begin(); jet
                   !=jets->end(); jet++) {
              trackSum = 0;
              for (reco::TrackCollection::const_iterator track = tracks->begin(); track
                     !=tracks->end(); track++) {
                if (deltaR(jet->eta(), jet->phi(), track->eta(), track->phi())
                    < 0.5)
                  trackSum += track->pt();
              }
              TrackJetKT[nJetsKT] = trackSum;
              nJetsKT++;
            }
          }
        
        if(addHitDetInfo)
          {  
            fwlite::Handle<reco::TrackCollection> tracks;
            tracks.getByLabel(ev, "generalTracks");
            
            nTracks = 0;
            
            if (!tracks.isValid())
              continue;
            
            for (reco::TrackCollection::const_iterator track = tracks->begin(); track
                   !=tracks->end(); track++) {
              fwlite::Handle<vector <reco::CaloCluster> > CaloClusters;
              CaloClusters.getByLabel(ev, "hybridSuperClusters");  //add instance
              //CaloClusters.getByLabel(ev, "multi5x5BasicClusters");
              
              //fwlite::Handle<vector <CaloTower> > CaloTowers;
              //CaloTowers.getByLabel(ev, "towerMaker");
              
              double DR = 100;
              double DetE = 0;
              double DetEECAL = 0;
              
              if (CaloClusters.isValid()){   //always fails
                //if (CaloTowers.isValid()){
                for (vector<reco::CaloCluster>::const_iterator CaloCluster = CaloClusters->begin(); CaloCluster
                       !=CaloClusters->end(); CaloCluster++) {
                  //   int DetectorId = track->outerDetId();      //unfortuately only the Tracker DetId...
                  // /DetId *dId = new DetId(DetectorId);
                  //cout<<dId->det()<<endl;
                  
                  
                  if(deltaR( CaloCluster->eta(), CaloCluster->phi(), track->outerEta(), track->phi()) < DR)
                    {
                      DetE = CaloCluster->energy();
                      if(CaloCluster->caloID().detector(reco::CaloID::DET_ECAL_BARREL) ||CaloCluster->caloID().detector(reco::CaloID::DET_ECAL_ENDCAP))
                        DetEECAL = CaloCluster->energy();
                    }
                }
              }
              
              TrackDetE[nTracks] = DetE;
              TrackDetEECAL[nTracks] = DetEECAL;
              
              
              nTracks++;
              if(nTracks > 2990) cout<<"To many Tracks!!!!!!!!!"<<endl;
            }
          }
        
        if(diTrackMass)
          {
            fwlite::Handle<reco::TrackCollection> tracks;
            tracks.getByLabel(ev, "generalTracks");
            
            DiTrackMass = 0;
            
            if (!tracks.isValid())
              continue;
            
            for (reco::TrackCollection::const_iterator track = tracks->begin(); track
                   !=tracks->end(); track++) {
              if(!(track->d0() < 2 && track->dz() < 10 && track->numberOfValidHits() > 7)) continue;
              for (reco::TrackCollection::const_iterator track1 = tracks->begin(); track1
                     !=tracks->end(); track1++) {
                if(!(track1->d0() < 2 && track1->dz() < 10 && track1->numberOfValidHits() > 7) || track1 == track) continue;
                
                TLorentzVector tr1(track->px(), track->py(), track->pz(), track->p() );
                TLorentzVector tr2(track1->px(), track1->py(),  track1->pz(),  track1->p() );
                
                double result = (tr1 + tr2).Mag();
                
                if(result > DiTrackMass) DiTrackMass = result;
                
              }
            }
            
          }
        
        if(addTrigger)
          {
            
            Trigger = false;
            Trigger1 = false;
            
            fwlite::Handle< L1GlobalTriggerReadoutRecord> gtReadoutRecord;
            gtReadoutRecord.getByLabel(ev, "gtDigis");  
            //edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
            //event.getByLabel( edm::InputTag("gtDigis"), gtReadoutRecord);
            
            if(!gtReadoutRecord.isValid()) continue;
            
            
            const TechnicalTriggerWord&  technicalTriggerWordBeforeMask  = gtReadoutRecord->technicalTriggerWord();
            
            
            //std::vector<bool>& technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
            
            bool techTrigger0 = technicalTriggerWordBeforeMask.at(0);
            bool techTrigger40 = technicalTriggerWordBeforeMask.at(40);
            bool techTrigger41 = technicalTriggerWordBeforeMask.at(41);
            bool techTrigger36 = technicalTriggerWordBeforeMask.at(36);
            bool techTrigger37 = technicalTriggerWordBeforeMask.at(37);
            bool techTrigger38 = technicalTriggerWordBeforeMask.at(38);
            bool techTrigger39 = technicalTriggerWordBeforeMask.at(39);
            
            if(techTrigger0 && (techTrigger40 || techTrigger41)) Trigger = true;
            if(Trigger && !techTrigger36 && !techTrigger37 && !techTrigger38 && !techTrigger39) Trigger1 = true;
          }

        if( addTower) {
          fwlite::Handle<CaloTowerCollection> towers;
          towers.getByLabel(ev, "towerMaker");
          
          int jtow = 0;
          NobjTowCal=towers->size();
          for(CaloTowerCollection::const_iterator tow = towers->begin();
              tow != towers->end(); ++tow, ++jtow){
            towet [jtow] = tow->et();
            toweta[jtow] = tow->eta();
            towphi[jtow] = tow->phi();
            towen [jtow] = tow->energy();
            towem [jtow] = tow->emEnergy();
            towhd [jtow] = tow->hadEnergy();
            towoe [jtow] = tow->outerEnergy();
            towid_phi[jtow] = tow->id().iphi();
            towid_eta[jtow] = tow->id().ieta();
            towid [jtow] = tow->id().rawId();
          }
        }
        friendTree->Fill();
        
        if( eventNum != 0 && eventNum%tenth == 0) cout<<"Processed "<<eventNum <<" of "<<ev.size()<<" events. "<<endl;
        
      }
      f->cd();
      friendTree->Write();
    }
    f->Close();
    ((TChain*) this->At(currChain))->AddFriend(friendTreeName.c_str(),
                                               friendFileName.c_str());
    
  }


  timer.Stop();
  if(verbose) timer.Print();
       
  return this->GetEntries();

}
//------------------------------------------------------------------
//Draw given Chain with given variable and given cuts
int PlotTool::plot(int chainIndex, string histName, string cutName,
                int nEntries, string drwOpt) {

        if (chainIndex >= this->GetEntries() )
                return -1;

        TStopwatch timer;
        if(verbose) {
          cout<<"Plot: "<<histName<<" "<<cutName<<" from chain "<<this->At(chainIndex)->GetName()<<endl;
          timer.Start();
        }

        int currN = nEntries;
        if (nEntries < 0)
                currN = ((TChain*) this->At(chainIndex))->GetEntries(); //nEntries<0 : all entries are plotted!

        //++++ Create and name Canvases according to global variables +++++++++++++
        ostringstream currHistName;
        if (samePad_trees)
                currHistName<<((TChain*) this->At(chainIndex))->GetName()<<":";
        if (samePad_vars)
                currHistName<<histName;
        if (samePad_cuts)
                currHistName<<"{"<<cutName<<"}";

        ostringstream currPadName;
        if (!samePad_trees)
                currPadName<<((TChain*) this->At(chainIndex))->GetName()<<":";
        if (!samePad_vars)
                currPadName<<histName;
        if (!samePad_cuts)
                currPadName<<"{"<<cutName<<"}";

        ostringstream currCanvName;
        if (!sameCanv_trees && !samePad_trees)
                currCanvName<<((TChain*) this->At(chainIndex))->GetName()<<":";
        if (!sameCanv_vars && !samePad_vars)
                currCanvName<<histName;
        if (!sameCanv_cuts && !samePad_cuts)
                currCanvName<<"{"<<cutName<<"}";

        if ( (sameCanv_trees || samePad_trees) && (sameCanv_vars || samePad_vars)
                        && (sameCanv_cuts || samePad_cuts))
                currCanvName<<"All";

        string currOpt = drwOpt;
        bool useSubPads = (sameCanv_trees && !samePad_trees) || (sameCanv_vars
                        && !samePad_vars) || (sameCanv_cuts && !samePad_cuts);

        if (useSubPads) {
                if ( !canvases_[currCanvName.str()] || !gROOT->GetListOfCanvases()->FindObject(currCanvName.str().c_str()) ) {
                        canvases_[currCanvName.str()] = new TCanvas(currCanvName.str().c_str(), currCanvName.str().c_str());
                } else if (canvases_[currCanvName.str()]->GetNumber() >= 0) {
                        canvases_[currCanvName.str()]->SetNumber(canvases_[currCanvName.str()]->GetNumber()+1);
                }
        }

        if ( !pads_[currPadName.str()] || !gROOT->GetListOfCanvases()->FindObject(currPadName.str().c_str()) ) {
                pads_[currPadName.str()] = new TCanvas(currPadName.str().c_str(), currPadName.str().c_str());
                //              if (logY)
                //                      pads_[currPadName.str()]->SetLogy(1);
        } else {
                pads_[currPadName.str()]->cd();
                currOpt += "SAMES";
                if (useSubPads)
                        canvases_[currCanvName.str()]->SetNumber(canvases_[currCanvName.str()]->GetNumber()-1);
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


        //Draw histogram:
        string cutNameWithGlobal = cutName;
        if(cutName!="" && globalCuts!="")    cutNameWithGlobal += "&&";
        if(globalCuts!="")                   cutNameWithGlobal += globalCuts;

        int draw_err = ((TChain*) this->At(chainIndex))->Draw(histName.c_str(), cutNameWithGlobal.c_str(),
                        currOpt.c_str(), currN);
        if (draw_err < 0)
                return draw_err;

        //++++ Fix for histos with no entries +++++++++++++
        if (draw_err == 0) {
                if (currOpt.find("SAMES") == string::npos) {
                        pads_[currPadName.str()]->Close();
                        pads_.erase(currPadName.str() );
                        if (useSubPads)
                                canvases_[currCanvName.str()]->SetNumber(canvases_[currCanvName.str()]->GetNumber()-1);
                }
                cout<< "Warning: "<<currHistName.str().c_str()<<" in "<<currPadName.str()
                                <<" has no entries and is not drawn!"<<endl;
                return draw_err;
        }
        //++++++++++++++++++++++++++++++++++++++++++++++++++

        ((TH1F*) pads_[currPadName.str()]->GetPrimitive("htemp"))->SetName(currHistName.str().c_str()); // Set Name of histogram

        // Update for 2D histograms
        if( pads_[currPadName.str()]->GetPrimitive(currHistName.str().c_str())->InheritsFrom("TH2") ) {
                        pads_[currPadName.str()]->Draw();
                        setCanvas( pads_[currPadName.str()] );
        }


        if(verbose && draw_err > 0) {
          timer.Stop();
          cout<<"Done: Selected "<<draw_err<<" objects in "<< currN <<" processed events."<<endl;
          timer.Print();
        }
        return draw_err;

}
//------------------------------------------------------------------
//Standard loop to draw all chains and multiple variables and cuts
int PlotTool::loop(vector<string> _histName, vector<string> _cutName,
                int nEntries, string drwOpt, bool correspond) {

        int numProcessed = 0;

        if (correspond == false) {

                for (int i=0; i<this->GetEntries(); ++i) {
                        for (vector<string>::iterator it1 =_histName.begin(); it1
                                        != _histName.end(); ++it1) {
                                for (vector<string>::iterator it2 =_cutName.begin(); it2
                                                != _cutName.end(); ++it2) {

                                        numProcessed += plot(i, *it1, *it2, nEntries, drwOpt);
                                

                                }
                        }
                }
        } else {

                if (_histName.size() != _cutName.size() )
                        return -1;

                for (int i=0; i<this->GetEntries(); ++i) {
                        for (vector<string>::iterator it1 =_histName.begin(), it2 =
                                        _cutName.begin(); it1 != _cutName.end()&&it2
                                        != _cutName.end(); ++it1, ++it2) {

                                numProcessed += plot(i, *it1, *it2, nEntries, drwOpt);

                        }
                }
        }

        updatePads();
        fillCanvases();

        return numProcessed;
}

//Helper methods to allow for using simple strings as input
int PlotTool::loop(string histName, string cutName, int nEntries, string drwOpt) {

        vector<string> _histName;
        _histName.push_back(histName);
        vector<string> _cutName;
        _cutName.push_back(cutName);

        return loop(_histName, _cutName, nEntries, drwOpt);

}

int PlotTool::loop(vector<string> _histName, string cutName, int nEntries,
                string drwOpt) {

        vector<string> _cutName;
        _cutName.push_back(cutName);

        return loop(_histName, _cutName, nEntries, drwOpt);

}

int PlotTool::loop(string histName, vector<string> _cutName, int nEntries,
                string drwOpt) {

        vector<string> _histName;
        _histName.push_back(histName);

        return loop(_histName, _cutName, nEntries, drwOpt);
}
//------------------------------------------------------------------
//redraw the canvas to make changes in style visible
int PlotTool::updatePads() {

        for (map< string, TCanvas* >::iterator it=pads_.begin() ; it != pads_.end(); ++it) {
                if ( ! gROOT->GetListOfCanvases()->FindObject((*it).first.c_str() ) ) {
                  pads_.erase(it);
                  continue;
                }
                if( (*it).second->GetListOfPrimitives()->GetEntries() == 0 ) {
                  (*it).second->Close();
                  pads_.erase(it);
                  continue;
                }
                (*it).second->Draw();
                setCanvas((*it).second);

        }

        return pads_.size();

}

int PlotTool::fillCanvases() {

        for (map< string, TCanvas* >::iterator it=canvases_.begin() ; it
                        != canvases_.end(); ++it) {
                string canvName = (*it).first;
                if (gROOT->GetListOfCanvases()->FindObject(canvName.c_str()) ) {

                        int numP = (*it).second->GetNumber()+1;

                        if (numP <= 0)
                                continue;

                        int x = int( sqrt(numP) );
                        int y = x;
                        if (x*y < numP)
                                x += 1;
                        if (x*y < numP)
                                x += 1;
                        if (x*y < numP) {
                                x -= 1;
                                y += 1;
                        }


                        (*it).second->Divide(x, y);
                        int padIndex = 1;
                        for (map< string, TCanvas* >::iterator it2=pads_.begin() ; it2
                                        != pads_.end(); ++it2) {
                                string padName = (*it2).first;
                                if (gROOT->GetListOfCanvases()->FindObject(padName.c_str() ) ) {
                string n1 = canvName.substr(0,canvName.find(":"));
                string n2("");
                if(canvName.find("{") != string::npos ) n2 = canvName.substr(canvName.find("{"),canvName.find("}"));

                if ( (padName.find(canvName.c_str()) != string::npos ) || (padName.find(n1) != string::npos && padName.find(n2) != string::npos) || !strcmp( canvName.c_str(), "All") ) {
                                                (*it).second->cd(padIndex);
                                                (*it2).second->DrawClonePad();
                                                (*it2).second->Close();
                                                pads_.erase(it2);
                                                ++padIndex;
                                        }
                                } else
                                        pads_.erase(it2);
                        }
                        (*it).second->SetNumber(-1);

                } else
                        canvases_.erase(it);
        }

        return canvases_.size();

}
//------------------------------------------------------------------
void PlotTool::setCanvas(TCanvas* thisCanvas) {

        TH1* thisHist = 0;
        TPaveStats* thisStatsBox = 0;
        TPaletteAxis* palette =0;
        TLegend* thisLeg = 0;
        int counter =0;
        double maxEntry=0;

        ((TFrame*) thisCanvas->GetFrame())->Delete();
        thisCanvas->GetFrame()->SetLineWidth(gStyle->GetLineWidth() );

        thisCanvas->SetLeftMargin(0.2);
        thisCanvas->SetRightMargin(0.06);
        thisCanvas->SetBottomMargin(0.2);
        thisCanvas->SetTopMargin(0.1);

        if (logY)
                thisCanvas->SetLogy(1);
        else
                thisCanvas->SetLogy(0);

        if (showLegend)
                thisLeg = new TLegend();

        for (int i = 0; i != thisCanvas->GetListOfPrimitives()->GetSize(); ++i) {

                if ( !thisCanvas->GetListOfPrimitives()->At(i)->InheritsFrom("TH1"))
                        continue;

                thisHist = ((TH1*) thisCanvas->GetListOfPrimitives()->At(i));
                setColor(thisHist, counter);
                setMathLabels(thisHist);

                if (thisHist->GetMaximum() > maxEntry)
                        maxEntry = thisHist->GetMaximum();

                thisStatsBox = (TPaveStats*) thisHist->GetListOfFunctions()->FindObject("stats");
                if (thisStatsBox)
                        setStats(thisCanvas, thisStatsBox, thisHist, counter);

                palette = (TPaletteAxis*) thisHist->GetListOfFunctions()->FindObject("palette");
                if (palette)
                        setPalette(thisCanvas, palette);

                if (thisLeg)
                        thisLeg->AddEntry(thisHist,thisHist->GetName())->SetTextSize(0.04);

                ++counter;

        }

        if (maxEntry != 0)
                setHistMax(thisCanvas, maxEntry);

        thisCanvas->cd();
        thisCanvas->Update();

        if (thisLeg)
                setLegend(thisCanvas, thisLeg, counter);

}
//------------------------------------------------------------------
//private helper classes to set the canvas and hist style
void PlotTool::setStats(TCanvas* thisCanvas, TPaveStats* thisStatsBox,
                TH1* thisHist, int counter) {

        if (thisCanvas->GetRightMargin() < .2)
                thisCanvas->SetRightMargin(.2);

        thisStatsBox->SetLineColor(thisHist->GetLineColor());
        thisStatsBox->SetX2NDC(1.);
        thisStatsBox->SetY2NDC(1-thisCanvas->GetTopMargin()-0.16*counter);
        thisStatsBox->SetX1NDC(1-thisCanvas->GetRightMargin()+0.01);
        thisStatsBox->SetY1NDC(thisStatsBox->GetY2NDC()-.15);

}

void PlotTool::setMathLabels(TH1* thisHist) {

        string t = thisHist ->GetTitle();
        string x = thisHist->GetXaxis()->GetTitle();
        string y = thisHist->GetYaxis()->GetTitle();

//      if (x.find("__") != string::npos)
//        x = x.substr(0,x.find("__"));
//      if (y.find("__") != string::npos)
//        y = y.substr(0,y.find("__"));

        if (x.find(".phi()") != string::npos)
                x.replace(x.find(".phi()"), 6, " #phi");
        if (y.find(".phi()") != string::npos)
                y.replace(y.find(".phi()"), 6, " #phi");

        if (x.find(".eta()") != string::npos)
                x.replace(x.find(".eta()"), 6, " #eta");
        if (y.find(".eta()") != string::npos)
                y.replace(y.find(".eta()"), 6, " #eta");

        if (x.find(".theta()") != string::npos)
                x.replace(x.find(".theta()"), 6, " #theta");
        if (y.find(".theta()") != string::npos)
                y.replace(y.find(".theta()"), 6, " #theta");

        if (x.find(".chi2()") != string::npos)
                x.replace(x.find(".chi2()"), 7, " #chi^{2}");
        if (y.find(".chi2()") != string::npos)
                y.replace(y.find(".chi2()"), 7, " #chi^{2}");

        if (x.find(".pt()") != string::npos)
                x.replace(x.find(".pt()"), 5, " p_{T} [GeV]");
        if (y.find(".pt()") != string::npos)
                y.replace(y.find(".pt()"), 5, " p_{T} [GeV]");

        if (x.find(".et()") != string::npos)
                x.replace(x.find(".et()"), 5, " E_{T} [GeV]");
        if (y.find(".et()") != string::npos)
                y.replace(y.find(".et()"), 5, " E_{T} [GeV]");

        //splitlines for many cuts
        string test1= "{" + globalCuts + "}";
        string test2= "&&" + globalCuts;

        if(t.find(test1) != string::npos) t.replace(t.find(test1),test1.length(),"");
        if(t.find(test2) != string::npos) t.replace(t.find(test2),test2.length(),"");

        if (t.find("{") != string::npos && t.find("#splitline") == string::npos) {
          t.replace(t.find_last_of("{"), 1, "}{");
          string t_old = t;
          t = "#splitline{";
          t +=  t_old;
        }

        thisHist ->SetTitle(t.c_str());
        thisHist->GetXaxis()->SetTitle(x.c_str());
        thisHist->GetYaxis()->SetTitle(y.c_str());

}

void PlotTool::setColor(TH1* thisHist, int counter) {

        if (counter == 0) {

                thisHist->SetLineColor(kRed+1);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kRed+1);
                //else 
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kRed+1);
        } else if (counter == 1) {
                thisHist->SetLineColor(kBlue+1);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kBlue+1);
                //else 
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kBlue+1);

        } else if (counter == 2) {
                thisHist->SetLineColor(kGreen+2);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kGreen+2);
                //else
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kGreen+2);

        } else if (counter == 3) {
                thisHist->SetLineColor(kMagenta+2);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kMagenta+2);
                //else
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kMagenta+2);

        } else if (counter == 4) {
                thisHist->SetLineColor(kCyan+2);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kCyan+2);
                //else
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kCyan+2);

        } else {
                thisHist->SetLineColor(kBlack);
                //if( gStyle->GetHistFillStyle() ) thisHist->SetFillColor  (kBlack);
                //else
                thisHist->SetFillStyle(0);
                thisHist->SetMarkerColor(kBlack);

        }

}

void PlotTool::setPalette(TCanvas* thisCanvas, TPaletteAxis* palette) {

        palette->SetLabelSize(0.045);
        if (thisCanvas->GetRightMargin() < .15)
                thisCanvas->SetRightMargin(.15);
        palette->SetX1NDC(1-thisCanvas->GetRightMargin()+0.01);
        palette->SetY1NDC(thisCanvas->GetBottomMargin());
        palette->SetX2NDC(palette->GetX1NDC()+0.05);
        palette->SetY2NDC(1-thisCanvas->GetTopMargin());

}

void PlotTool::setHistMax(TCanvas* thisCanvas, double maxEntry) {

        TH1* thisHist = 0;

        for (int i = 0; i != thisCanvas->GetListOfPrimitives()->GetSize(); ++i) {

                if ( !thisCanvas->GetListOfPrimitives()->At(i)->InheritsFrom("TH1") )
                        continue;
                thisHist = ((TH1*) thisCanvas->GetListOfPrimitives()->At(i));
                if (thisHist->GetMaximum() < maxEntry) {
                        double minEntry=thisHist->GetBinContent(thisHist->GetMinimumBin());
                        if (logY) {
                                int bin = (thisHist->GetMinimumBin()+1);
                                minEntry=thisHist->GetBinContent(bin);
                        }
                        thisHist->GetYaxis()->SetRangeUser(minEntry, maxEntry*1.2);
                }
                break;

        }

}

void PlotTool::setLegend(TCanvas* thisCanvas, TLegend* thisLeg, int counter) {

        thisLeg->SetFillStyle(0);
        thisLeg->SetX1NDC(1-thisCanvas->GetRightMargin()-0.5);
        thisLeg->SetY1NDC(1-thisCanvas->GetTopMargin()-0.05*counter);
        thisLeg->SetX2NDC(1-thisCanvas->GetRightMargin());
        thisLeg->SetY2NDC(1-thisCanvas->GetTopMargin());
        thisLeg->SetEntrySeparation(0.5);
        thisLeg->Draw("NDC");

}

void PlotTool::createColors() {
        const Int_t NRGBs = 5;
        const Int_t NCont = 255;

        Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
        Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
        Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
        Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
        TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
        gStyle->SetNumberContours(NCont);

}

//------------------------------------------------------------------

int PlotTool::saveCanvases(string name, string path) {

        TSystemDirectory d("", path.c_str());
        if (!d.GetListOfFiles())
                return -1;

        string namePs = path +"/" + name;
        string nameRt = path +"/" + name;

        if (name.find(".ps") == string::npos)
                namePs += ".ps";

        if (name.find(".root") == string::npos)
                nameRt += ".root";

        int savedCanvs =0;


        TIter next(gROOT->GetListOfCanvases());
        TCanvas* canv;
        while ((canv=(TCanvas*)next())) {
          string modName = (*canv).GetName();
          if(modName.find("{") != string::npos) modName = modName.substr(0,modName.find("{"));
          (*canv).SetTitle( modName.c_str() );

        }


        TPostScript ps(namePs.c_str(),112);
        TFile rt(nameRt.c_str(),"recreate");

        next.Reset();
        while ((canv=(TCanvas*)next())) {
          ps.NewPage();


          (*canv).Write();
          (*canv).Draw();
          (*canv).Update();


          ++savedCanvs;

        }

        ps.Close();
        rt.Close();


        if(verbose && savedCanvs) {
          cout<<"Saved file "<<rt.GetName()<<" with "<<savedCanvs<<" canvases."<<endl;
          cout<<"Saved file "<<ps.GetName()<<" with "<<savedCanvs<<" pages."<<endl;
}
        return savedCanvs;

}
//------------------------------------------------------------------

int PlotTool::clearCanvases()
{

  while(gROOT->GetListOfCanvases()->GetEntries()) ((TCanvas*) gROOT->GetListOfCanvases()->At(0))->Close();
  pads_.clear();
  autoVars.clear();
  return pads_.size() + gROOT->GetListOfCanvases()->GetEntries();

}
//------------------------------------------------------------------
int PlotTool::setVariables(string label)
{

        for (int i=0; i< this->GetEntries(); ++i) {
          cout<<"--------------------------------"<<endl;
          cout<<((TChain*) this->At(i))->GetName()<<endl;
          cout<<"------------"<<endl;

          TList* leaves = (TList*) ((TChain*) this->At(i))->GetListOfLeaves();
          for (int j=0; j< leaves->GetEntries(); ++j) {
            TString leafName ( leaves->At(j)->GetName() );
            if(! leafName.EndsWith(".") ) continue;
            if(! leafName.Contains(label.c_str() ) ) continue;

            TClass cl( ( (TLeafElement*) leaves->At(j) )->GetTypeName() );
            cout<<"++++++"<<endl;
            cout<< leafName.Data() <<endl;
            for(int k=0;k<cl.GetListOfAllPublicMethods()->GetEntries();++k) {
              string typeName ( ((TMethod*) cl.GetListOfAllPublicMethods()->At(k))->GetReturnTypeName() );
              string methName ( cl.GetListOfAllPublicMethods()->At(k)->GetName() );
              if( methName != "product") continue;
                cout<< typeName <<endl;
                string _type;
                TString testString( typeName.c_str() );
                if( testString.BeginsWith("vector<") ) {
                  _type = typeName.substr(typeName.find("<")+1,typeName.find(">")-typeName.find("<")-1);
                  string _varSize = leafName.Data();
                  _varSize += "obj@.size()";
                  autoVars.push_back( _varSize );
                  }
                else _type = typeName.substr(0,typeName.find("*"));
                  TClass _cl( _type.c_str() );
                  for(int l=0;l<_cl.GetListOfAllPublicMethods()->GetEntries();++l) {
                    string _typeName ( ((TMethod*) _cl.GetListOfAllPublicMethods()->At(l))->GetReturnTypeName() );
                    string _methName ( _cl.GetListOfAllPublicMethods()->At(l)->GetName() );
                    //              if(_typeName.find("void") != string::npos ) continue;
                    //              cout<<"   "<<_typeName<<"  "<<_methName<<endl;

                    cout<<" "<<_typeName<<"  "<<_methName<<endl;
                    if(_methName.find("operator")==string::npos&&(_typeName=="float"||_typeName=="double"||_typeName=="int"||_typeName=="unsigned int"||_typeName=="bool"||_typeName=="unsigned short"||_typeName=="unsigned long"||_typeName=="unsigned long long")) {

                      cout<<"-->   "<<_typeName<<"  "<<_methName<<endl;

                      string  _varName = leafName.Data();
                      _varName += "obj.";
                      _varName += _methName;
                      _varName += "()";
                      autoVars.push_back( _varName );
                    }
                  }
                  
              
            }

          }

        }

        return autoVars.size();

}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//Draw efficiencies
int  PlotTool::plotEff(int chainIndex, string histName, string cutName, int nEntries, double fitXmin, double fitXmax, string fitFormula)
{

  if( chainIndex >= this->GetEntries() ) return -1;

  int currN = nEntries;
  if(nEntries < 0) currN = ((TChain*) this->At(chainIndex))->GetEntries(); //nEntries<0 : all entries are plotted!

  //++++ Create and name Canvases according to global variables +++++++++++++
  ostringstream currHistName;
  if( samePad_trees) currHistName<<((TChain*) this->At(chainIndex))->GetName()<<":";
  if( samePad_vars)  currHistName<<histName;
  if( samePad_cuts)  currHistName<<"{"<<cutName<<"}";

  ostringstream currPadName;
  if(! samePad_trees) currPadName<<((TChain*) this->At(chainIndex))->GetName()<<":";
  if(! samePad_vars)  currPadName<<histName;
  if(! samePad_cuts)  currPadName<<"{"<<cutName<<"}";


  //Draw total histogram:
  if( ! pads_["total"] || ! gROOT->FindObject("total") ) {
    pads_["total"] = new TCanvas("total", "total");
  } else {
    pads_["total"]->cd();
  }
  ostringstream bins_total;
  bins_total<<histName;
  bins_total<<">>total(100,0,1000)";


  int draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_total.str().c_str(), "", "", currN);
  if( draw_err <= 0 ) return draw_err;


  TH1* total = ((TH1*) pads_["total"]->GetPrimitive("total"));

  ostringstream bins_bkg;
  bins_bkg<<histName;
  bins_bkg<<">>bkg(";
  bins_bkg<<total->GetNbinsX();
  bins_bkg<<",";
  bins_bkg<<total->GetXaxis()->GetXmin();
  bins_bkg<<",";
  bins_bkg<<total->GetXaxis()->GetXmax();
  bins_bkg<<")";
 
  //Draw bkg histogram: 
  if( ! pads_["bkg"] || ! gROOT->FindObject("bkg") ) {
    pads_["bkg"] = new TCanvas("bkg", currPadName.str().c_str());
  } else {
    pads_["bkg"]->cd();
  }

  draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_bkg.str().c_str(), cutName.c_str(),"" , currN);
  if( draw_err <= 0 ) return draw_err;


  TH1* bkg = ((TH1*) pads_["bkg"]->GetPrimitive("bkg"));

  //Draw pass histogram: 
  ostringstream bins_pass;
  bins_pass<<histName;
  bins_pass<<">>pass(";
  bins_pass<<total->GetNbinsX();
  bins_pass<<",";
  bins_pass<<total->GetXaxis()->GetXmin();
  bins_pass<<",";
  bins_pass<<total->GetXaxis()->GetXmax();
  bins_pass<<")";

  ostringstream cut_pass;
  cut_pass<<"!(";
  cut_pass<<cutName;
  cut_pass<<")";
  

  if( ! pads_["pass"] || ! gROOT->FindObject("pass") ) {
    pads_["pass"] = new TCanvas("pass", currPadName.str().c_str());
  } else {
    pads_["pass"]->cd();
  }

  draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_pass.str().c_str(),cut_pass.str().c_str() ,"" , currN);
  if( draw_err <= 0 ) return draw_err;


  TH1* pass = ((TH1*) pads_["pass"]->GetPrimitive("pass"));


  currPadName<<"Eff";
  //Draw Efficiency Graph:
  if( ! pads_["EffGraph"] || ! gROOT->FindObject("EffGraph") ) {
    pads_["EffGraph"] = new TCanvas("EffGraph", currPadName.str().c_str());
  } else {
    pads_["EffGraph"]->cd();
  }

  TGraphAsymmErrors* EffGraph = new TGraphAsymmErrors(bkg, total);
  EffGraph->SetName(currHistName.str().c_str());

  TF1* reverse = new TF2("reverse","1/y-1",total->GetXaxis()->GetXmin(),total->GetXaxis()->GetXmax());
  EffGraph->Apply(reverse);
  EffGraph->Draw("A*");

  TF1* fitfunc = new TF1("fitfunc",fitFormula.c_str(),fitXmin,fitXmax);
  EffGraph->Fit("fitfunc","R+");


  //Save fit function

  ostringstream savefuncName;
  savefuncName<<histName;
  savefuncName<<"_";

  if(cutName.find("<") != std::string::npos ) cutName.replace(cutName.find("<"),1,"_");
  if(cutName.find(">") != std::string::npos ) cutName.replace(cutName.find(">"),1,"_");
  savefuncName<<cutName;


  TFile file("ABCDFunctions.root","UPDATE");
  TF1* savefunc = new TF1(savefuncName.str().c_str(), fitfunc->GetExpFormula().Data(),0,10000);
  savefunc->SetParameters( fitfunc->GetParameters() );
  savefunc->SetParErrors ( fitfunc->GetParErrors() );

  //Show results
  if( ! pads_["abcd"] || ! gROOT->FindObject("abcd") ) {
    pads_["abcd"] = new TCanvas("abcd", "abcd");
  } else {
    pads_["abcd"]->cd();
  }

  bkg->Multiply(savefunc);
  bkg->Draw();

  //  total->Add(bkg,-1);
  pass->SetLineColor(kRed);
  pass->Draw("sames,e");
  pads_["abcd"]->SetLogy();


  savefunc->Write();
  //  file.Close();


  return draw_err;

}
//------------------------------------------------------------------
Revision:	1.15
Committed:	Wed Dec 16 16:33:28 2009 UTC (15 years, 4 months ago) by thomsen
Content type:	text/plain
Branch:	MAIN
Changes since 1.14:	+169 -36 lines
Log Message:	Trigger & TrackDetE