ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/tschum/FWlite_Analysis/PlotTool.cc
(Generate patch)

Comparing UserCode/tschum/FWlite_Analysis/PlotTool.cc (file contents):
Revision 1.7 by tschum, Fri Nov 27 17:39:46 2009 UTC vs.
Revision 1.17 by thomsen, Mon Feb 8 16:51:26 2010 UTC

# Line 15 | Line 15 | PlotTool::PlotTool() {
15          showLegend = false;
16          logY = true;
17          addTrackJets = true;
18 +        addHitDetInfo = true;
19 +        diTrackMass = true;
20 +        addTrigger = true;
21 +        addTrackJets = false;
22 +        addEventInfo = false;
23 +        addTower = false;
24 +        addJetContent = false;
25 +        verbose = true;
26 +
27 +        globalCuts="";
28 +
29 +        recreateTree = false; //force recreation of friend tree
30  
31   }
32  
# Line 22 | Line 34 | PlotTool::PlotTool() {
34   //Fill object PlotTool with Chains constructed from files from given source
35  
36   int PlotTool::init(string fileName, string dirPath, string treeName,
37 <                string fileLabel) {
38 <        this->New(this->GetEntries() );
39 <        int currChain = this->GetEntries() - 1;
40 <
41 <        if (currChain < 0)
42 <                return currChain;
43 <
44 <        fileNames.clear();
45 <
46 <        //check if alternative label is different from default(searchString)
47 <        if (fileLabel=="") {
48 <                fileLabel=fileName;
49 <        }
50 <        ((TChain*) this->At(currChain))->SetName(fileLabel.c_str());
51 <        TSystemDirectory dir("sourceDir", dirPath.c_str());
52 <        TList *files = dir.GetListOfFiles();
53 <        if (files) {
54 <                TIter next(files);
55 <                TSystemFile *file;
56 <                TString fname;
57 <                string filePath;
58 <
59 <                cout<<"Open"<<dirPath.c_str()<<" Search for .root files that contain: "
60 <                                <<fileName.c_str()<<endl;
61 <
62 <                while ((file=(TSystemFile*)next())) {
63 <                        fname = file->GetName();
64 <                        if (!fname.EndsWith(".root"))
65 <                                continue;
66 <                        if (!fname.Contains(fileName.c_str()))
67 <                                continue;
37 >                   string fileLabel) {
38 >  this->New(this->GetEntries() );
39 >  int currChain = this->GetEntries() - 1;
40 >
41 >  if (currChain < 0)
42 >    return currChain;
43 >
44 >  TStopwatch timer;
45 >  timer.Start();
46 >
47 >
48 >  fileNames.clear();
49 >
50 >  //check if alternative label is different from default(searchString)
51 >  if (fileLabel=="") {
52 >    fileLabel=fileName;
53 >  }
54 >  ((TChain*) this->At(currChain))->SetName(fileLabel.c_str());
55 >
56 >  TList *files = new TList();
57 >  TSystemFile* sysFile = 0;
58 >  if(fileName.find(".") != string::npos) {
59 >    ifstream f(fileName.c_str());
60 >    if( ! f.is_open() ) return -1;
61 >    string line;
62 >
63 >    while (!f.eof()) {
64 >      getline(f,line);
65 >      sysFile = new TSystemFile(line.c_str(),dirPath.c_str());
66 >      files->Add(sysFile);
67 >    }
68 >
69 >
70 >  } else {
71 >    TSystemDirectory dir("sourceDir", dirPath.c_str());
72 >    files = dir.GetListOfFiles();
73 >  }
74 >  if (files->GetEntries()>0) {
75 >    TIter next(files);
76 >    TSystemFile *file;
77 >    TString fname;
78 >    string filePath;
79 >
80 >    if(verbose && fileName.find(".") == string::npos )  cout<<"Open"<<dirPath.c_str()<<" Search for .root files that contain: "
81 >                                                            <<fileName.c_str()<<endl;
82 >    if(verbose && fileName.find(".") != string::npos )  cout<<"Open"<<dirPath.c_str()<<" Search lines with .root in: "
83 >                                                            <<fileName.c_str()<<endl;
84 >
85 >
86 >    while ((file=(TSystemFile*)next())) {
87 >      fname = file->GetName();
88 >      if (!fname.EndsWith(".root"))
89 >        continue;
90 >      if (!fname.Contains(fileName.c_str()) && fileName.find(".") == string::npos )
91 >        continue;
92 >
93 >      filePath = dirPath;
94 >      filePath += fname.Data();
95 >      //fwlite::ChainEvent to lop over events, jets, etc
96 >      fileNames.push_back(filePath.c_str());
97 >
98 >      if (((TChain*) this->At(currChain))->AddFile(filePath.c_str(), -1,
99 >                                                   treeName.c_str()) )
100 >        if(verbose)     cout<<"Chained "<<((TChain*) this->At(currChain))->GetNtrees()<<" file(s) with "<<((TChain*) this->At(currChain))->GetEntries()<<" events."<<endl;
101 >        else
102 >          return -1;
103  
104 <                        filePath = dirPath;
105 <                        filePath += fname.Data();
106 <                        //fwlite::ChainEvent to lop over events, jets, etc
107 <                        fileNames.push_back(filePath.c_str());
108 <
109 <                        if (((TChain*) this->At(currChain))->AddFile(filePath.c_str(), -1,
110 <                                        treeName.c_str()) )
111 <                                cout<<"Chained "<<((TChain*) this->At(currChain))->GetNtrees()<<" file(s) with "<<((TChain*) this->At(currChain))->GetEntries()<<" events."<<endl;
112 <                        else
113 <                                return -1;
104 >    }
105 >  } else
106 >    return -1;
107 >
108 >  for (int i=0; i<((TChain*) this->At(currChain))->GetListOfBranches()->GetEntries(); ++i) {
109 >
110 >    string s(((TChain*) this->At(currChain))->GetListOfBranches()->At(i)->GetName());
111 >
112 >    string branch_alias = s;
113 >    string branch_name = s;
114 >    if (s.find(".", s.size()-1) != string::npos)
115 >      branch_name += "obj";
116 >
117 >    size_t a = s.find("_");
118 >    if (a != string::npos) {
119 >      size_t b = s.find("_", a+1);
120 >      if (b != string::npos) {
121 >        size_t c = s.find("_", b+1);
122 >        if (c != string::npos) {
123 >          string _prod     =s.substr(0,a);
124 >          string _label    =s.substr(a+1, b-a-1);
125 >          string _instance =s.substr(b+1, c-b-1);
126 >          branch_alias = _label;
127 >          if(_instance.length() > 0 )  {
128 >            branch_alias += "_";
129 >            branch_alias += _instance;
130 >            ((TChain*) this->At(currChain))->SetAlias(_instance.c_str(),
131 >                                                      branch_name.c_str());
132 >          }
133 >          string branch_alias_full = _prod + "_" + branch_alias;
134 >          ((TChain*) this->At(currChain))->SetAlias(branch_alias_full.c_str(),
135 >                                                    branch_name.c_str());
136 >        }
137 >      }
138 >    }
139  
140 <                }
141 <        } else
70 <                return -1;
140 >    ((TChain*) this->At(currChain))->SetAlias(branch_alias.c_str(),
141 >                                              branch_name.c_str());
142  
143 <        for (int i=0; i<((TChain*) this->At(currChain))->GetListOfBranches()->GetEntries(); ++i) {
143 >  }
144  
145 <                string s(((TChain*) this->At(currChain))->GetListOfBranches()->At(i)->GetName());
145 >  if( addEventInfo || addTrackJets || addTower || addTrackJets || addHitDetInfo || diTrackMass || addTrigger || addJetContent) {
146 >
147 >    if(verbose)   cout<<"add friend tree with additional variables..."<<endl;
148 >    
149 >    //make file for tree friends (adding additional, computed branches)
150 >    string friendFileName("/scratch/hh/current/cms/user/");
151 >    friendFileName += gSystem->GetUserInfo()->fUser;
152 >    friendFileName += "/temp/";
153 >    friendFileName += fileLabel;
154 >    friendFileName += ".root";
155 >    string fileOpt = "update";
156 >    if( recreateTree )    fileOpt = "recreate";
157 >
158 >    TFile *f = new TFile(friendFileName.c_str(),fileOpt.c_str());
159 >    
160 >    if( f->IsZombie() ) return -1;
161 >    
162 >    string friendTreeName("friendTree");
163 >    //char number = currChain;
164 >    //friendTreeName += number;
165 >    //TTree *friendTree = new TTree("friendTree","friendTree");
166 >
167 >    if(! f->FindKey(friendTreeName.c_str())) {
168 >      if(verbose)   cout<<"calculating additional variables..."<<endl;
169 >
170 >      TTree *friendTree = new TTree(friendTreeName.c_str(),friendTreeName.c_str());
171 >      fwlite::ChainEvent ev(fileNames);
172 >      int _event, _run, _lumi;
173 >
174 >      //tower data
175 >      const int kMAX = 10000;
176 >      int NobjTowCal;
177 >      towet  = new float [ kMAX ];
178 >      toweta = new float [ kMAX ];
179 >      towphi = new float [ kMAX ];
180 >      towen  = new float [ kMAX ];
181 >      towem  = new float [ kMAX ];
182 >      towhd  = new float [ kMAX ];
183 >      towoe  = new float [ kMAX ];
184 >      towid_phi = new int[ kMAX ];
185 >      towid_eta = new int[ kMAX ];
186 >      towid     = new int[ kMAX ];
187 >      
188 >      int nJetsKT;
189 >      TrackJetKT = new float [100];
190 >      int nTracks;
191 >      TrackDetE = new float [3000];
192 >      TrackDetEECAL = new float [3000];
193 >      float DiTrackMass;
194 >      bool Trigger;
195 >      bool Trigger1;
196 >
197 >
198 >      int NTowLeadJet;
199 >      TowLeadJetE = new float[1000];
200 >      TowLeadJetDPhi = new float[1000];
201 >      TowLeadJetDEta = new float[1000];
202 >      TowLeadJetEMF = new float[1000];
203 >      int NTrackLeadJet;
204 >      TrackLeadJetE = new float[1000];
205 >      TrackLeadJetDPhi = new float[1000];
206 >      TrackLeadJetDEta = new float[1000];
207 >      
208 >      
209 >      if( addEventInfo ) {
210 >
211 >
212 >        friendTree->Branch("event", &_event, "event/I");
213 >        friendTree->Branch("run", &_run, "run/I");
214 >        friendTree->Branch("lumi", &_lumi, "lumi/I");
215 >
216 >      }
217 >
218 >      if( addTrackJets ) {
219 >
220 >
221 >        friendTree->Branch("nJetsKT", &nJetsKT, "nJetsKT/I");
222 >        friendTree->Branch("TrackJetKT", TrackJetKT, "TrackJetKT[nJetsKT]/F");
223 >      }
224 >
225 >      if( addTower) {
226 >
227 >        // CaloTower branches
228 >        friendTree->Branch( "NobjTowCal",&NobjTowCal,"NobjTowCal/I"            );
229 >        friendTree->Branch( "TowId",     towid,      "TowId[NobjTowCal]/I"     );
230 >        friendTree->Branch( "TowId_phi", towid_phi,  "TowId_phi[NobjTowCal]/I" );
231 >        friendTree->Branch( "TowId_eta", towid_eta,  "TowId_eta[NobjTowCal]/I" );
232 >        friendTree->Branch( "TowEt",     towet,      "TowEt[NobjTowCal]/F"     );
233 >        friendTree->Branch( "TowEta",    toweta,     "TowEta[NobjTowCal]/F"    );
234 >        friendTree->Branch( "TowPhi",    towphi,     "TowPhi[NobjTowCal]/F"    );
235 >        friendTree->Branch( "TowE",      towen,      "TowE[NobjTowCal]/F"      );
236 >        friendTree->Branch( "TowEm",     towem,      "TowEm[NobjTowCal]/F"     );
237 >        friendTree->Branch( "TowHad",    towhd,      "TowHad[NobjTowCal]/F"    );
238 >        friendTree ->Branch( "TowOE",     towoe,      "TowOE[NobjTowCal]/F"     );
239 >
240 >
241 >      }
242 >
243 >      if( addTrackJets){            
244 >        friendTree->Branch("nJetsKT", &nJetsKT, "nJetsKT/I");
245 >        friendTree->Branch("TrackJetKT", TrackJetKT, "TrackJetKT[nJetsKT]/F");
246 >      }
247 >      
248 >      if( addHitDetInfo){    
249 >        friendTree->Branch("nTracks", &nTracks, "nTracks/I");
250 >        friendTree->Branch("TrackDetE", TrackDetE, "TrackDetE[nTracks]/F");
251 >        friendTree->Branch("TrackDetEECAL", TrackDetEECAL, "TrackDetEECAL[nTracks]/F");
252 >      }
253 >      if(diTrackMass)  friendTree->Branch("DiTrackMass", &DiTrackMass, "DiTrackMass/F");
254 >      if(addTrigger)  {
255 >        friendTree->Branch("Trigger", &Trigger, "Trigger/O");
256 >        friendTree->Branch("Trigger1", &Trigger1, "Trigger1/O");
257 >      }
258 >
259 >      if(addJetContent)
260 >        {
261 >        friendTree->Branch( "NTowLeadJet",&NTowLeadJet,"NTowLeadJet/I"           );
262 >        friendTree->Branch( "TowLeadJetE",TowLeadJetE,"TowLeadJetE[NTowLeadJet]/F" );
263 >        friendTree->Branch( "TowLeadJetDPhi",TowLeadJetDPhi,"TowLeadJetDPhi[NTowLeadJet]/F" );
264 >        friendTree->Branch( "TowLeadJetDEta",TowLeadJetDEta,"TowLeadJetDEta[NTowLeadJet]/F" );
265 >        friendTree->Branch( "TowLeadJetEMF",TowLeadJetEMF,"TowLeadJetEMF[NTowLeadJet]/F"    );
266 >        friendTree->Branch( "NTrackLeadJet",&NTrackLeadJet,"NTrackLeadJet/I"           );
267 >        friendTree->Branch( "TrackLeadJetE",TrackLeadJetE,"TrackLeadJetE[NTrackLeadJet]/F"  );
268 >        friendTree->Branch( "TrackLeadJetDPhi",TrackLeadJetDPhi,"TrackLeadJetDPhi[NTrackLeadJet]/F"  );
269 >        friendTree->Branch( "TrackLeadJetDEta",TrackLeadJetDEta,"TrackLeadJetDEta[NTrackLeadJet]/F"   );
270  
271 <                string branch_alias = s;
77 <                string branch_name = s;
78 <                if (s.find(".", s.size()-1) != string::npos)
79 <                        branch_name += "obj";
80 <
81 <                size_t a = s.find("_");
82 <                if (a != string::npos) {
83 <                  size_t b = s.find("_", a+1);
84 <                  if (b != string::npos) {
85 <                    size_t c = s.find("_", b+1);
86 <                    if (c != string::npos) {
87 <                      string _label    =s.substr(a+1, b-a-1);
88 <                      string _instance =s.substr(b+1, c-b-1);
89 <                      branch_alias = _label;
90 <                      if(_instance.length() > 0 )  {
91 <                        branch_alias += "_";
92 <                        branch_alias += _instance;
93 <                      }
94 <                    }
95 <                  }
96 <                }
271 >        }
272  
273 <                ((TChain*) this->At(currChain))->SetAlias(branch_alias.c_str(),
274 <                                branch_name.c_str());
273 >      int tenth = ev.size() / 10;
274 >      int eventNum =0;
275  
276 <        }
276 >      for (ev.toBegin(); !ev.atEnd(); ++ev, ++eventNum) {
277  
278 <        // add branch with track  Jets
278 >        if( addEventInfo ) {
279  
105        if( addTrackJets ) {
280  
281 <          //make file for tree friends (adding additional, computed branches)
282 <          string friendFileName("/scratch/hh/lustre/cms/user/thomsen/temp/");
283 <          friendFileName += fileName;
110 <          friendFileName += ".root";
111 <          TFile *f = new TFile(friendFileName.c_str(),"recreate");
112 <          string friendTreeName("friendTree");
113 <          //char number = currChain;
114 <          //friendTreeName += number;
115 <          //TTree *friendTree = new TTree("friendTree","friendTree");
116 <          TTree *friendTree = new TTree(friendTreeName.c_str(),friendTreeName.c_str());
281 >          if(eventNum==0) {
282 >            //    fwlite::Handle<edm::TriggerResults> hltHandle;
283 >            //    hltHandle.getByLabel(ev, "TriggerResults");
284  
285 <          int nJetsKT;
119 <          TrackJetKT = new float [100];
120 <          fwlite::ChainEvent ev(fileNames);
285 >            //    fwlite::Handle<edm::TriggerResults> hTriggerResults;
286  
287 <          friendTree->Branch("nJetsKT", &nJetsKT, "nJetsKT/I");
288 <          friendTree->Branch("TrackJetKT", TrackJetKT, "TrackJetKT[nJetsKT]/F");
287 >            //    hTriggerResults.getByLabel(ev, "TriggerResults","","TEST");
288 >            //    fwlite::TriggerNames const&  triggerNames = ev.triggerNames(*hTriggerResults);
289  
290 <          cout<<"calculating additional variables..."<<endl;
291 <          for (ev.toBegin(); !ev.atEnd(); ++ev) {
290 >            //    //      std::vector<std::string> const& names = triggerNames.triggerNames();
291 >            //    for (unsigned i = 0; i < triggerNames.size(); ++i) {
292 >            //    std::cout << i << "  " << triggerNames.triggerName(i) << std::endl;
293 >            //    }
294 >          }
295 >          _event = ev.id().event();
296 >          _run   = ev.id().run();
297 >          _lumi = ev.luminosityBlock();
298 >
299 >        }
300  
301 +        if( addTrackJets)
302 +          {
303              fwlite::Handle<reco::CaloJetCollection> jets;
304 <            jets.getByLabel(ev, "antikt5CaloJets");
304 >            jets.getByLabel(ev, "ak5CaloJets");
305 >            
306 >            fwlite::Handle<reco::TrackCollection> tracks;
307 >            tracks.getByLabel(ev, "generalTracks");
308 >            
309  
310 +            if (jets.isValid() && tracks.isValid()){
311 +              double trackSum;
312 +              nJetsKT = 0;
313 +              for (reco::CaloJetCollection::const_iterator jet = jets->begin(); jet
314 +                     !=jets->end(); jet++) {
315 +                trackSum = 0;
316 +                for (reco::TrackCollection::const_iterator track = tracks->begin(); track
317 +                       !=tracks->end(); track++) {
318 +                  if (deltaR(jet->eta(), jet->phi(), track->eta(), track->phi())
319 +                      < 0.5)
320 +                    trackSum += track->pt();
321 +                }
322 +                TrackJetKT[nJetsKT] = trackSum;
323 +                nJetsKT++;
324 +              }
325 +            }
326 +          }
327 +        
328 +
329 +        if(addHitDetInfo)
330 +          {  
331 +            fwlite::Handle<reco::TrackCollection> tracks;
332 +            tracks.getByLabel(ev, "generalTracks");
333 +            
334 +            nTracks = 0;
335 +            
336 +              if (tracks.isValid()){
337 +                
338 +                for (reco::TrackCollection::const_iterator track = tracks->begin(); track
339 +                       !=tracks->end(); track++) {
340 +                  fwlite::Handle<vector <reco::CaloCluster> > CaloClusters;
341 +                  CaloClusters.getByLabel(ev, "hybridSuperClusters","hybridBarrelBasicClusters");//::hybridBarrelBasicClusters");  //add instance
342 +                  //CaloClusters.getByLabel(ev, "multi5x5BasicClusters");
343 +                  
344 +                  //fwlite::Handle<vector <CaloTower> > CaloTowers;
345 +                  //CaloTowers.getByLabel(ev, "towerMaker");
346 +                  
347 +                  double DR = 100;
348 +                  double DetE = 0;
349 +                  double DetEECAL = 0;    
350 +                  
351 +                  if (CaloClusters.isValid()){
352 +                    //if (CaloTowers.isValid()){
353 +                    for (vector<reco::CaloCluster>::const_iterator CaloCluster = CaloClusters->begin(); CaloCluster
354 +                           !=CaloClusters->end(); CaloCluster++) {
355 +                      //   int DetectorId = track->outerDetId();      //unfortuately only the Tracker DetId...
356 +                      // /DetId *dId = new DetId(DetectorId);
357 +                      //cout<<dId->det()<<endl;
358 +                      
359 +                      
360 +                      if(deltaR( CaloCluster->eta(), CaloCluster->phi(), track->outerEta(), track->phi()) < DR)
361 +                        {
362 +                          DR = deltaR( CaloCluster->eta(), CaloCluster->phi(), track->outerEta(), track->phi());
363 +                          DetE = CaloCluster->energy();
364 +                          if(CaloCluster->caloID().detector(reco::CaloID::DET_ECAL_BARREL) ||CaloCluster->caloID().detector(reco::CaloID::DET_ECAL_ENDCAP))
365 +                            DetEECAL = CaloCluster->energy();
366 +                        }
367 +                    }
368 +                  }
369 +                                
370 +                  TrackDetE[nTracks] = DetE;
371 +                  TrackDetEECAL[nTracks] = DetEECAL;
372 +                  
373 +                  nTracks++;
374 +                  if(nTracks > 2990) cout<<"To many Tracks!!!!!!!!!"<<endl;
375 +                }
376 +              }
377 +          }
378 +        
379 +        if(diTrackMass)
380 +          {
381              fwlite::Handle<reco::TrackCollection> tracks;
382              tracks.getByLabel(ev, "generalTracks");
383 +            
384 +            DiTrackMass = 0;
385 +            
386 +            if (tracks.isValid()){
387 +            
388 +              for (reco::TrackCollection::const_iterator track = tracks->begin(); track
389 +                     !=tracks->end(); track++) {
390 +                if(!(track->d0() < 2 && track->dz() < 10 && track->numberOfValidHits() > 7)) continue;
391 +                for (reco::TrackCollection::const_iterator track1 = tracks->begin(); track1
392 +                       !=tracks->end(); track1++) {
393 +                  if(!(track1->d0() < 2 && track1->dz() < 10 && track1->numberOfValidHits() > 7) || track1 == track) continue;
394 +                  
395 +                  TLorentzVector tr1(track->px(), track->py(), track->pz(), track->p() );
396 +                  TLorentzVector tr2(track1->px(), track1->py(),  track1->pz(),  track1->p() );
397 +                  
398 +                  double result = (tr1 + tr2).Mag();
399 +                  
400 +                  if(result > DiTrackMass) DiTrackMass = result;
401 +                  
402 +                }
403 +              }
404 +            }
405 +            
406 +          }
407 +        
408 +        if(addTrigger)
409 +          {
410 +            
411 +            Trigger = false;
412 +            Trigger1 = false;
413 +            
414 +            fwlite::Handle< L1GlobalTriggerReadoutRecord> gtReadoutRecord;
415 +            gtReadoutRecord.getByLabel(ev, "gtDigis");  
416 +            //edm::Handle< L1GlobalTriggerReadoutRecord > gtReadoutRecord;
417 +            //event.getByLabel( edm::InputTag("gtDigis"), gtReadoutRecord);
418 +            
419 +            if(gtReadoutRecord.isValid()){
420 +            
421 +            
422 +              const TechnicalTriggerWord&  technicalTriggerWordBeforeMask  = gtReadoutRecord->technicalTriggerWord();
423 +            
424 +            
425 +              //std::vector<bool>& technicalTriggerWordBeforeMask = gtReadoutRecord->technicalTriggerWord();
426 +              
427 +              bool techTrigger0 = technicalTriggerWordBeforeMask.at(0);
428 +              bool techTrigger40 = technicalTriggerWordBeforeMask.at(40);
429 +              bool techTrigger41 = technicalTriggerWordBeforeMask.at(41);
430 +              bool techTrigger36 = technicalTriggerWordBeforeMask.at(36);
431 +              bool techTrigger37 = technicalTriggerWordBeforeMask.at(37);
432 +              bool techTrigger38 = technicalTriggerWordBeforeMask.at(38);
433 +              bool techTrigger39 = technicalTriggerWordBeforeMask.at(39);
434 +              
435 +              if(techTrigger0 && (techTrigger40 || techTrigger41)) Trigger = true;
436 +              if(Trigger && !techTrigger36 && !techTrigger37 && !techTrigger38 && !techTrigger39) Trigger1 = true;
437 +            }
438 +          }
439 +
440 +        if( addTower) {
441 +          fwlite::Handle<CaloTowerCollection> towers;
442 +          towers.getByLabel(ev, "towerMaker");
443 +          
444 +          int jtow = 0;
445 +          NobjTowCal=towers->size();
446 +          for(CaloTowerCollection::const_iterator tow = towers->begin();
447 +              tow != towers->end(); ++tow, ++jtow){
448 +            towet [jtow] = tow->et();
449 +            toweta[jtow] = tow->eta();
450 +            towphi[jtow] = tow->phi();
451 +            towen [jtow] = tow->energy();
452 +            towem [jtow] = tow->emEnergy();
453 +            towhd [jtow] = tow->hadEnergy();
454 +            towoe [jtow] = tow->outerEnergy();
455 +            towid_phi[jtow] = tow->id().iphi();
456 +            towid_eta[jtow] = tow->id().ieta();
457 +            towid [jtow] = tow->id().rawId();
458 +          }
459 +        }
460 +        if(addJetContent) {
461 +          NTowLeadJet = 0;
462 +          NTrackLeadJet = 0;
463 +          fwlite::Handle<reco::CaloJetCollection> jets;
464 +          jets.getByLabel(ev, "ak5CaloJets");
465 +          
466 +          if(jets.isValid() && jets->size()){
467 +            double jetPhi =  jets->begin()->phi();
468 +            double jetEta =  jets->begin()->eta();
469 +            double jetE   =  jets->begin()->energy();
470 +            
471 +            fwlite::Handle<CaloTowerCollection> towers;
472 +            towers.getByLabel(ev, "towerMaker");
473 +            
474 +            if(towers.isValid()){
475 +              for(CaloTowerCollection::const_iterator tow = towers->begin();
476 +                  tow != towers->end(); ++tow){
477 +                if(deltaR(jetEta, jetPhi, tow->eta(), tow->phi()) > 0.8) continue;
478 +
479 +                TowLeadJetE [NTowLeadJet] = tow->energy()/jetE;
480 +                TowLeadJetDPhi [NTowLeadJet] = deltaPhi(jetPhi,tow->phi());
481 +                TowLeadJetDEta [NTowLeadJet] = abs(jetEta - tow->eta());
482 +                TowLeadJetEMF [NTowLeadJet] = tow->emEnergy()/tow->hadEnergy();
483 +                NTowLeadJet ++;
484 +              }
485 +            }
486  
487 <            if (!jets.isValid())
488 <              continue;
489 <            if (!tracks.isValid())
490 <              continue;
138 <            double trackSum;
139 <            nJetsKT = 0;
140 <            for (reco::CaloJetCollection::const_iterator jet = jets->begin(); jet
141 <                   !=jets->end(); jet++) {
142 <              trackSum = 0;
487 >            fwlite::Handle<reco::TrackCollection> tracks;
488 >            tracks.getByLabel(ev, "generalTracks");
489 >            if(tracks.isValid()){
490 >              
491                for (reco::TrackCollection::const_iterator track = tracks->begin(); track
492                       !=tracks->end(); track++) {
493 <                if (deltaR(jet->eta(), jet->phi(), track->eta(), track->phi())
494 <                    < 0.5)
495 <                  trackSum += track->pt();
493 >                if(!(track->d0() < 2 && track->dz() < 10 && track->numberOfValidHits() > 7)) continue;
494 >                if(deltaR(jetEta, jetPhi, track->eta(), track->phi()) > 0.8) continue;
495 >                
496 >                TrackLeadJetE[NTrackLeadJet] = track->p()/jetE;
497 >                TrackLeadJetDPhi[NTrackLeadJet] = deltaPhi(jetPhi,track->phi());
498 >                TrackLeadJetDEta[NTrackLeadJet] = abs(jetEta-track->eta());
499 >                
500 >                NTrackLeadJet ++;
501                }
149              TrackJetKT[nJetsKT] = trackSum;
150              nJetsKT++;
502              }
152            friendTree->Fill();
503            }
154          f->cd();
155          friendTree->Write();
156          f->Close();
157          ((TChain*) this->At(currChain))->AddFriend(friendTreeName.c_str(),
158                                                     friendFileName.c_str());
504          }
505 <        return this->GetEntries();
505 >        friendTree->Fill();
506 >        
507 >        if( eventNum != 0 && eventNum%tenth == 0) cout<<"Processed "<<eventNum <<" of "<<ev.size()<<" events. "<<endl;
508 >        
509 >      }
510 >      f->cd();
511 >      friendTree->Write();
512 >    }
513 >    f->Close();
514 >    ((TChain*) this->At(currChain))->AddFriend(friendTreeName.c_str(),
515 >                                               friendFileName.c_str());
516 >    
517 >  }
518 >
519 >
520 >  timer.Stop();
521 >  if(verbose) timer.Print();
522 >      
523 >  return this->GetEntries();
524  
525   }
526   //------------------------------------------------------------------
# Line 168 | Line 531 | int PlotTool::plot(int chainIndex, strin
531          if (chainIndex >= this->GetEntries() )
532                  return -1;
533  
534 +        TStopwatch timer;
535 +        if(verbose) {
536 +          cout<<"Plot: "<<histName<<" "<<cutName<<" from chain "<<this->At(chainIndex)->GetName()<<endl;
537 +          timer.Start();
538 +        }
539 +
540          int currN = nEntries;
541          if (nEntries < 0)
542                  currN = ((TChain*) this->At(chainIndex))->GetEntries(); //nEntries<0 : all entries are plotted!
# Line 215 | Line 584 | int PlotTool::plot(int chainIndex, strin
584  
585          if ( !pads_[currPadName.str()] || !gROOT->GetListOfCanvases()->FindObject(currPadName.str().c_str()) ) {
586                  pads_[currPadName.str()] = new TCanvas(currPadName.str().c_str(), currPadName.str().c_str());
587 <                if (logY)
588 <                        pads_[currPadName.str()]->SetLogy(1);
587 >                //              if (logY)
588 >                //                      pads_[currPadName.str()]->SetLogy(1);
589          } else {
590                  pads_[currPadName.str()]->cd();
591                  currOpt += "SAMES";
# Line 227 | Line 596 | int PlotTool::plot(int chainIndex, strin
596  
597  
598          //Draw histogram:
599 <        int draw_err = ((TChain*) this->At(chainIndex))->Draw(histName.c_str(), cutName.c_str(),
599 >        string cutNameWithGlobal = cutName;
600 >        if(cutName!="" && globalCuts!="")    cutNameWithGlobal += "&&";
601 >        if(globalCuts!="")                   cutNameWithGlobal += globalCuts;
602 >
603 >        int draw_err = ((TChain*) this->At(chainIndex))->Draw(histName.c_str(), cutNameWithGlobal.c_str(),
604                          currOpt.c_str(), currN);
605          if (draw_err < 0)
606                  return draw_err;
# Line 248 | Line 621 | int PlotTool::plot(int chainIndex, strin
621  
622          ((TH1F*) pads_[currPadName.str()]->GetPrimitive("htemp"))->SetName(currHistName.str().c_str()); // Set Name of histogram
623  
624 +        // Update for 2D histograms
625 +        if( pads_[currPadName.str()]->GetPrimitive(currHistName.str().c_str())->InheritsFrom("TH2") ) {
626 +                        pads_[currPadName.str()]->Draw();
627 +                        setCanvas( pads_[currPadName.str()] );
628 +        }
629 +
630  
631 +        if(verbose && draw_err > 0) {
632 +          timer.Stop();
633 +          cout<<"Done: Selected "<<draw_err<<" objects in "<< currN <<" processed events."<<endl;
634 +          timer.Print();
635 +        }
636          return draw_err;
637  
638   }
# Line 268 | Line 652 | int PlotTool::loop(vector<string> _histN
652                                                  != _cutName.end(); ++it2) {
653  
654                                          numProcessed += plot(i, *it1, *it2, nEntries, drwOpt);
655 <                                        updatePads();
655 >                                
656  
657                                  }
658                          }
# Line 289 | Line 673 | int PlotTool::loop(vector<string> _histN
673                  }
674          }
675  
676 <        //  updatePads();
676 >        updatePads();
677          fillCanvases();
678  
679          return numProcessed;
# Line 330 | Line 714 | int PlotTool::loop(string histName, vect
714   int PlotTool::updatePads() {
715  
716          for (map< string, TCanvas* >::iterator it=pads_.begin() ; it != pads_.end(); ++it) {
717 <                if (gROOT->GetListOfCanvases()->FindObject((*it).first.c_str() ) ) {
718 <                        (*it).second->Draw();
719 <                        setCanvas((*it).second);
720 <                } else
721 <                        pads_.erase(it);
717 >                if ( ! gROOT->GetListOfCanvases()->FindObject((*it).first.c_str() ) ) {
718 >                  pads_.erase(it);
719 >                  continue;
720 >                }
721 >                if( (*it).second->GetListOfPrimitives()->GetEntries() == 0 ) {
722 >                  (*it).second->Close();
723 >                  pads_.erase(it);
724 >                  continue;
725 >                }
726 >                (*it).second->Draw();
727 >                setCanvas((*it).second);
728 >
729          }
730  
731          return pads_.size();
# Line 345 | Line 736 | int PlotTool::fillCanvases() {
736  
737          for (map< string, TCanvas* >::iterator it=canvases_.begin() ; it
738                          != canvases_.end(); ++it) {
739 <                const char* canvName = (*it).first.c_str();
740 <                if (gROOT->GetListOfCanvases()->FindObject(canvName) ) {
739 >                string canvName = (*it).first;
740 >                if (gROOT->GetListOfCanvases()->FindObject(canvName.c_str()) ) {
741  
742                          int numP = (*it).second->GetNumber()+1;
743  
# Line 364 | Line 755 | int PlotTool::fillCanvases() {
755                                  y += 1;
756                          }
757  
758 +
759                          (*it).second->Divide(x, y);
760                          int padIndex = 1;
761                          for (map< string, TCanvas* >::iterator it2=pads_.begin() ; it2
762                                          != pads_.end(); ++it2) {
763                                  string padName = (*it2).first;
764                                  if (gROOT->GetListOfCanvases()->FindObject(padName.c_str() ) ) {
765 <                                        if ( (padName.find(canvName) != string::npos ) || !strcmp(
766 <                                                        canvName, "All") ) {
765 >                string n1 = canvName.substr(0,canvName.find(":"));
766 >                string n2("");
767 >                if(canvName.find("{") != string::npos ) n2 = canvName.substr(canvName.find("{"),canvName.find("}"));
768 >
769 >                if ( (padName.find(canvName.c_str()) != string::npos ) || (padName.find(n1) != string::npos && padName.find(n2) != string::npos) || !strcmp( canvName.c_str(), "All") ) {
770                                                  (*it).second->cd(padIndex);
771                                                  (*it2).second->DrawClonePad();
772                                                  (*it2).second->Close();
# Line 475 | Line 870 | void PlotTool::setMathLabels(TH1* thisHi
870          string x = thisHist->GetXaxis()->GetTitle();
871          string y = thisHist->GetYaxis()->GetTitle();
872  
873 <        if (t.find(".phi()") != string::npos)
874 <                t.replace(t.find(".phi()"), 6, " #phi");
873 > //      if (x.find("__") != string::npos)
874 > //        x = x.substr(0,x.find("__"));
875 > //      if (y.find("__") != string::npos)
876 > //        y = y.substr(0,y.find("__"));
877 >
878          if (x.find(".phi()") != string::npos)
879                  x.replace(x.find(".phi()"), 6, " #phi");
880          if (y.find(".phi()") != string::npos)
881                  y.replace(y.find(".phi()"), 6, " #phi");
882  
485        if (t.find(".eta()") != string::npos)
486                t.replace(t.find(".eta()"), 6, " #eta");
883          if (x.find(".eta()") != string::npos)
884                  x.replace(x.find(".eta()"), 6, " #eta");
885          if (y.find(".eta()") != string::npos)
886                  y.replace(y.find(".eta()"), 6, " #eta");
887  
888 <        if (t.find(".pt()") != string::npos)
889 <                t.replace(t.find(".pt()"), 5, " p_{T}");
888 >        if (x.find(".theta()") != string::npos)
889 >                x.replace(x.find(".theta()"), 6, " #theta");
890 >        if (y.find(".theta()") != string::npos)
891 >                y.replace(y.find(".theta()"), 6, " #theta");
892 >
893 >        if (x.find(".chi2()") != string::npos)
894 >                x.replace(x.find(".chi2()"), 7, " #chi^{2}");
895 >        if (y.find(".chi2()") != string::npos)
896 >                y.replace(y.find(".chi2()"), 7, " #chi^{2}");
897 >
898          if (x.find(".pt()") != string::npos)
899 <                x.replace(x.find(".pt()"), 5, " p_{T}");
899 >                x.replace(x.find(".pt()"), 5, " p_{T} [GeV]");
900          if (y.find(".pt()") != string::npos)
901 <                y.replace(y.find(".pt()"), 5, " p_{T}");
901 >                y.replace(y.find(".pt()"), 5, " p_{T} [GeV]");
902 >
903 >        if (x.find(".et()") != string::npos)
904 >                x.replace(x.find(".et()"), 5, " E_{T} [GeV]");
905 >        if (y.find(".et()") != string::npos)
906 >                y.replace(y.find(".et()"), 5, " E_{T} [GeV]");
907 >
908 >        //splitlines for many cuts
909 >        string test1= "{" + globalCuts + "}";
910 >        string test2= "&&" + globalCuts;
911 >
912 >        if(t.find(test1) != string::npos) t.replace(t.find(test1),test1.length(),"");
913 >        if(t.find(test2) != string::npos) t.replace(t.find(test2),test2.length(),"");
914 >
915 >        if (t.find("{") != string::npos && t.find("#splitline") == string::npos) {
916 >          t.replace(t.find_last_of("{"), 1, "}{");
917 >          string t_old = t;
918 >          t = "#splitline{";
919 >          t +=  t_old;
920 >        }
921  
922          thisHist ->SetTitle(t.c_str());
923          thisHist->GetXaxis()->SetTitle(x.c_str());
# Line 612 | Line 1035 | void PlotTool::createColors() {
1035  
1036   //------------------------------------------------------------------
1037  
1038 < int PlotTool::saveCanvases(string type, string path) {
1038 > int PlotTool::saveCanvases(string name, string path) {
1039  
1040          TSystemDirectory d("", path.c_str());
1041          if (!d.GetListOfFiles())
1042                  return -1;
620        if (type.find(".") == string::npos)
621                return -1;
1043  
1044 <        int savedFiles =0;
1044 >        string namePs = path +"/" + name;
1045 >        string nameRt = path +"/" + name;
1046 >
1047 >        if (name.find(".ps") == string::npos)
1048 >                namePs += ".ps";
1049 >
1050 >        if (name.find(".root") == string::npos)
1051 >                nameRt += ".root";
1052 >
1053 >        int savedCanvs =0;
1054 >
1055  
1056          TIter next(gROOT->GetListOfCanvases());
1057          TCanvas* canv;
1058          while ((canv=(TCanvas*)next())) {
1059 +          string modName = (*canv).GetName();
1060 +          if(modName.find("{") != string::npos) modName = modName.substr(0,modName.find("{"));
1061 +          (*canv).SetTitle( modName.c_str() );
1062 +
1063 +        }
1064 +
1065 +
1066 +        TPostScript ps(namePs.c_str(),112);
1067 +        TFile rt(nameRt.c_str(),"recreate");
1068 +
1069 +        next.Reset();
1070 +        while ((canv=(TCanvas*)next())) {
1071 +          ps.NewPage();
1072 +
1073  
1074 <                string s = "";
1075 <                s += path;
1076 <                s += canv->GetName();
632 <                s += type;
1074 >          (*canv).Write();
1075 >          (*canv).Draw();
1076 >          (*canv).Update();
1077  
1078 <                canv->SaveAs(s.c_str());
1079 <                ++savedFiles;
1078 >
1079 >          ++savedCanvs;
1080  
1081          }
1082  
1083 <        return savedFiles;
1083 >        ps.Close();
1084 >        rt.Close();
1085 >
1086  
1087 +        if(verbose && savedCanvs) {
1088 +          cout<<"Saved file "<<rt.GetName()<<" with "<<savedCanvs<<" canvases."<<endl;
1089 +          cout<<"Saved file "<<ps.GetName()<<" with "<<savedCanvs<<" pages."<<endl;
1090   }
1091 < //------------------------------------------------------------------
643 < void PlotTool::showChainInfo()
644 < {
1091 >        return savedCanvs;
1092  
1093 + }
1094 + //------------------------------------------------------------------
1095  
1096 <  cout<<endl;
1097 <  cout<<"****** Show Chain Information:"<<endl;
1096 > int PlotTool::clearCanvases()
1097 > {
1098  
1099 +  while(gROOT->GetListOfCanvases()->GetEntries()) ((TCanvas*) gROOT->GetListOfCanvases()->At(0))->Close();
1100 +  pads_.clear();
1101 +  autoVars.clear();
1102 +  return pads_.size() + gROOT->GetListOfCanvases()->GetEntries();
1103  
1104 <  this->ls();
1104 > }
1105 > //------------------------------------------------------------------
1106 > int PlotTool::setVariables(string label)
1107 > {
1108  
653  cout<<endl;
654  cout<<"We have "<<this->GetEntries()<<" TChains with following aliases set:"<<endl;
1109          for (int i=0; i< this->GetEntries(); ++i) {
1110 <          cout<<endl;
1110 >          cout<<"--------------------------------"<<endl;
1111            cout<<((TChain*) this->At(i))->GetName()<<endl;
1112 <          ((TChain*) this->At(i))->GetListOfAliases()->ls();
1112 >          cout<<"------------"<<endl;
1113 >
1114 >          TList* leaves = (TList*) ((TChain*) this->At(i))->GetListOfLeaves();
1115 >          for (int j=0; j< leaves->GetEntries(); ++j) {
1116 >            TString leafName ( leaves->At(j)->GetName() );
1117 >            if(! leafName.EndsWith(".") ) continue;
1118 >            if(! leafName.Contains(label.c_str() ) ) continue;
1119 >
1120 >            TClass cl( ( (TLeafElement*) leaves->At(j) )->GetTypeName() );
1121 >            cout<<"++++++"<<endl;
1122 >            cout<< leafName.Data() <<endl;
1123 >            for(int k=0;k<cl.GetListOfAllPublicMethods()->GetEntries();++k) {
1124 >              string typeName ( ((TMethod*) cl.GetListOfAllPublicMethods()->At(k))->GetReturnTypeName() );
1125 >              string methName ( cl.GetListOfAllPublicMethods()->At(k)->GetName() );
1126 >              if( methName != "product") continue;
1127 >                cout<< typeName <<endl;
1128 >                string _type;
1129 >                TString testString( typeName.c_str() );
1130 >                if( testString.BeginsWith("vector<") ) {
1131 >                  _type = typeName.substr(typeName.find("<")+1,typeName.find(">")-typeName.find("<")-1);
1132 >                  string _varSize = leafName.Data();
1133 >                  _varSize += "obj@.size()";
1134 >                  autoVars.push_back( _varSize );
1135 >                  }
1136 >                else _type = typeName.substr(0,typeName.find("*"));
1137 >                  TClass _cl( _type.c_str() );
1138 >                  for(int l=0;l<_cl.GetListOfAllPublicMethods()->GetEntries();++l) {
1139 >                    string _typeName ( ((TMethod*) _cl.GetListOfAllPublicMethods()->At(l))->GetReturnTypeName() );
1140 >                    string _methName ( _cl.GetListOfAllPublicMethods()->At(l)->GetName() );
1141 >                    //              if(_typeName.find("void") != string::npos ) continue;
1142 >                    //              cout<<"   "<<_typeName<<"  "<<_methName<<endl;
1143 >
1144 >                    cout<<" "<<_typeName<<"  "<<_methName<<endl;
1145 >                    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")) {
1146 >
1147 >                      cout<<"-->   "<<_typeName<<"  "<<_methName<<endl;
1148 >
1149 >                      string  _varName = leafName.Data();
1150 >                      _varName += "obj.";
1151 >                      _varName += _methName;
1152 >                      _varName += "()";
1153 >                      autoVars.push_back( _varName );
1154 >                    }
1155 >                  }
1156 >                  
1157 >              
1158 >            }
1159 >
1160 >          }
1161  
1162          }
1163  
1164 +        return autoVars.size();
1165 +
1166 + }
1167 + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1168 + //Draw efficiencies
1169 + int  PlotTool::plotEff(int chainIndex, string histName, string cutName, int nEntries, double fitXmin, double fitXmax, string fitFormula)
1170 + {
1171 +
1172 +  if( chainIndex >= this->GetEntries() ) return -1;
1173 +
1174 +  int currN = nEntries;
1175 +  if(nEntries < 0) currN = ((TChain*) this->At(chainIndex))->GetEntries(); //nEntries<0 : all entries are plotted!
1176 +
1177 +  //++++ Create and name Canvases according to global variables +++++++++++++
1178 +  ostringstream currHistName;
1179 +  if( samePad_trees) currHistName<<((TChain*) this->At(chainIndex))->GetName()<<":";
1180 +  if( samePad_vars)  currHistName<<histName;
1181 +  if( samePad_cuts)  currHistName<<"{"<<cutName<<"}";
1182 +
1183 +  ostringstream currPadName;
1184 +  if(! samePad_trees) currPadName<<((TChain*) this->At(chainIndex))->GetName()<<":";
1185 +  if(! samePad_vars)  currPadName<<histName;
1186 +  if(! samePad_cuts)  currPadName<<"{"<<cutName<<"}";
1187 +
1188 +
1189 +  //Draw total histogram:
1190 +  if( ! pads_["total"] || ! gROOT->FindObject("total") ) {
1191 +    pads_["total"] = new TCanvas("total", "total");
1192 +  } else {
1193 +    pads_["total"]->cd();
1194 +  }
1195 +  ostringstream bins_total;
1196 +  bins_total<<histName;
1197 +  bins_total<<">>total(100,0,1000)";
1198 +
1199 +
1200 +  int draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_total.str().c_str(), "", "", currN);
1201 +  if( draw_err <= 0 ) return draw_err;
1202 +
1203 +
1204 +  TH1* total = ((TH1*) pads_["total"]->GetPrimitive("total"));
1205 +
1206 +  ostringstream bins_bkg;
1207 +  bins_bkg<<histName;
1208 +  bins_bkg<<">>bkg(";
1209 +  bins_bkg<<total->GetNbinsX();
1210 +  bins_bkg<<",";
1211 +  bins_bkg<<total->GetXaxis()->GetXmin();
1212 +  bins_bkg<<",";
1213 +  bins_bkg<<total->GetXaxis()->GetXmax();
1214 +  bins_bkg<<")";
1215 +
1216 +  //Draw bkg histogram:
1217 +  if( ! pads_["bkg"] || ! gROOT->FindObject("bkg") ) {
1218 +    pads_["bkg"] = new TCanvas("bkg", currPadName.str().c_str());
1219 +  } else {
1220 +    pads_["bkg"]->cd();
1221 +  }
1222 +
1223 +  draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_bkg.str().c_str(), cutName.c_str(),"" , currN);
1224 +  if( draw_err <= 0 ) return draw_err;
1225 +
1226  
1227 +  TH1* bkg = ((TH1*) pads_["bkg"]->GetPrimitive("bkg"));
1228 +
1229 +  //Draw pass histogram:
1230 +  ostringstream bins_pass;
1231 +  bins_pass<<histName;
1232 +  bins_pass<<">>pass(";
1233 +  bins_pass<<total->GetNbinsX();
1234 +  bins_pass<<",";
1235 +  bins_pass<<total->GetXaxis()->GetXmin();
1236 +  bins_pass<<",";
1237 +  bins_pass<<total->GetXaxis()->GetXmax();
1238 +  bins_pass<<")";
1239 +
1240 +  ostringstream cut_pass;
1241 +  cut_pass<<"!(";
1242 +  cut_pass<<cutName;
1243 +  cut_pass<<")";
1244 +  
1245 +
1246 +  if( ! pads_["pass"] || ! gROOT->FindObject("pass") ) {
1247 +    pads_["pass"] = new TCanvas("pass", currPadName.str().c_str());
1248 +  } else {
1249 +    pads_["pass"]->cd();
1250 +  }
1251 +
1252 +  draw_err = ((TChain*) this->At(chainIndex))->Draw(bins_pass.str().c_str(),cut_pass.str().c_str() ,"" , currN);
1253 +  if( draw_err <= 0 ) return draw_err;
1254 +
1255 +
1256 +  TH1* pass = ((TH1*) pads_["pass"]->GetPrimitive("pass"));
1257 +
1258 +
1259 +  currPadName<<"Eff";
1260 +  //Draw Efficiency Graph:
1261 +  if( ! pads_["EffGraph"] || ! gROOT->FindObject("EffGraph") ) {
1262 +    pads_["EffGraph"] = new TCanvas("EffGraph", currPadName.str().c_str());
1263 +  } else {
1264 +    pads_["EffGraph"]->cd();
1265 +  }
1266 +
1267 +  TGraphAsymmErrors* EffGraph = new TGraphAsymmErrors(bkg, total);
1268 +  EffGraph->SetName(currHistName.str().c_str());
1269 +
1270 +  TF1* reverse = new TF2("reverse","1/y-1",total->GetXaxis()->GetXmin(),total->GetXaxis()->GetXmax());
1271 +  EffGraph->Apply(reverse);
1272 +  EffGraph->Draw("A*");
1273 +
1274 +  TF1* fitfunc = new TF1("fitfunc",fitFormula.c_str(),fitXmin,fitXmax);
1275 +  EffGraph->Fit("fitfunc","R+");
1276 +
1277 +
1278 +  //Save fit function
1279 +
1280 +  ostringstream savefuncName;
1281 +  savefuncName<<histName;
1282 +  savefuncName<<"_";
1283 +
1284 +  if(cutName.find("<") != std::string::npos ) cutName.replace(cutName.find("<"),1,"_");
1285 +  if(cutName.find(">") != std::string::npos ) cutName.replace(cutName.find(">"),1,"_");
1286 +  savefuncName<<cutName;
1287 +
1288 +
1289 +  TFile file("ABCDFunctions.root","UPDATE");
1290 +  TF1* savefunc = new TF1(savefuncName.str().c_str(), fitfunc->GetExpFormula().Data(),0,10000);
1291 +  savefunc->SetParameters( fitfunc->GetParameters() );
1292 +  savefunc->SetParErrors ( fitfunc->GetParErrors() );
1293 +
1294 +  //Show results
1295 +  if( ! pads_["abcd"] || ! gROOT->FindObject("abcd") ) {
1296 +    pads_["abcd"] = new TCanvas("abcd", "abcd");
1297 +  } else {
1298 +    pads_["abcd"]->cd();
1299 +  }
1300 +
1301 +  bkg->Multiply(savefunc);
1302 +  bkg->Draw();
1303 +
1304 +  //  total->Add(bkg,-1);
1305 +  pass->SetLineColor(kRed);
1306 +  pass->Draw("sames,e");
1307 +  pads_["abcd"]->SetLogy();
1308 +
1309 +
1310 +  savefunc->Write();
1311 +  //  file.Close();
1312 +
1313 +
1314 +  return draw_err;
1315 +
1316 + }
1317 + //------------------------------------------------------------------
1318 +
1319 +
1320 + //Eta-Phi to Det-komponent matching
1321 +
1322 + int PlotTool::getEtaBin(double eta){
1323 +  for(int i=-41;i<42;i++){
1324 +    double low = getEtaFromBin(i,true);
1325 +    double high = getEtaFromBin(i,false);
1326 +    if(eta >= low && eta <= high) return i;
1327 +  }
1328 +  cout<<"eta-bin not found"<<endl;
1329 +  return 0;
1330 + }
1331 +
1332 + double PlotTool::getEtaFromBin(int etaBin, bool lowerEdge){
1333 +  // return eta bin - eta edge mappting
1334 +  switch(etaBin){
1335 +  case -41: return (lowerEdge ? -5.191 : -4.889); break;
1336 +  case -40: return (lowerEdge ? -4.889 : -4.716); break;
1337 +  case -39: return (lowerEdge ? -4.716 : -4.538); break;
1338 +  case -38: return (lowerEdge ? -4.538 : -4.363); break;
1339 +  case -37: return (lowerEdge ? -4.363 : -4.191); break;
1340 +  case -36: return (lowerEdge ? -4.191 : -4.013); break;
1341 +  case -35: return (lowerEdge ? -4.013 : -3.839); break;
1342 +  case -34: return (lowerEdge ? -3.839 : -3.664); break;
1343 +  case -33: return (lowerEdge ? -3.664 : -3.489); break;
1344 +  case -32: return (lowerEdge ? -3.489 : -3.314); break;
1345 +  case -31: return (lowerEdge ? -3.314 : -3.139); break;
1346 +  case -30: return (lowerEdge ? -3.139 : -2.964); break;
1347 +  case -29: return (lowerEdge ? -2.964 : -2.853); break;
1348 +  case -28: return (lowerEdge ? -2.853 :  -2.65); break;
1349 +  case -27: return (lowerEdge ?  -2.65 :   -2.5); break;
1350 +  case -26: return (lowerEdge ?   -2.5 : -2.322); break;
1351 +  case -25: return (lowerEdge ? -2.322 : -2.172); break;
1352 +  case -24: return (lowerEdge ? -2.172 : -2.043); break;
1353 +  case -23: return (lowerEdge ? -2.043 :  -1.93); break;
1354 +  case -22: return (lowerEdge ?  -1.93 :  -1.83); break;
1355 +  case -21: return (lowerEdge ?  -1.83 :  -1.74); break;
1356 +  case -20: return (lowerEdge ?  -1.74 : -1.653); break;
1357 +  case -19: return (lowerEdge ? -1.653 : -1.566); break;
1358 +  case -18: return (lowerEdge ? -1.566 : -1.479); break;
1359 +  case -17: return (lowerEdge ? -1.479 : -1.392); break;
1360 +  case -16: return (lowerEdge ? -1.392 : -1.305); break;
1361 +  case -15: return (lowerEdge ? -1.305 : -1.218); break;
1362 +  case -14: return (lowerEdge ? -1.218 : -1.131); break;
1363 +  case -13: return (lowerEdge ? -1.131 : -1.044); break;
1364 +  case -12: return (lowerEdge ? -1.044 : -0.957); break;
1365 +  case -11: return (lowerEdge ? -0.957 : -0.879); break;
1366 +  case -10: return (lowerEdge ? -0.879 : -0.783); break;
1367 +  case  -9: return (lowerEdge ? -0.783 : -0.696); break;
1368 +  case  -8: return (lowerEdge ? -0.696 : -0.609); break;
1369 +  case  -7: return (lowerEdge ? -0.609 : -0.522); break;
1370 +  case  -6: return (lowerEdge ? -0.522 : -0.435); break;
1371 +  case  -5: return (lowerEdge ? -0.435 : -0.348); break;
1372 +  case  -4: return (lowerEdge ? -0.348 : -0.261); break;
1373 +  case  -3: return (lowerEdge ? -0.261 : -0.174); break;
1374 +  case  -2: return (lowerEdge ? -0.174 : -0.087); break;
1375 +  case  -1: return (lowerEdge ? -0.087 :      0); break;
1376 +  case  +1: return (lowerEdge ?      0 :  0.087); break;
1377 +  case  +2: return (lowerEdge ?  0.087 :  0.174); break;
1378 +  case  +3: return (lowerEdge ?  0.174 :  0.261); break;
1379 +  case  +4: return (lowerEdge ?  0.261 :  0.348); break;
1380 +  case  +5: return (lowerEdge ?  0.348 :  0.435); break;
1381 +  case  +6: return (lowerEdge ?  0.435 :  0.522); break;
1382 +  case  +7: return (lowerEdge ?  0.522 :  0.609); break;
1383 +  case  +8: return (lowerEdge ?  0.609 :  0.696); break;
1384 +  case  +9: return (lowerEdge ?  0.696 :  0.783); break;
1385 +  case +10: return (lowerEdge ?  0.783 :  0.879); break;
1386 +  case +11: return (lowerEdge ?  0.879 :  0.957); break;
1387 +  case +12: return (lowerEdge ?  0.957 :  1.044); break;
1388 +  case +13: return (lowerEdge ?  1.044 :  1.131); break;
1389 +  case +14: return (lowerEdge ?  1.131 :  1.218); break;
1390 +  case +15: return (lowerEdge ?  1.218 :  1.305); break;
1391 +  case +16: return (lowerEdge ?  1.305 :  1.392); break;
1392 +  case +17: return (lowerEdge ?  1.392 :  1.479); break;
1393 +  case +18: return (lowerEdge ?  1.479 :  1.566); break;
1394 +  case +19: return (lowerEdge ?  1.566 :  1.653); break;
1395 +  case +20: return (lowerEdge ?  1.653 :   1.74); break;
1396 +  case +21: return (lowerEdge ?   1.74 :   1.83); break;
1397 +  case +22: return (lowerEdge ?   1.83 :   1.93); break;
1398 +  case +23: return (lowerEdge ?   1.93 :  2.043); break;
1399 +  case +24: return (lowerEdge ?  2.043 :  2.172); break;
1400 +  case +25: return (lowerEdge ?  2.172 :  2.322); break;
1401 +  case +26: return (lowerEdge ?  2.322 :    2.5); break;
1402 +  case +27: return (lowerEdge ?    2.5 :   2.65); break;
1403 +  case +28: return (lowerEdge ?   2.65 :  2.853); break;
1404 +  case +29: return (lowerEdge ?  2.853 :  2.964); break;
1405 +  case +30: return (lowerEdge ?  2.964 :  3.139); break;
1406 +  case +31: return (lowerEdge ?  3.139 :  3.314); break;
1407 +  case +32: return (lowerEdge ?  3.314 :  3.489); break;
1408 +  case +33: return (lowerEdge ?  3.489 :  3.664); break;
1409 +  case +34: return (lowerEdge ?  3.664 :  3.839); break;
1410 +  case +35: return (lowerEdge ?  3.839 :  4.013); break;
1411 +  case +36: return (lowerEdge ?  4.013 :  4.191); break;
1412 +  case +37: return (lowerEdge ?  4.191 :  4.363); break;
1413 +  case +38: return (lowerEdge ?  4.363 :  4.538); break;
1414 +  case +39: return (lowerEdge ?  4.538 :  4.716); break;
1415 +  case +40: return (lowerEdge ?  4.716 :  4.889); break;
1416 +  case +41: return (lowerEdge ?  4.889 :  5.191); break;
1417 +    //something went wrong;
1418 +  default : return -1; break;
1419 +  }
1420 + }
1421 +
1422 + int PlotTool::getPhiBin(double phi){
1423 +
1424 +  return 0;
1425 + }
1426 +        
1427 + double PlotTool::getPhiFromBin(int phiBin, bool lowerEdge){
1428  
1429 +  return 0;
1430   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines