StRoot: StEmcPool/StPhotonAnalysis/dAuBackground.cxx Source File

00001 #include <stdio.h>
00002 #include <Riostream.h>
00003 #include <TFile.h>
00004 #include <TTree.h>
00005 #include <TH2F.h>
00006 #include <TPostScript.h>
00007 #include <TCanvas.h>
00008 #include <TVector3.h>
00009 #include <TPad.h>
00010 #include <TStyle.h>
00011 #include <TString.h>
00012 #include <TPaveText.h>
00013 #include <TMath.h>
00014 
00015 #include <StEmcPool/StPhotonCommon/MyEvent.h>
00016 #include <StEmcPool/StPhotonCommon/MyPoint.h>
00017 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
00018 
00019 #include "AnaCuts.h"
00020 
00021 #include "dAuBackground.h"
00022 
00023 ClassImp(dAuBackground)
00024 
00025 Int_t isBad[2400];
00026 
00027 dAuBackground::dAuBackground()
00028 {
00029     cuts = 0;
00030 
00031   gStyle->SetPadGridX(kFALSE);
00032   gStyle->SetPadGridY(kFALSE);
00033 
00034   for(Int_t i=0;i<2400;i++)
00035     {
00036       isBad[i]=0;
00037     }
00038   Int_t badIds[]={82,156,240,264,733,756,818,853,931,1289,1292,1355,1359,1527,1759,1943,2011,2025,2090,2152,2369,-1};
00039   Int_t j=0;
00040   while(badIds[j]>0)
00041     {
00042       isBad[badIds[j]-1]=1;
00043       j++;
00044     }
00045 
00046     cuts = new AnaCuts();
00047   cout<<"constructed!!"<<endl;
00048 }
00049 dAuBackground::~dAuBackground()
00050 {
00051     if (cuts) delete cuts; cuts = 0;
00052   cout<<"destructed!!"<<endl;
00053 }
00054 
00055 void dAuBackground::run(const char *file)
00056 {
00057   //0-20
00058   TH1F *ratio_mb_020=new TH1F("ratio_mb_020","2001&&2003, FTPC-Au:0-20%", 25,0.,1.);
00059   TH1F *ratio_ht1_020=new TH1F("ratio_ht1_020","2201, FTPC-Au:0-20%", 25,0.,1.);
00060   TH1F *ratio_ht2_020=new TH1F("ratio_ht2_020","2202, FTPC-Au:0-20%", 25,0.,1.);
00061   TH1F *ratio_sn_020=new TH1F("ratio_sn_020","2002, FTPC-Au:0-20%", 25,0.,1.);
00062 
00063   //20-40
00064   TH1F *ratio_mb_2040=new TH1F("ratio_mb_2040","2001&&2003, FTPC-Au:20-40%", 25,0.,1.);
00065   TH1F *ratio_ht1_2040=new TH1F("ratio_ht1_2040","2201, FTPC-Au:20-40%", 25,0.,1.);
00066   TH1F *ratio_ht2_2040=new TH1F("ratio_ht2_2040","2202, FTPC-Au:20-40%", 25,0.,1.);
00067   TH1F *ratio_sn_2040=new TH1F("ratio_sn_2040","2002, FTPC-Au:20-40%", 25,0.,1.);
00068 
00069   //40-100
00070   TH1F *ratio_mb_40100=new TH1F("ratio_mb_40100","2001&&2003, FTPC-Au:40-100%", 25,0.,1.);
00071   TH1F *ratio_ht1_40100=new TH1F("ratio_ht1_40100","2201, FTPC-Au:40-100%", 25,0.,1.);
00072   TH1F *ratio_ht2_40100=new TH1F("ratio_ht2_40100","2202, FTPC-Au:40-100%", 25,0.,1.);
00073   TH1F *ratio_sn_40100=new TH1F("ratio_sn_40100","2002, FTPC-Au:40-100%", 25,0.,1.);
00074 
00075   TH1F *hfrac_mb=new TH1F("frac_mb","2001&&2003, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.05);
00076   TH1F *hfrac_ht1=new TH1F("frac_ht1","2201, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.5);
00077   TH1F *hfrac_ht2=new TH1F("frac_ht2","2202, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.5);
00078   TH1F *hfrac_sn=new TH1F("frac_sn","2000, I[0.8,1.0]/I[0.0,1.0], 40-100%",50,0.0,0.05);
00079 
00080   //totals
00081   //0-20
00082   TH1F *Ratio_mb_020=new TH1F("Ratio_mb_020","2001&&2003, FTPC-Au:0-20%", 25,0.,1.);
00083   TH1F *Ratio_ht1_020=new TH1F("Ratio_ht1_020","2201, FTPC-Au:0-20%", 25,0.,1.);
00084   TH1F *Ratio_ht2_020=new TH1F("Ratio_ht2_020","2202, FTPC-Au:0-20%", 25,0.,1.);
00085   TH1F *Ratio_sn_020=new TH1F("Ratio_sn_020","2002, FTPC-Au:0-20%", 25,0.,1.);
00086 
00087   //20-40
00088   TH1F *Ratio_mb_2040=new TH1F("Ratio_mb_2040","2001&&2003, FTPC-Au:20-40%", 25,0.,1.);
00089   TH1F *Ratio_ht1_2040=new TH1F("Ratio_ht1_2040","2201, FTPC-Au:20-40%", 25,0.,1.);
00090   TH1F *Ratio_ht2_2040=new TH1F("Ratio_ht2_2040","2202, FTPC-Au:20-40%", 25,0.,1.);
00091   TH1F *Ratio_sn_2040=new TH1F("Ratio_sn_2040","2002, FTPC-Au:20-40%", 25,0.,1.);
00092 
00093   //40-100
00094   TH1F *Ratio_mb_40100=new TH1F("Ratio_mb_40100","2001&&2003, FTPC-Au:40-100%", 25,0.,1.);
00095   TH1F *Ratio_ht1_40100=new TH1F("Ratio_ht1_40100","2201, FTPC-Au:40-100%", 25,0.,1.);
00096   TH1F *Ratio_ht2_40100=new TH1F("Ratio_ht2_40100","2202, FTPC-Au:40-100%", 25,0.,1.);
00097   TH1F *Ratio_sn_40100=new TH1F("Ratio_sn_40100","2002, FTPC-Au:40-100%", 25,0.,1.);
00098 
00099   TH2F *EvsE_mb=new TH2F("EvsE_mb","2001&&2003, neutral Energy vs TPC pT",240,0.,80,240,0.,80.);
00100   TH2F *EvsE_ht1=new TH2F("EvsE_ht1","2201, neutral Energy vs TPC pT",124,0.,80,240,0.,80.);
00101   TH2F *EvsE_ht2=new TH2F("EvsE_ht2","2202, neutral Energy vs TPC pT",124,0.,80,240,0.,80.);
00102   TH2F *EvsE_sn=new TH2F("EvsE_sn","2002, neutral Energy vs TPC pT",240,0.,80,240,0.,80.);
00103 
00104   //spectra
00105   //0-20
00106   TH1F *ptspec_mb_020=new TH1F("ptspec_mb_020","2001&&2003, pT spectrum neutral points, 0-20%",30,0.,15.);
00107   TH1F *ptspec_ht1_020=new TH1F("ptspec_ht1_020","2201, pT spectrum neutral points, 0-20%",30,0.,15.);
00108   TH1F *ptspec_ht2_020=new TH1F("ptspec_ht2_020","2202, pT spectrum neutral points, 0-20%",30,0.,15.);
00109   TH1F *ptspec_sn_020=new TH1F("ptspec_sn_020","2002, pT spectrum neutral points, 0-20%",30,0.,15.);
00110 
00111   TH1F *ptspec_mb_cut_020=new TH1F("ptspec_mb_cut_020","2001&&2003, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
00112   TH1F *ptspec_ht1_cut_020=new TH1F("ptspec_ht1_cut_020","2201, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
00113   TH1F *ptspec_ht2_cut_020=new TH1F("ptspec_ht2_cut_020","2202, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
00114   TH1F *ptspec_sn_cut_020=new TH1F("ptspec_sn_cut_020","2002, pT spectrum neutral points after cut, 0-20%",30,0.,15.);
00115 
00116   TH1F *ptspec_mb_bg_020=new TH1F("ptspec_mb_bg_020","2001&&2003, pT spectrum neutral points, 0-20%",30,0.,15.);
00117   TH1F *ptspec_ht1_bg_020=new TH1F("ptspec_ht1_bg_020","2201, pT spectrum neutral points, 0-20%",30,0.,15.);
00118   TH1F *ptspec_ht2_bg_020=new TH1F("ptspec_ht2_bg_020","2202, pT spectrum neutral points, 0-20%",30,0.,15.);
00119   TH1F *ptspec_sn_bg_020=new TH1F("ptspec_sn_bg_020","2002, pT spectrum neutral points, 0-20%",30,0.,15.);
00120 
00121   //20-40
00122   TH1F *ptspec_mb_2040=new TH1F("ptspec_mb_2040","2001&&2003, pT spectrum neutral points, 20-40%",30,0.,15.);
00123   TH1F *ptspec_ht1_2040=new TH1F("ptspec_ht1_2040","2201, pT spectrum neutral points, 20-40%",30,0.,15.);
00124   TH1F *ptspec_ht2_2040=new TH1F("ptspec_ht2_2040","2202, pT spectrum neutral points, 20-40%",30,0.,15.);
00125   TH1F *ptspec_sn_2040=new TH1F("ptspec_sn_2040","2002, pT spectrum neutral points, 20-40%",30,0.,15.);
00126 
00127   TH1F *ptspec_mb_cut_2040=new TH1F("ptspec_mb_cut_2040","2001&&2003, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
00128   TH1F *ptspec_ht1_cut_2040=new TH1F("ptspec_ht1_cut_2040","2201, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
00129   TH1F *ptspec_ht2_cut_2040=new TH1F("ptspec_ht2_cut_2040","2202, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
00130   TH1F *ptspec_sn_cut_2040=new TH1F("ptspec_sn_cut_2040","2002, pT spectrum neutral points after cut, 20-40%",30,0.,15.);
00131 
00132   TH1F *ptspec_mb_bg_2040=new TH1F("ptspec_mb_bg_2040","2001&&2003, pT spectrum neutral points, 20-40%",30,0.,15.);
00133   TH1F *ptspec_ht1_bg_2040=new TH1F("ptspec_ht1_bg_2040","2201, pT spectrum neutral points, 20-40%",30,0.,15.);
00134   TH1F *ptspec_ht2_bg_2040=new TH1F("ptspec_ht2_bg_2040","2202, pT spectrum neutral points, 20-40%",30,0.,15.);
00135   TH1F *ptspec_sn_bg_2040=new TH1F("ptspec_sn_bg_2040","2002, pT spectrum neutral points, 20-40%",30,0.,15.);
00136 
00137   //40-100
00138   TH1F *ptspec_mb_40100=new TH1F("ptspec_mb_40100","2001&&2003, pT spectrum neutral points, 40-100%",30,0.,15.);
00139   TH1F *ptspec_ht1_40100=new TH1F("ptspec_ht1_40100","2201, pT spectrum neutral points, 40-100%",30,0.,15.);
00140   TH1F *ptspec_ht2_40100=new TH1F("ptspec_ht2_40100","2202, pT spectrum neutral points, 40-100%",30,0.,15.);
00141   TH1F *ptspec_sn_40100=new TH1F("ptspec_sn_40100","2002, pT spectrum neutral points, 40-100%",30,0.,15.);
00142 
00143   TH1F *ptspec_mb_cut_40100=new TH1F("ptspec_mb_cut_40100","2001&&2003, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
00144   TH1F *ptspec_ht1_cut_40100=new TH1F("ptspec_ht1_cut_40100","2201, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
00145   TH1F *ptspec_ht2_cut_40100=new TH1F("ptspec_ht2_cut_40100","2202, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
00146   TH1F *ptspec_sn_cut_40100=new TH1F("ptspec_sn_cut_40100","2002, pT spectrum neutral points after cut, 40-100%",30,0.,15.);
00147 
00148   TH1F *ptspec_mb_bg_40100=new TH1F("ptspec_mb_bg_40100","2001&&2003, pT spectrum neutral points, 40-100%",30,0.,15.);
00149   TH1F *ptspec_ht1_bg_40100=new TH1F("ptspec_ht1_bg_40100","2201, pT spectrum neutral points, 40-100%",30,0.,15.);
00150   TH1F *ptspec_ht2_bg_40100=new TH1F("ptspec_ht2_bg_40100","2202, pT spectrum neutral points, 40-100%",30,0.,15.);
00151   TH1F *ptspec_sn_bg_40100=new TH1F("ptspec_sn_bg_40100","2002, pT spectrum neutral points, 40-100%",30,0.,15.);
00152 
00153 
00154   TH1F *ftpc_mb=new TH1F("ftpc_mb","2001&&2003, ftpc-Au mult.",60,0.,60.);
00155   TH1F *ftpc_ht1=new TH1F("ftpc_ht1","2201, ftpc-Au mult.",60,0.,60.);
00156   TH1F *ftpc_ht2=new TH1F("ftpc_ht2","2202, ftpc-Au mult.",60,0.,60.);
00157   TH1F *ftpc_sn=new TH1F("ftpc_sn","2002, ftpc-Au mult.",60,0.,60.);
00158 
00159   TH1F *ftpc_mb_cut=new TH1F("ftpc_mb_cut","2001&&2003, ftpc-Au mult. for r<.8",60,0.,60.);
00160   TH1F *ftpc_ht1_cut=new TH1F("ftpc_ht1_cut","2201, ftpc-Au mult. for r<.8",60,0.,60.);
00161   TH1F *ftpc_ht2_cut=new TH1F("ftpc_ht2_cut","2202, ftpc-Au mult. for r<.8",60,0.,60.);
00162   TH1F *ftpc_sn_cut=new TH1F("ftpc_sn_cut","2002, ftpc-Au mult. for r<.8",60,0.,60.);
00163 
00164   TH1F *tpc_mb=new TH1F("tpc_mb","2001&&2003, tpc mult.",100,0.,100.);
00165   TH1F *tpc_ht1=new TH1F("tpc_ht1","2201, tpc mult.",100,0.,100.);
00166   TH1F *tpc_ht2=new TH1F("tpc_ht2","2202, tpc mult.",100,0.,100.);
00167   TH1F *tpc_sn=new TH1F("tpc_sn","2002, tpc mult.",100,0.,100.);
00168 
00169   TH1F *tpc_mb_cut=new TH1F("tpc_mb_cut","2001&&2003, tpc mult. for r<.8",100,0.,100.);
00170   TH1F *tpc_ht1_cut=new TH1F("tpc_ht1_cut","2201, tpc mult. for r<.8",100,0.,100.);
00171   TH1F *tpc_ht2_cut=new TH1F("tpc_ht2_cut","2202, tpc mult. for r<.8",100,0.,100.);
00172   TH1F *tpc_sn_cut=new TH1F("tpc_sn_cut","2002, tpc mult. for r<.8",100,0.,100.);
00173 
00174   TH1F *bemc_mb=new TH1F("bemc_mb","2001&&2003, bemc mult.",60,0.,60.);
00175   TH1F *bemc_ht1=new TH1F("bemc_ht1","2201, bemc mult.",60,0.,60.);
00176   TH1F *bemc_ht2=new TH1F("bemc_ht2","2202, bemc mult.",60,0.,60.);
00177   TH1F *bemc_sn=new TH1F("bemc_sn","2002, bemc mult.",60,0.,60.);
00178 
00179   TH1F *bemc_mb_cut=new TH1F("bemc_mb_cut","2001&&2003, bemc mult. for r<.8",60,0.,60.);
00180   TH1F *bemc_ht1_cut=new TH1F("bemc_ht1_cut","2201, bemc mult. for r<.8",60,0.,60.);
00181   TH1F *bemc_ht2_cut=new TH1F("bemc_ht2_cut","2202, bemc mult. for r<.8",60,0.,60.);
00182   TH1F *bemc_sn_cut=new TH1F("bemc_sn_cut","2002, bemc mult. for r<.8",60,0.,60.);
00183 
00184   TH2F *etaphi_fill3085_sn=new TH2F("etaphi_fill3085_sn","eta/phi for 3085 sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
00185   TH2F *etaphi_fill3085_ht2=new TH2F("etaphi_fill3085_ht2","eta/phi for 3085 ht-2 40-100%",500,0.0,1.0,1000,-3.2,3.2);
00186   TH2F *etaphi_fill3088_sn=new TH2F("etaphi_fill3088_sn","eta/phi for 3088 sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
00187   TH2F *etaphi_fill3088_ht2=new TH2F("etaphi_fill3088_ht2","eta/phi for 3088 ht-2 40-100%",500,0.0,1.0,1000,-3.2,3.2);
00188   TH2F *etaphi_fill3088b_ht2=new TH2F("etaphi_fill3088b_ht2","eta/phi for 3088 ht-2 40-100% and R<0.7",500,0.0,1.0,1000,-3.2,3.2);
00189 
00190   TH2F *etaphi_sn=new TH2F("etaphi_sn","eta/phi sn 40-100%",500,0.0,1.0,1000,-3.2,3.2);
00191 
00192 
00193   TH1F *h_trigger=new TH1F("h_trigger","trigger ids after cut",2400,0.5,2400.5);
00194 
00195   ofstream fout;
00196   fout.open("bgfrac.dat");
00197   fout<<"pagina\trunid\tmb\tnmb\tht1\tnht1\tht2\tnht2\tsn\tnsn"<<endl<<endl;
00198 
00199   TFile *inf=new TFile(file,"OPEN");
00200   TTree *eventTree=(TTree*)inf->Get("mEventTree");
00201   
00202   MyEvent *ev=new MyEvent();
00203   eventTree->SetBranchAddress("branch",&ev);
00204 
00205   TPostScript *ps=new TPostScript("dAu_bgperrun.ps",111);
00206   TCanvas *c9=new TCanvas("c9","c9",800,600);
00207   c9->Divide(4,4);
00208 
00209   Int_t entry=0;
00210   Int_t runprev=0;
00211   Int_t pagina=0;
00212   while(eventTree->GetEntry(entry))
00213     {
00214       //if(ev->runId()<4061003)
00215       //{
00216       //entry++;
00217       // continue;
00218       //}
00219       //if(ev->runId()>4061025)
00220       //{
00221       // cout<<"end of this fill!"<<endl;
00222       // break;
00223       //}
00224 
00225       if(ev->runId()!=runprev && runprev!=0)
00226         {
00227           if(ev->runId()<runprev)
00228             {
00229               cout<<"Something wrong!"<<endl;
00230               return;
00231             }
00232           cout<<endl<<"new run: "<<ev->runId()<<endl<<endl;
00233           pagina++;
00234 
00235           Float_t total_mb=ratio_mb_40100->Integral();
00236           Float_t bg_mb=ratio_mb_40100->Integral(21,25);
00237           Float_t total_ht1=ratio_ht1_40100->Integral();
00238           Float_t bg_ht1=ratio_ht1_40100->Integral(21,25);
00239           Float_t total_ht2=ratio_ht2_40100->Integral();
00240           Float_t bg_ht2=ratio_ht2_40100->Integral(21,25);
00241           Float_t total_sn=ratio_sn_40100->Integral();
00242           Float_t bg_sn=ratio_sn_40100->Integral(21,25);
00243 
00244           Float_t frac_mb=total_mb>0 ? bg_mb/total_mb : 0.;
00245           Float_t frac_ht1=total_ht1>0 ? bg_ht1/total_ht1 : 0.;
00246           Float_t frac_ht2=total_ht2>0 ? bg_ht2/total_ht2 : 0.;
00247           Float_t frac_sn=total_sn>0 ? bg_sn/total_sn : 0.;
00248           fout<<pagina<<"\t"<<runprev<<"\t"<<frac_mb<<"\t"<<total_mb<<"\t"<<frac_ht1<<"\t"<<total_ht1<<"\t"<<frac_ht2<<"\t"<<total_ht2<<"\t"<<frac_sn<<"\t"<<total_sn<<endl;
00249 
00250           hfrac_mb->Fill(frac_mb,total_mb);//use weight
00251           hfrac_ht1->Fill(frac_ht1,total_ht1);
00252           hfrac_ht2->Fill(frac_ht2,total_ht2);
00253           hfrac_sn->Fill(frac_sn,total_sn);
00254 
00255           //create new ps page and plot previous run
00256           ps->NewPage();
00257           c9->cd(1);
00258           ratio_mb_020->Draw();
00259           c9->cd(2);
00260           ratio_ht1_020->Draw();
00261           c9->cd(3);
00262           ratio_ht2_020->Draw();
00263           c9->cd(4);
00264           ratio_sn_020->Draw();
00265           c9->cd(5);
00266           ratio_mb_2040->Draw();
00267           c9->cd(6);
00268           ratio_ht1_2040->Draw();
00269           c9->cd(7);
00270           ratio_ht2_2040->Draw();
00271           c9->cd(8);
00272           ratio_sn_2040->Draw();
00273           c9->cd(9);
00274           ratio_mb_40100->Draw();
00275           c9->cd(10);
00276           ratio_ht1_40100->Draw();
00277           c9->cd(11);
00278           ratio_ht2_40100->Draw();
00279           c9->cd(12);
00280           ratio_sn_40100->Draw();
00281           c9->cd(13);
00282           TPaveText *pave=new TPaveText(.1,.1,.9,.9);
00283           char label[20];
00284           sprintf(label,"run %d",runprev);
00285           pave->AddText(label);
00286           pave->Draw();
00287           c9->cd(0);
00288           c9->Update();
00289           
00290           //reset hists
00291           ratio_mb_020->Reset();
00292           ratio_ht1_020->Reset();
00293           ratio_ht2_020->Reset();
00294           ratio_sn_020->Reset();
00295           ratio_mb_2040->Reset();
00296           ratio_ht1_2040->Reset();
00297           ratio_ht2_2040->Reset();
00298           ratio_sn_2040->Reset();
00299           ratio_mb_40100->Reset();
00300           ratio_ht1_40100->Reset();
00301           ratio_ht2_40100->Reset();
00302           ratio_sn_40100->Reset();
00303 
00304           //set runprev
00305           runprev=ev->runId();
00306           
00307         }
00308       else if(ev->runId()!=runprev)
00309         {
00310           cout<<endl<<"first run: "<<ev->runId()<<endl<<endl;
00311           runprev=ev->runId();    
00312         }
00313 
00314       Int_t nTrig=0;
00315       if(ev->runId()==0&&ev->trigger()==0)
00316         {
00317           nTrig=1;
00318           if(ev->highTowerAdc()>cuts->ht1AdcCUT) nTrig+=2;
00319           if(ev->highTowerAdc()>cuts->ht2AdcCUT) nTrig+=4;        
00320           ev->setTrigger(nTrig);
00321         }
00322 
00323 
00324       Bool_t eventOK=kTRUE;
00325       TVector3 vPos=ev->vertex();
00326       //apply EVENT cuts
00327       if(TMath::Abs(vPos.Z())>cuts->vertexzCUT) eventOK=kFALSE;
00328       if(TMath::Abs(vPos.X())>cuts->vertexxCUT) eventOK=kFALSE;
00329       if(TMath::Abs(vPos.Y())>cuts->vertexyCUT) eventOK=kFALSE;
00330       if(vPos.Mag()<.000001) eventOK=kFALSE;
00331       if(ev->trigger()<1 && ev->getMcTrack()->id()==0) eventOK=kFALSE;
00332 
00333       if(ev->trigger()&2||ev->trigger()&4)
00334         {
00335           if(ev->highTowerEnergy()<0.0001) eventOK=kFALSE;
00336           if(isBad[ev->highTowerId()-1]==1) eventOK=kFALSE;
00337         }
00338 
00339 
00340       if(!eventOK)
00341         {
00342           entry++;
00343           continue;
00344         }
00345 
00346       if(ev->trigger()&2||ev->trigger()&4) h_trigger->Fill(ev->highTowerId());
00347 
00348       Int_t t=ev->trigger();
00349       Float_t en=ev->energyInBarrel();
00350       Float_t pt=ev->momentumInTpc();
00351       Float_t ratio=TMath::Abs(en+pt)>0 ? en/(en+pt) : -.1;
00352             
00353       if(ev->refMult()>17)
00354         {
00355           if(t&1) {ratio_mb_020->Fill(ratio); Ratio_mb_020->Fill(ratio);}
00356           if(t&2) {ratio_ht1_020->Fill(ratio); Ratio_ht1_020->Fill(ratio);}
00357           if(t&4) {ratio_ht2_020->Fill(ratio); Ratio_ht2_020->Fill(ratio);}
00358           if(t&8) {ratio_sn_020->Fill(ratio); Ratio_sn_020->Fill(ratio);}
00359         }
00360       else if(ev->refMult()>11)
00361         {
00362           if(t&1) {ratio_mb_2040->Fill(ratio); Ratio_mb_2040->Fill(ratio);}
00363           if(t&2) {ratio_ht1_2040->Fill(ratio); Ratio_ht1_2040->Fill(ratio);}
00364           if(t&4) {ratio_ht2_2040->Fill(ratio); Ratio_ht2_2040->Fill(ratio);}
00365           if(t&8) {ratio_sn_2040->Fill(ratio); Ratio_sn_2040->Fill(ratio);}
00366         }
00367       else if(ev->refMult()>=0)
00368         {
00369 
00370           if(t&1) {ratio_mb_40100->Fill(ratio); Ratio_mb_40100->Fill(ratio);}
00371           if(t&2) {ratio_ht1_40100->Fill(ratio); Ratio_ht1_40100->Fill(ratio);}
00372           if(t&4) {ratio_ht2_40100->Fill(ratio); Ratio_ht2_40100->Fill(ratio);}
00373           if(t&8) {ratio_sn_40100->Fill(ratio); Ratio_sn_40100->Fill(ratio);}
00374         }
00375       
00376       if(t&1){
00377         tpc_mb->Fill(ev->goodPrimaries());
00378         ftpc_mb->Fill(ev->refMult());
00379         bemc_mb->Fill(ev->numberOfPoints());
00380         EvsE_mb->Fill(ev->momentumInTpc(),ev->energyInBarrel());
00381         if(ratio<cuts->ratioCUT){
00382           tpc_mb_cut->Fill(ev->goodPrimaries());
00383           ftpc_mb_cut->Fill(ev->refMult());
00384           bemc_mb_cut->Fill(ev->numberOfPoints());
00385         }
00386       }
00387       if(t&2){
00388         tpc_ht1->Fill(ev->goodPrimaries());
00389         ftpc_ht1->Fill(ev->refMult());
00390         bemc_ht1->Fill(ev->numberOfPoints());
00391         EvsE_ht1->Fill(ev->momentumInTpc(),ev->energyInBarrel());
00392         if(ratio<cuts->ratioCUT){
00393           tpc_ht1_cut->Fill(ev->goodPrimaries());
00394           ftpc_ht1_cut->Fill(ev->refMult());
00395           bemc_ht1_cut->Fill(ev->numberOfPoints());
00396         }
00397       } 
00398       if(t&4){
00399         tpc_ht2->Fill(ev->goodPrimaries());
00400         ftpc_ht2->Fill(ev->refMult());
00401         bemc_ht2->Fill(ev->numberOfPoints());
00402         EvsE_ht2->Fill(ev->momentumInTpc(),ev->energyInBarrel());
00403         if(ratio<cuts->ratioCUT){
00404           tpc_ht2_cut->Fill(ev->goodPrimaries());
00405           ftpc_ht2_cut->Fill(ev->refMult());
00406           bemc_ht2_cut->Fill(ev->numberOfPoints());
00407         }
00408       } 
00409       if(t&8){
00410         tpc_sn->Fill(ev->goodPrimaries());
00411         ftpc_sn->Fill(ev->refMult());
00412         bemc_sn->Fill(ev->numberOfPoints());
00413         EvsE_sn->Fill(ev->momentumInTpc(),ev->energyInBarrel());
00414         if(ratio<cuts->ratioCUT){
00415           tpc_sn_cut->Fill(ev->goodPrimaries());
00416           ftpc_sn_cut->Fill(ev->refMult());
00417           bemc_sn_cut->Fill(ev->numberOfPoints());
00418         }
00419       }
00420  
00421       //get neutral points
00422       MyPoint *p;
00423       TClonesArray *clA=(TClonesArray*)ev->getPointArray();
00424       for(Int_t j=0;j<clA->GetEntries();j++)
00425         {
00426           p=(MyPoint*)clA->At(j);
00427           TVector3 pPos=p->position();
00428           TVector3 pMom=pPos-vPos;
00429           pMom.SetMag(p->energy());
00430         
00431           if(p->distanceToTrack()<cuts->neutralCUT) continue;
00432           if(p->nHitsEta()<cuts->etaHitsCUT) continue;
00433           if(p->nHitsPhi()<cuts->phiHitsCUT) continue;
00434           if(pPos.PseudoRapidity()<cuts->etaMinCUT) continue;
00435           if(pPos.PseudoRapidity()>cuts->etaMaxCUT) continue;
00436 
00437           Float_t pT=pMom.Pt();
00438           if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT)
00439             {
00440               if(ev->refMult()>17)
00441                 {
00442                   if(t&1)
00443                     {
00444                       ptspec_mb_020->Fill(pT);
00445                       if(ratio>=cuts->ratioCUT) 
00446                         ptspec_mb_bg_020->Fill(pT);
00447                       else
00448                         ptspec_mb_cut_020->Fill(pT);
00449                     }
00450                   if(t&2)
00451                     {
00452                       ptspec_ht1_020->Fill(pT);
00453                       if(ratio>=cuts->ratioCUT)
00454                         ptspec_ht1_bg_020->Fill(pT);
00455                       else
00456                         ptspec_ht1_cut_020->Fill(pT);
00457                     }
00458                   if(t&4)
00459                     {
00460                       ptspec_ht2_020->Fill(pT);
00461                       if(ratio>=cuts->ratioCUT)
00462                         ptspec_ht2_bg_020->Fill(pT);
00463                       else
00464                         ptspec_ht2_cut_020->Fill(pT);
00465                     }
00466                   if(t&8)
00467                     {
00468                       ptspec_sn_020->Fill(pT);
00469                       if(ratio>=cuts->ratioCUT)
00470                         ptspec_sn_bg_020->Fill(pT);
00471                       else
00472                         ptspec_sn_cut_020->Fill(pT);
00473                     }
00474                 }
00475               else if(ev->refMult()>11)
00476                 {
00477                   if(t&1)
00478                     {
00479                       ptspec_mb_2040->Fill(pT);
00480                       if(ratio>=cuts->ratioCUT)
00481                         ptspec_mb_bg_2040->Fill(pT);
00482                       else
00483                         ptspec_mb_cut_2040->Fill(pT);
00484                     }
00485                   if(t&2)
00486                     {
00487                       ptspec_ht1_2040->Fill(pT);
00488                       if(ratio>=cuts->ratioCUT)
00489                         ptspec_ht1_bg_2040->Fill(pT);
00490                       else
00491                         ptspec_ht1_cut_2040->Fill(pT);
00492                     }
00493                   if(t&4)
00494                     {
00495                       ptspec_ht2_2040->Fill(pT);
00496                       if(ratio>=cuts->ratioCUT)
00497                         ptspec_ht2_bg_2040->Fill(pT);
00498                       else
00499                         ptspec_ht2_cut_2040->Fill(pT);
00500                     }
00501                   if(t&8)
00502                     {
00503                       ptspec_sn_2040->Fill(pT);
00504                       if(ratio>=cuts->ratioCUT)
00505                         ptspec_sn_bg_2040->Fill(pT);
00506                       else
00507                         ptspec_sn_cut_2040->Fill(pT);
00508                     }
00509                 }
00510               else if(ev->refMult()>=0)
00511                 {
00512                   if(t&1)
00513                     {
00514                       ptspec_mb_40100->Fill(pT);
00515                       if(ratio>=cuts->ratioCUT)
00516                         ptspec_mb_bg_40100->Fill(pT);
00517                       else
00518                         ptspec_mb_cut_40100->Fill(pT);
00519                     }
00520                   if(t&2)
00521                     {
00522                       ptspec_ht1_40100->Fill(pT);
00523                       if(ratio>=cuts->ratioCUT)
00524                         ptspec_ht1_bg_40100->Fill(pT);
00525                       else
00526                         ptspec_ht1_cut_40100->Fill(pT);
00527                     }
00528                   if(t&4)
00529                     {
00530                       ptspec_ht2_40100->Fill(pT);
00531                       if(ratio<0.7) etaphi_fill3088b_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
00532                       if(ratio>=cuts->ratioCUT)
00533                         {
00534                         ptspec_ht2_bg_40100->Fill(pT);
00535                         if(ev->runId()>=4060055&&ev->runId()<=4060088) etaphi_fill3085_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
00536                         if(ev->runId()>=4061003&&ev->runId()<=4061025) etaphi_fill3088_ht2->Fill(pPos.PseudoRapidity(),pPos.Phi());
00537                         }
00538                       else
00539                         ptspec_ht2_cut_40100->Fill(pT);
00540                     }
00541                   if(t&8)
00542                     {
00543                       ptspec_sn_40100->Fill(pT);
00544                       if(ratio>=cuts->ratioCUT)
00545                         {
00546                           ptspec_sn_bg_40100->Fill(pT);
00547                           etaphi_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
00548                           if(ev->runId()>=4060055&&ev->runId()<=4060088) etaphi_fill3085_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
00549                           if(ev->runId()>=4061003&&ev->runId()<=4061025) etaphi_fill3088_sn->Fill(pPos.PseudoRapidity(),pPos.Phi());
00550                         }
00551                       else
00552                         ptspec_sn_cut_40100->Fill(pT);
00553                     }
00554                 }
00555 
00556             }
00557 
00558         }
00559       
00560       entry++;
00561     }
00562   
00563   fout.close();
00564   
00565   ps->Close();
00566 
00567   TString oname(file);
00568   oname.Prepend("outputBg113005_");
00569   TFile o(oname.Data(),"RECREATE");
00570 
00571   etaphi_fill3085_ht2->Write();
00572   etaphi_fill3088_ht2->Write();
00573   etaphi_fill3085_sn->Write();
00574   etaphi_fill3088_sn->Write();
00575   etaphi_fill3088b_ht2->Write();
00576   etaphi_sn->Write();
00577 
00578   hfrac_mb->Write();
00579   hfrac_ht1->Write();
00580   hfrac_ht2->Write();
00581   hfrac_sn->Write();  
00582   h_trigger->Write();
00583 
00584   Ratio_mb_020->Write();
00585   Ratio_ht1_020->Write();
00586   Ratio_ht2_020->Write();
00587   Ratio_sn_020->Write();
00588   Ratio_mb_2040->Write();
00589   Ratio_ht1_2040->Write();
00590   Ratio_ht2_2040->Write();
00591   Ratio_sn_2040->Write();
00592   Ratio_mb_40100->Write();
00593   Ratio_ht1_40100->Write();
00594   Ratio_ht2_40100->Write();
00595   Ratio_sn_40100->Write();
00596 
00597   EvsE_mb->Write();
00598   EvsE_ht1->Write();
00599   EvsE_ht2->Write();
00600   EvsE_sn->Write();
00601 
00602   ptspec_mb_020->Write();
00603   ptspec_mb_cut_020->Write();
00604   ptspec_mb_bg_020->Write();
00605   ptspec_ht1_020->Write();
00606   ptspec_ht1_cut_020->Write();
00607   ptspec_ht1_bg_020->Write();
00608   ptspec_ht2_020->Write();
00609   ptspec_ht2_cut_020->Write();
00610   ptspec_ht2_bg_020->Write();
00611   ptspec_sn_020->Write();
00612   ptspec_sn_cut_020->Write();
00613   ptspec_sn_bg_020->Write();
00614 
00615   ptspec_mb_2040->Write();
00616   ptspec_mb_cut_2040->Write();
00617   ptspec_mb_bg_2040->Write();
00618   ptspec_ht1_2040->Write();
00619   ptspec_ht1_cut_2040->Write();
00620   ptspec_ht1_bg_2040->Write();
00621   ptspec_ht2_2040->Write();
00622   ptspec_ht2_cut_2040->Write();
00623   ptspec_ht2_bg_2040->Write();
00624   ptspec_sn_2040->Write();
00625   ptspec_sn_cut_2040->Write();
00626   ptspec_sn_bg_2040->Write();
00627 
00628   ptspec_mb_40100->Write();
00629   ptspec_mb_cut_40100->Write();
00630   ptspec_mb_bg_40100->Write();
00631   ptspec_ht1_40100->Write();
00632   ptspec_ht1_cut_40100->Write();
00633   ptspec_ht1_bg_40100->Write();
00634   ptspec_ht2_40100->Write();
00635   ptspec_ht2_cut_40100->Write();
00636   ptspec_ht2_bg_40100->Write();
00637   ptspec_sn_40100->Write();
00638   ptspec_sn_cut_40100->Write();
00639   ptspec_sn_bg_40100->Write();
00640 
00641   tpc_mb->Write();
00642   tpc_mb_cut->Write();
00643   tpc_ht1->Write();
00644   tpc_ht1_cut->Write();
00645   tpc_ht2->Write();
00646   tpc_ht2_cut->Write();
00647   tpc_sn->Write();
00648   tpc_sn_cut->Write();
00649   ftpc_mb->Write();
00650   ftpc_mb_cut->Write();
00651   ftpc_ht1->Write();
00652   ftpc_ht1_cut->Write();
00653   ftpc_ht2->Write();
00654   ftpc_ht2_cut->Write();
00655   ftpc_sn->Write();
00656   ftpc_sn_cut->Write();
00657   bemc_mb->Write();
00658   bemc_mb_cut->Write();
00659   bemc_ht1->Write();
00660   bemc_ht1_cut->Write();
00661   bemc_ht2->Write();
00662   bemc_ht2_cut->Write();
00663   bemc_sn->Write();
00664   bemc_sn_cut->Write();
00665 
00666   o.Close();
00667 
00668 }
00669 
00670 
00671 void dAuBackground::runSim(const char *file)
00672 { 
00673   TH1F *ratio_pythia=new TH1F("ratio_pythia","pythia, E over pT", 25,0.,1.);
00674   TH1F *energy_pythia=new TH1F("energy_pythia","energy spectrum",200,0.,100.);
00675   TH1F *tpcpt_pythia=new TH1F("tpcpt_pythia","TPC pT spectrum",200,0.,100.);
00676   TH1F *spec_pythia=new TH1F("spec_pythia","pT spectrum neutral points",30,0.,15.);
00677   TH1F *spec_pythia_bg=new TH1F("spec_pythia_bg","pT spectrum neutral points",30,0.,15.);
00678 
00679   TH1F *ratio_pythia_ht1=new TH1F("ratio_pythia_ht1","pythia, E over pT ht1", 25,0.,1.);
00680   TH1F *energy_pythia_ht1=new TH1F("energy_pythia_ht1","energy spectrum ht1",200,0.,100.);
00681   TH1F *tpcpt_pythia_ht1=new TH1F("tpcpt_pythia_ht1","TPC pT spectrum ht1",200,0.,100.);
00682   TH1F *spec_pythia_ht1=new TH1F("spec_pythia_ht1","pT spectrum neutral points ht1",30,0.,15.);
00683   TH1F *spec_pythia_bg_ht1=new TH1F("spec_pythia_bg_ht1","pT spectrum neutral points ht1",30,0.,15.);
00684   TH1F *spec_pythia_ht1_cut=new TH1F("spec_pythia_ht1_cut","pT spectrum neutral points ht1 after cut",30,0.,15.);
00685 
00686   TH1F *pythia_tpc_ht1=new TH1F("pythia_tpc_ht1","tpc mult. ht-1 sim",100,0.,100.);
00687   TH1F *pythia_tpc_ht1_cut=new TH1F("pythia_tpc_ht1_cut","tpc mult. ht-1 sim after cut",100,0.,100.);
00688   TH1F *pythia_bemc_ht1=new TH1F("pythia_bemc_ht1","bemc mult. ht-1 sim",60,0.,60.);  
00689   TH1F *pythia_bemc_ht1_cut=new TH1F("pythia_bemc_ht1_cut","bemc mult. ht-1 sim after cut",60,0.,60.);
00690 
00691   TFile *inf=new TFile(file,"OPEN");
00692   TTree *eventTree=(TTree*)inf->Get("mEventTree");
00693   MyEvent *ev=new MyEvent();
00694   eventTree->SetBranchAddress("branch",&ev);
00695 
00696   TPostScript *ps=new TPostScript("dAu_bgPythia.ps",111);
00697   TCanvas *c9=new TCanvas("c9","c9",800,1000);
00698   c9->Divide(2,4);  
00699 
00700   Int_t entry=0;
00701   while(eventTree->GetEntry(entry))
00702     {
00703       if(entry%10000==0) cout<<entry<<endl;
00704 
00705       Bool_t eventOK=kTRUE;
00706       TVector3 vPos=ev->vertex();
00707 
00708       //apply EVENT cuts
00709       if(TMath::Abs(vPos.Z())>cuts->vertexzCUT) eventOK=kFALSE;
00710       if(TMath::Abs(vPos.X())>cuts->vertexxCUT) eventOK=kFALSE;
00711       if(TMath::Abs(vPos.Y())>cuts->vertexyCUT) eventOK=kFALSE;
00712       if(vPos.Mag()<.000001) eventOK=kFALSE;
00713 
00714       if(!eventOK)
00715         {
00716           entry++;
00717           continue;
00718         }
00719   
00720       Float_t en=ev->energyInBarrel();
00721       Float_t pt=ev->momentumInTpc();
00722       Float_t ratio=TMath::Abs(en+pt)>0 ? en/(en+pt) : 2.;
00723 
00724       ratio_pythia->Fill(ratio);
00725       energy_pythia->Fill(en);
00726       tpcpt_pythia->Fill(pt);
00727       if(ev->highTowerAdc()>cuts->ht1AdcCUT)
00728         {
00729           pythia_tpc_ht1->Fill(ev->goodPrimaries());
00730           pythia_bemc_ht1->Fill(ev->numberOfPoints());
00731           ratio_pythia_ht1->Fill(ratio);
00732           energy_pythia_ht1->Fill(en);
00733           tpcpt_pythia_ht1->Fill(pt);
00734           if(ratio<cuts->ratioCUT)
00735             {
00736               pythia_tpc_ht1_cut->Fill(ev->goodPrimaries());
00737               pythia_bemc_ht1_cut->Fill(ev->numberOfPoints());
00738             }
00739         }
00740       
00741 
00742       //get neutral points
00743       MyPoint *p;
00744       TClonesArray *clA=(TClonesArray*)ev->getPointArray();
00745       for(Int_t j=0;j<clA->GetEntries();j++)
00746         {
00747           p=(MyPoint*)clA->At(j);
00748           TVector3 pPos=p->position();
00749           TVector3 pMom=pPos-vPos;
00750           pMom.SetMag(p->energy());
00751 
00752           if(p->distanceToTrack()<cuts->neutralCUT) continue;
00753           if(p->nHitsEta()<cuts->etaHitsCUT) continue;
00754           if(p->nHitsPhi()<cuts->phiHitsCUT) continue;
00755           if(pPos.PseudoRapidity()<cuts->etaMinCUT) continue;
00756           if(pPos.PseudoRapidity()>cuts->etaMaxCUT) continue;
00757 
00758           Float_t pT=pMom.Pt();
00759           if(pMom.PseudoRapidity()>cuts->rapidityMinCUT && pMom.PseudoRapidity()<cuts->rapidityMaxCUT)
00760             {
00761               spec_pythia->Fill(pT);
00762               if(ratio>=cuts->ratioCUT) spec_pythia_bg->Fill(pT);
00763               if(ev->highTowerAdc()>cuts->ht1AdcCUT)
00764                 {
00765                   spec_pythia_ht1->Fill(pT);
00766                   if(ratio>=cuts->ratioCUT) spec_pythia_bg_ht1->Fill(pT);
00767                   else spec_pythia_ht1_cut->Fill(pT);
00768                 }
00769             }
00770 
00771 
00772 
00773         }
00774     
00775       entry++;
00776     }
00777 
00778   TFile f("hijing_tryout.root","RECREATE");
00779   spec_pythia_ht1->Write();
00780   spec_pythia_ht1_cut->Write();
00781   spec_pythia_bg_ht1->Write();
00782   ratio_pythia_ht1->Write();
00783   
00784   pythia_tpc_ht1->Write();
00785   pythia_tpc_ht1_cut->Write();
00786   pythia_bemc_ht1->Write();
00787   pythia_bemc_ht1_cut->Write();
00788   f.Close();
00789 
00790   ps->NewPage();
00791   c9->cd(1);
00792   ratio_pythia->Draw();
00793   c9->cd(2);
00794   ratio_pythia_ht1->Draw();
00795   c9->cd(3);
00796   energy_pythia->Draw();
00797   gPad->SetLogy();
00798   c9->cd(4);
00799   energy_pythia_ht1->Draw();
00800   gPad->SetLogy();
00801   c9->cd(5);
00802   tpcpt_pythia->Draw();
00803   gPad->SetLogy();
00804   c9->cd(6);
00805   tpcpt_pythia_ht1->Draw();
00806   gPad->SetLogy();
00807   c9->cd(7);
00808   spec_pythia->Draw();
00809   spec_pythia_bg->SetLineColor(6);
00810   spec_pythia_bg->Draw("same");
00811   gPad->SetLogy();
00812   c9->cd(8);
00813   spec_pythia_ht1->Draw();
00814   spec_pythia_bg_ht1->SetLineColor(6);
00815   spec_pythia_bg_ht1->Draw("same");
00816   gPad->SetLogy();
00817   c9->cd(0);
00818   c9->Update();
00819   ps->Close();
00820 }