ROCsoft: beamtime/cern-oct11/go4/AnalysisMacros/electron_pion_separation.cxx (r4864/r3051) Source File

Go to the documentation of this file.
00001 #ifndef __CINT__
00002 #include "Riostream.h"
00003 #include "TRint.h"
00004 #include "TROOT.h"
00005 #include "TStyle.h"
00006 #include "Riostream.h"
00007 
00008 
00009 #include "TH2.h"
00010 #include "TH1.h"
00011 #include "TF1.h"
00012 
00013 #include "TGraph.h"
00014 #include "TMultiGraph.h"
00015 #include "TCanvas.h"
00016 #include "TProfile.h"
00017 #include "TStopwatch.h"
00018 
00019 #include "TFile.h"
00020 #include "TTree.h"
00021 #include "TBranch.h"
00022 #include "TMath.h"
00023 #include "TLegend.h"
00024 #include "TColor.h"
00025 #include "TText.h"
00026 #include "TKey.h"
00027 #include "TLatex.h"
00028 #include "TCanvas.h"
00029 #include "TLine.h"
00030 
00031 #endif
00032 
00033 using namespace std;
00034 
00035 void Statusbar(Int_t i, Int_t n) {
00036   if (int(i * 100 / float(n)) - int((i-1) * 100 / float(n)) >= 1) {
00037     if (int(i * 100 / float(n)) == 1 || i == 1 || i == 0) 
00038       cout << "[" << flush;
00039     cout << "-" << flush;
00040     if (int(i * 10 / float(n)) - int((i-1) * 10 / float(n)) >= 1) 
00041       cout << "|";
00042     if (int(i * 100 / float(n)) >=99) 
00043       cout << "]" <<endl;
00044   }
00045 }
00046 
00047 void electron_pion_separation(TString ifname, Bool_t weightedProbability=false, Bool_t useAntiProbability=false, Bool_t debug=false){
00048   TString folder = "data2011/Analysis/";
00049   cout << ifname << " " << weightedProbability << " " << useAntiProbability << " " << debug << endl;
00050   TString picEx = ".pdf";
00051   Float_t minProb = 1e-6;
00052   Float_t maxProb = 1;
00053   if (weightedProbability){
00054     minProb = 1e-2;
00055     maxProb =10;
00056   }
00057   TFile inputFile(folder + TString(ifname(0,ifname.Sizeof()-5)).Append(".root"),"READ");
00058   if (inputFile.IsOpen()) {
00059     //cout << "Error: " << TString(folder + TString(ifname(0,ifname.Sizeof()-5)).Append(".root")).Data() << " not found." << endl;
00060     //return;
00061     //} 
00062     TFile outputFile(folder + "EP_" + TString(ifname(0,ifname.Sizeof()-5)).Append(".root"),"RECREATE","output of electron_pion_separation.cxx");
00063     if (outputFile.IsOpen()) {
00064       
00065       Int_t nEvents = 100000;
00066       //if (debug) nEvents = 1000;
00067       
00068       const Int_t nChambers = 12;
00069       const Int_t usedSusibos = 8;
00070       const Int_t recoSets = 6;
00071       const Int_t rebinSteps = 1;//5;
00072       const Int_t nMerge = 2;
00073       const Int_t nFilter = 4;
00074       
00075       TString mergeName[2] = {"NotMerged","Merged"};
00076       TString filterName[4] = {"wOwU","woOwU","woOwoU","wOwoU"}; // w == with; wo == without; O == overflow; U == undershoot
00077       Int_t usedSuId[8] = {                    2,         12,        6,        18,         17,         16,        4,        5  };
00078       //TString recoSetName[9] = {"ClusterInteg","3PadInteg","1PadInteg","ClusterIntegWindow","3PadIntegWindow","1PadIntegWindow","ClusterAmpl","3PadAmpl","1PadAmpl"};
00079       TString recoSetName[6] = {"ClusterInteg","3PadInteg","ClusterIntegWindow","3PadIntegWindow","ClusterAmpl","3PadAmpl"};
00080       
00081       TH1D* hElectronSpectra;//[usedSusibos][recoSets][nMerge][nFilter];      
00082       TH1D* hElectronProbability;//[usedSusibos][recoSets][nMerge][nFilter];
00083       TH1D* hPionSpectra;//[usedSusibos][recoSets][nMerge][nFilter];
00084       TH1D* hPionProbability;//[usedSusibos][recoSets][nMerge][nFilter];      
00085       TProfile* hEfficiency = new TProfile("temp","temp",12,0.5,12.5);//[usedSusibos][recoSets][rebinSteps][nMerge][nFilter];
00086       TLegend* leg[usedSusibos][recoSets][nMerge][nFilter];
00087       TString cHist, cTitle;
00088       TH1D *sLikelihoodPi[nChambers+1];//[usedSusibos][recoSets][nChambers][rebinSteps][nMerge][nFilter];
00089       TH1D *sLikelihoodEl[nChambers+1];//[usedSusibos][recoSets][nChambers][rebinSteps][nMerge][nFilter];
00090       
00091       for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
00092         cHist.Form("%s_Pi_Prob_temp_%i",TString(ifname(0,ifname.Sizeof()-5)).Data(),Chamberloop);
00093         sLikelihoodPi[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
00094         cHist.Form("%s_El_Prob_temp_%i",TString(ifname(0,ifname.Sizeof()-5)).Data(),Chamberloop);
00095         sLikelihoodEl[Chamberloop] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
00096       }
00097       
00098       
00099       TString title, newtitle, EP_temp;
00100       Double_t MeanElectronSpectra(0), MeanPionSpectra(0);
00101       const Int_t rebin = 5;//10;//4;
00102 
00103       Int_t xRes(2400), yRes(800);
00104       if (debug) {
00105         xRes = 1200;
00106         yRes = 400;
00107       }
00108       TCanvas *c = new TCanvas("c","c",xRes,yRes);
00109       c->Divide(2,1);
00110   
00111       xRes = 1200;
00112       yRes = 800;
00113       if (debug) {
00114         xRes = 600;
00115         yRes = 400;
00116       }  
00117       TCanvas *cE = new TCanvas("cE","cE",xRes,yRes);
00118       cE->Divide(1,1);
00119   
00120       xRes = 12000;
00121       yRes = 4000;
00122       if (debug) {
00123         xRes = 1200;
00124         yRes = 400;
00125       }  
00126       TLegend *lCutPiEff = new TLegend(0.15,0.15,0.5,0.5);
00127       TCanvas *cL = new TCanvas("cL","cL",xRes,yRes);
00128       cL->Divide(nChambers,recoSets);
00129       TCanvas *test;
00130       TH1I *Ptest;
00131       TH1I *Etest;
00132       if (debug) {
00133         test = new TCanvas("test","test",1200,800);
00134         test->Divide(2,2); 
00135       }
00136       
00137       // -------------------------------------------------------------------   clone spectra from readtree_spadic() output file
00138       for (Int_t suid = 0; suid < usedSusibos; suid++){
00139         for (Int_t rs = 0; rs < recoSets; rs++) {
00140           for (Int_t imerge = 0; imerge < nMerge; imerge++) {
00141             for (Int_t ifilter = 0; ifilter < nFilter; ifilter++) {
00142               //cout << "Search spectra" << endl;
00143 
00144               title.Form("preCompSpectra%sEl_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00145               newtitle.Form("Spectra%sEl_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00146               hElectronSpectra = (TH1D*)inputFile.Get(title)->Clone(newtitle);
00147               if (rs < 4) {
00148                 //hElectronSpectra[suid][rs]->Rebin(rebin);
00149                 hElectronSpectra->GetXaxis()->SetRangeUser(0,15000);
00150                 hElectronSpectra->GetYaxis()->SetRangeUser(1e-6,1);
00151               }
00152               else {
00153                 hElectronSpectra->GetXaxis()->SetRangeUser(0,1000);
00154                 hElectronSpectra->GetYaxis()->SetRangeUser(1e-6,1);
00155               }
00156               hElectronSpectra->SetLineColor(2);
00157               title.Form("preCompSpectra%sPi_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00158               newtitle.Form("Spectra%sPi_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00159               hPionSpectra = (TH1D*)inputFile.Get(title)->Clone(newtitle);
00160               if (rs < 4) {
00161                 //hPionSpectra[suid][rs]->Rebin(rebin);
00162                 hPionSpectra->GetXaxis()->SetRangeUser(0,15000);
00163                 hPionSpectra->GetYaxis()->SetRangeUser(1e-6,1);
00164               }
00165               else {
00166                 hPionSpectra->GetXaxis()->SetRangeUser(0,1000);
00167                 hPionSpectra->GetYaxis()->SetRangeUser(1e-6,1);
00168               }
00169      
00170               leg[suid][rs][imerge][ifilter] = new TLegend(0.5,0.6,0.85,0.85);
00171               leg[suid][rs][imerge][ifilter]->SetLineColor(0);
00172               leg[suid][rs][imerge][ifilter]->SetLineStyle(0);
00173               leg[suid][rs][imerge][ifilter]->SetFillColor(0);
00174               leg[suid][rs][imerge][ifilter]->SetFillStyle(0);
00175               leg[suid][rs][imerge][ifilter]->SetTextSize(0.035);
00176               if (hElectronSpectra->GetEntries() > 0)
00177                 MeanElectronSpectra = hElectronSpectra->GetMean(1);
00178               if (hPionSpectra->GetEntries() > 0)
00179                 MeanPionSpectra = hPionSpectra->GetMean(1);
00180               EP_temp.Form("#LTQ#GT_{e^{-}} = %.3f",MeanElectronSpectra);
00181               leg[suid][rs][imerge][ifilter]->AddEntry(hElectronSpectra,EP_temp,"l");
00182               EP_temp.Form("#LTQ#GT_{#pi^{-}} = %.3f",MeanPionSpectra);
00183               leg[suid][rs][imerge][ifilter]->AddEntry(hPionSpectra,EP_temp,"l");
00184               if (MeanPionSpectra != 0)
00185                 EP_temp.Form("#LTQ#GT_{e^{-}}/#LTQ#GT_{#pi^{-}} = %.3f",MeanElectronSpectra/MeanPionSpectra);
00186               else
00187                 EP_temp.Form("#LTe^{-}/#LT#pi^{-} = %.3f",0.0);
00188               leg[suid][rs][imerge][ifilter]->AddEntry((TObject*)0,EP_temp,"");
00189 
00190               newtitle.Form("Probability%sPi_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00191               hPionProbability = (TH1D*)hPionSpectra->Clone(newtitle);
00192               hPionProbability->SetXTitle("Deposited charge Q [a.u.]");
00193               hPionProbability->SetYTitle("dN/dQ [a.u.]");
00194               //printf("    %i Pions\n",hPionProbability->GetEntries());
00195 
00196               newtitle.Form("Probability%sEl_%i_%s_%s",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data());
00197               hElectronProbability =  (TH1D*)hElectronSpectra->Clone(newtitle);
00198               hElectronProbability->SetXTitle("Deposited charge Q [a.u.]");
00199               hElectronProbability->SetYTitle("dN/dQ [a.u.]");
00200               //printf("    %i Electrons\n",hElectronProbability->GetEntries());
00201 
00202               // normalization to probability to have at least one particle
00203               if (weightedProbability) {
00204                 hPionProbability->Reset();
00205                 hElectronProbability->Reset();
00206                 for (Int_t ibin = 1; ibin < hElectronSpectra->GetNbinsX(); ibin++) {
00207                   Double_t elecP = hElectronSpectra->GetBinContent(ibin);
00208                   Double_t pionP = hPionSpectra->GetBinContent(ibin);
00209                   Double_t sum =  elecP + pionP;
00210                   if (sum != 0) {
00211                     if (debug) printf("E:%.5f  P:%.5f  S:%.5f  E:%.5f  P:%.5f S:%.5f\n",elecP,pionP,sum,elecP/sum,pionP/sum,elecP/sum+pionP/sum);
00212                     hElectronProbability->SetBinContent(ibin, elecP / sum);
00213                     hPionProbability->SetBinContent(ibin, pionP / sum);
00214                   }
00215                   else {
00216                     hElectronProbability->SetBinContent(ibin, 0);
00217                     hPionProbability->SetBinContent(ibin, 0);
00218                   }
00219                 }
00220               }
00221               //cout << "Have spectra" << endl;
00222               for (Int_t irebin = 0; irebin < rebinSteps; irebin++) {
00223 
00224                 printf("<<<<<<<<<<<<  Susibo %i  %s  >>>>>>>>>>>>>>  %s %s %s <<<<<<<<<<<<<<< Rebin %i \n",usedSuId[suid],recoSetName[rs].Data(), ifname.Data(),mergeName[imerge].Data(),filterName[ifilter].Data(), irebin);
00225 
00226                 title.Form("EfficiencySuID%i_%s_%s_%s_Bin%i",usedSuId[suid],recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data(),irebin);
00227                 hEfficiency->Reset();
00228                 hEfficiency->SetNameTitle(title,title);
00229                 hEfficiency->GetYaxis()->SetRangeUser(0.01,100.5);
00230                 hEfficiency->SetMarkerStyle(20);
00231                 hEfficiency->SetXTitle("Number of layer");
00232                 hEfficiency->SetYTitle("Pion efficiency [%]");
00233 
00234                 if (weightedProbability) {
00235                   //rebinSteps = 0;
00236                   //rebin = 1;
00237                 }
00238                 if (rebinSteps == 1) {
00239                   printf("%i nBins",hElectronProbability->GetNbinsX());
00240                   if (rs < 4) {
00241                     //hElectronProbability->Rebin(rebin*rebin);
00242                     //hPionProbability->Rebin(rebin*rebin);
00243                     for (Int_t i = 0; i < 3; i++) {
00244                       hElectronProbability->Rebin(rebin);
00245                       hPionProbability->Rebin(rebin);
00246                     }
00247                   }
00248                   else {
00249                     hElectronProbability->Rebin(rebin);
00250                     hPionProbability->Rebin(rebin);
00251                   }
00252                 }
00253                 printf("  %i nReBins\n",hElectronProbability->GetNbinsX());
00254                 
00255                 title.Form("_%s_Sus%i_%s_%s_Rebin%i",recoSetName[rs].Data(),usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data(),irebin);
00256                 if (debug) {
00257                   Ptest = (TH1I*)hPionSpectra->Clone("Ptest");
00258                   test->cd(4);
00259                   Ptest->DrawCopy();
00260                   Ptest->Reset();
00261                   Ptest->GetXaxis()->SetRangeUser(0,15000);
00262                   test->cd(2);
00263                   Ptest->Draw();
00264                   Etest = (TH1I*)hElectronSpectra->Clone("Etest");
00265                   test->cd(3);
00266                   Etest->DrawCopy();
00267                   Etest->Reset();
00268                   Etest->GetXaxis()->SetRangeUser(0,15000);
00269                   test->cd(1);
00270                   Etest->Draw();
00271                 }
00272                 for (Int_t Chamberloop = 1; Chamberloop <= nChambers; Chamberloop++){
00273                   cHist.Form("sLikelihoodPi%02dSpadic%i_%s_%s_%s_Bin%i",Chamberloop,suid,recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data(),irebin);
00274                   //sLikelihoodPi[suid][rs][Chamberloop][irebin][imerge][ifilter] = new TH1D(cHist,"e and pi Prob",220,-0.05,1.05);
00275                   sLikelihoodPi[Chamberloop]->Reset();
00276                   sLikelihoodPi[Chamberloop]->SetNameTitle(cHist,cHist);
00277                   sLikelihoodPi[Chamberloop]->SetXTitle("Likelihood");
00278                   sLikelihoodPi[Chamberloop]->SetYTitle("#");
00279                   cHist.Form("sLikelihoodEl%02dSpadic%i_%s_%s_%s_Bin%i",Chamberloop,suid,recoSetName[rs].Data(),mergeName[imerge].Data(),filterName[ifilter].Data(),irebin);
00280                   //sLikelihoodEl[suid][rs][Chamberloop][irebin][imerge][ifilter] = new TH1D(cHist,"e and p prob",220,-0.05,1.05);
00281                   sLikelihoodEl[Chamberloop]->Reset();
00282                   sLikelihoodEl[Chamberloop]->SetNameTitle(cHist,cHist);
00283                   sLikelihoodEl[Chamberloop]->SetXTitle("Likelihood");
00284                   sLikelihoodEl[Chamberloop]->SetYTitle("#");
00285 
00286                   Bool_t lowStatistic = false;
00287                   Int_t nElTries(0), nPiTries(0);
00288                   for(Int_t iE = 0;iE < nEvents;iE++){
00289                     if (debug) Statusbar(iE,nEvents);
00290                     Double_t P_e = 1;
00291                     Double_t P_pi = 1;
00292          
00293                     for(int iCh = 1; iCh <= Chamberloop; iCh++){
00294                       if (lowStatistic)
00295                         break;
00296                       // ladung würfeln
00297                       Double_t fEl = 0;
00298                       Double_t fPi = 0;
00299                       while (fEl == 0 || fPi == 0) { //right
00300                         //while (fEl == 0 && fPi == 0) { // wrong
00301                         Double_t chEl = hElectronSpectra->GetRandom();  
00302                         fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chEl));
00303                         fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chEl));
00304                         nElTries++;
00305                         if (nElTries/Float_t(Chamberloop*nEvents) > 15) {
00306                           printf("too low electron statistic\n");
00307                           lowStatistic = true;
00308                           break;
00309                         }
00310         
00311                       }
00312                       P_e = P_e * (fEl * 1e9);
00313                       P_pi = P_pi * (fPi * 1e9);
00314         
00315                     }      
00316                     Double_t L_e = 0;
00317                     if((P_e+P_pi)>0){
00318                       L_e = P_e /(P_e+P_pi);
00319                     }
00320                     sLikelihoodEl[Chamberloop]->Fill(L_e);
00321       
00322                     P_e = 1;
00323                     P_pi = 1;
00324           
00325                     Double_t L_p = 0;
00326                     if (useAntiProbability) {
00327                       for(Int_t iCh = 1;iCh<= Chamberloop;iCh++){
00328                         if (lowStatistic)
00329                           break;
00330                         // ladung würfeln
00331                         Double_t fEl = 0;
00332                         Double_t fPi = 0;
00333                         while (fEl == 0 || fPi == 0) { // right
00334                           //while (fEl == 0 && fPi == 0) { //wrong
00335                           Double_t chPi = hPionSpectra->GetRandom();    
00336                           fEl = hPionProbability->GetBinContent( hPionProbability->FindBin(chPi));
00337                           fPi = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chPi));
00338                           nPiTries++;
00339                           if (nPiTries/Float_t(Chamberloop*nEvents) > 15) {
00340                             printf("too low pion statistic\n");
00341                             lowStatistic = true;
00342                             break;
00343                           }
00344         
00345                         }
00346                         P_e = P_e * (fEl * 1e9);
00347                         P_pi = P_pi * (fPi * 1e9);
00348         
00349                       }      
00350                       L_p = 0;
00351                       if((P_e+P_pi)>0){
00352                         L_p = P_pi /(P_e+P_pi);
00353                       }  
00354                     } 
00355                     else {
00356                       for(Int_t iCh = 1;iCh<= Chamberloop;iCh++){
00357                         if (lowStatistic)
00358                           break;
00359                         // ladung würfeln
00360                         Double_t fEl = 0;
00361                         Double_t fPi = 0;
00362                         while (fEl == 0 || fPi == 0) { // right
00363                           // while (fEl == 0 && fPi == 0) { // wrong
00364                           Double_t chPi = hPionSpectra->GetRandom();    
00365                           fPi = hPionProbability->GetBinContent( hPionProbability->FindBin(chPi));
00366                           fEl = hElectronProbability->GetBinContent( hElectronProbability->FindBin(chPi));
00367                           nPiTries++;
00368                           if (nPiTries/Float_t(Chamberloop*nEvents) > 15) {
00369                             printf("too low pion statistic\n");
00370                             lowStatistic = true;
00371                             break;
00372                           }
00373         
00374                         }
00375                         P_e = P_e * (fEl * 1e9);
00376                         P_pi = P_pi * (fPi * 1e9);
00377               
00378                       }    
00379                       L_p = 0; 
00380                       if((P_e+P_pi)>0){
00381                         L_p = 1. - (P_pi /(P_e+P_pi));
00382                       }   
00383                     }
00384                     if (lowStatistic){
00385                       printf(" ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! too low statistic\n");
00386                       break;
00387                     }
00388                     else
00389                       sLikelihoodPi[Chamberloop]->Fill(L_p);
00390                   }
00391                   if (lowStatistic)
00392                     printf("%2i %.2f %.2f\n",Chamberloop,nElTries/Float_t(Chamberloop*nEvents),nPiTries/Float_t(Chamberloop*nEvents));
00393                   Double_t ElectronIntegral =  0.0;
00394                   Double_t PionIntegral = 0.0;
00395                   Int_t ElectronBin = 220;        
00396                   Double_t ElectronIntAll = sLikelihoodEl[Chamberloop]->Integral(0,220);
00397                   Double_t PionIntAll = sLikelihoodPi[Chamberloop]->Integral(0,220);
00398           
00399                   while (ElectronIntegral < 0.90000 && ElectronBin > 0){
00400                     ElectronBin -= 1;
00401                     ElectronIntegral = (sLikelihoodEl[Chamberloop]->Integral(ElectronBin,220))/ElectronIntAll;
00402                     PionIntegral     = (sLikelihoodPi[Chamberloop]->Integral(ElectronBin,220))/PionIntAll;
00403                   }
00404                   if (debug) {
00405                     printf("  layer %d\n",Chamberloop);
00406                     printf("    Pioneneffizienz     : %8.4f Bin %.3f\n",PionIntegral*100,(ElectronBin)/220.);
00407                     printf("    remaining Pieff.    : %8.4f Bin %.3f\n",PionIntegral*100*PionIntAll/(PionIntegral*PionIntAll+ElectronIntegral*ElectronIntAll),(ElectronBin)/220.);
00408                     printf("    Electroneneffizienz : %8.4f Bin %.3f\n\n",ElectronIntegral*100,(ElectronBin)/220.);
00409                   }
00410           
00411                   sLikelihoodEl[Chamberloop]->SetLineColor(2);
00412                   sLikelihoodEl[Chamberloop]->SetLineStyle(kDashed);
00413                   sLikelihoodEl[Chamberloop]->GetYaxis()->SetRangeUser(0.1,10*nEvents);
00414                   sLikelihoodPi[Chamberloop]->SetLineColor(1);
00415                   sLikelihoodPi[Chamberloop]->GetYaxis()->SetRangeUser(0.1,10*nEvents);
00416 
00417             
00418                   TLegend *legend = new TLegend(.15,.77,.4,.86);
00419                   legend->SetLineColor(0);
00420                   legend->SetFillColor(0);
00421                   legend->SetFillStyle(0);
00422                   //legend->AddEntry(sLikelihoodEl[Chamberloop],"electron","l");
00423                   //legend->AddEntry(sLikelihoodPi[Chamberloop],"pion","l");
00424                   //legend->Draw();
00425 
00426                   cTitle.Form("Chamberstack: %2d Pioneneffizienz:  %.2f  Elektroneneffizienz: %.2f Likelihood: %.2f",Chamberloop,PionIntegral*100,ElectronIntegral*100,(ElectronBin)/220.);
00427                   sLikelihoodEl[Chamberloop]->SetTitle(cTitle);
00428  
00429                   Double_t x1 = ElectronBin/220.;
00430                   Double_t y1 = 0.0;
00431                   Double_t x2 = ElectronBin/220.;
00432                   Double_t y2 = 
00433                     sLikelihoodEl[Chamberloop]->GetBinContent(sLikelihoodEl[Chamberloop]->FindBin(ElectronBin/220))
00434                     +sLikelihoodEl[Chamberloop]->GetBinContent(sLikelihoodEl[Chamberloop]->FindBin(ElectronBin/220+2*1.1/220))
00435                     +sLikelihoodEl[Chamberloop]->GetBinContent(sLikelihoodEl[Chamberloop]->FindBin(ElectronBin/220-2*1.1/220))
00436                     +sLikelihoodEl[Chamberloop]->GetBinContent(sLikelihoodEl[Chamberloop]->FindBin(ElectronBin/220+4*1.1/220))
00437                     +sLikelihoodEl[Chamberloop]->GetBinContent(sLikelihoodEl[Chamberloop]->FindBin(ElectronBin/220-4*1.1/220));
00438                   TLine *Likelihood = new TLine(x1,y1,x2,y2);
00439                   Likelihood->SetLineColor(kRed);   
00440             
00441                   // c1->Update();
00442                   hEfficiency->Fill(Chamberloop,PionIntegral*100.);
00443                   //cL->cd(rs * nChambers + Chamberloop)->SetLogy(1);
00444                   //sLikelihoodEl[Chamberloop]->DrawCopy("");   
00445                   //sLikelihoodPi[Chamberloop]->DrawCopy("same");
00446                   //legend->Draw("same");
00447                   //Likelihood->Draw("same");
00448                   //cL->Update();
00449                   if (Chamberloop == 12) {
00450                     if (hPionSpectra->GetEntries() > 0) c->cd(1)->SetLogy(1);
00451                     c->cd(1)->SetTickx(1);
00452                     c->cd(1)->SetTicky(1);
00453                     if (debug) {
00454                       hPionSpectra->DrawCopy();
00455                       hElectronSpectra->DrawCopy("same");
00456                     }
00457                     //if (debug) {
00458                     if (rs < 4) {
00459                       hPionProbability->GetXaxis()->SetRangeUser(0,15000);
00460                       hPionProbability->GetYaxis()->SetRangeUser(minProb,maxProb);
00461                     }
00462                     else {
00463                       hPionProbability->GetXaxis()->SetRangeUser(0,1000);
00464                       hPionProbability->GetYaxis()->SetRangeUser(minProb,maxProb);
00465                     } 
00466                     hPionProbability->DrawCopy();
00467                     if (rs < 4) {
00468                       hElectronProbability->GetXaxis()->SetRangeUser(0,15000);
00469                       hElectronProbability->GetYaxis()->SetRangeUser(minProb,maxProb);
00470                     }
00471                     else {
00472                       hElectronProbability->GetXaxis()->SetRangeUser(0,1000);
00473                       hElectronProbability->GetYaxis()->SetRangeUser(minProb,maxProb);
00474                     } 
00475                     hElectronProbability->DrawCopy("same");
00476                     //}
00477                     legend->Draw("same");
00478                     leg[suid][rs][imerge][ifilter]->Draw("same");
00479 
00480                     if (sLikelihoodEl[Chamberloop]->GetEntries() > 0) c->cd(2)->SetLogy(1);
00481                     c->cd(2)->SetTickx(1);
00482                     c->cd(2)->SetTicky(1);
00483                     c->cd(2)->SetLogy(1);
00484                     hEfficiency->DrawCopy("P,text0");
00485         
00486                     c->Update();
00487                     c->SaveAs(folder + "Pics/EP/" + TString(ifname(0,ifname.Sizeof()-5)) + title + picEx);
00488 
00489                     outputFile.cd();
00490                     c->Write(TString(ifname(0,ifname.Sizeof()-5)) + title, TObject::kOverwrite);
00491                   }
00492                   //outputFile.cd();
00493                   //sLikelihoodPi[suid][rs][Chamberloop][irebin][imerge][ifilter]->Write("", TObject::kOverwrite);
00494                   //sLikelihoodEl[suid][rs][Chamberloop][irebin][imerge][ifilter]->Write("", TObject::kOverwrite);
00495                 } // for Chamberloop
00496                 cE->cd(1)->SetLogy(1);
00497                 hEfficiency->SetLineColor(rs+1);
00498                 hEfficiency->SetMarkerColor(rs+1);
00499                 if (suid == 0 && irebin == 0 && imerge == 0 && ifilter == 0)
00500                   lCutPiEff->AddEntry(hEfficiency,recoSetName[rs].Data(),"l");
00501                 if (rs == 0)
00502                   hEfficiency->DrawCopy("P");
00503                 else {
00504                   hEfficiency->DrawCopy("P same");
00505                   lCutPiEff->Draw("same");
00506                 }
00507                 cE->Update();
00508                 outputFile.cd();
00509                 hEfficiency->Write("", TObject::kOverwrite);
00510                 hPionProbability->Write("", TObject::kOverwrite);
00511                 hElectronProbability->Write("", TObject::kOverwrite);
00512 
00513                 if (rebinSteps > 1) {
00514                   hElectronProbability->Rebin(rebin);
00515                   if (weightedProbability)
00516                     hElectronProbability->Scale(1./rebin);
00517                   hPionProbability->Rebin(rebin);
00518                   if (weightedProbability)
00519                     hPionProbability->Scale(1./rebin);
00520                 }
00521 
00522                 title.Form("Sus_%i_%s_%s_Rebin_%i",usedSuId[suid],mergeName[imerge].Data(),filterName[ifilter].Data(),irebin);
00523                 cE->SaveAs(folder + "Pics/EP/PionEff/" + TString(ifname(0,ifname.Sizeof()-5)) + title + "_PionEff" + picEx);
00524                 //cL->SaveAs(folder + "Pics/EP/Likelihood/" + TString(ifname(0,ifname.Sizeof()-5)) + title + "_Likelihood" + picEx);
00525                 outputFile.cd();
00526                 cE->Write(TString(ifname(0,ifname.Sizeof()-5)) + title + "_PionEff", TObject::kOverwrite);
00527                 //cL->Write(TString(ifname(0,ifname.Sizeof()-5)) + title + "_Likelihood", TObject::kOverwrite);
00528               } //for rebin
00529             } // for ifilter
00530           } // for imerge
00531         } // for recroset
00532       }// for susid
00533       
00534       outputFile.cd();
00535       if(outputFile.IsOpen())
00536         outputFile.Close();
00537       else 
00538         cout << "output already closed?!?" << endl;
00539       
00540     } // if output.isopen()
00541     
00542     inputFile.cd();
00543     if(inputFile.IsOpen())
00544       inputFile.Close();
00545     else 
00546       cout << "input already closed?!?" << endl;
00547     
00548   } // if input.isopen()
00549   else {    
00550     cout << "Error: " << TString(folder + TString(ifname(0,ifname.Sizeof()-5)).Append(".root")).Data() << " not found." << endl;
00551   } 
00552 
00553 }