StRoot: StHighptPool/macros/backgroundDca.C Source File

00001 #include "commonmacro/common.h"
00002 #include "common/Name.cc"
00003 #include "commonmacro/histutil.h"
00004 
00005 int show=0;
00006 
00007 void backgroundDca(const char* inName=
00008             "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
00009             const char* psDir="ps",
00010             int cut = 1,
00011             const char* outDir="./",
00012             const char* more = "west",
00013             float extraValue = 0)
00014 {
00015   cout << "--------------------------" << endl;
00016   cout << "in name=" << inName << endl
00017        << "ps dir=" << psDir << endl
00018        << "cut=" << cut << endl;
00019   cout << "--------------------------" << endl;
00020 
00021   TFile* inRoot[2]={0};
00022  
00023   // real
00024   inRoot[0] = new TFile(inName);
00025 
00026   // mc
00027   inRoot[1] = new TFile(more);
00028 
00029   if(!inRoot[0] || !inRoot[1]){
00030     cout << "cannot find input file" << endl; return;
00031   }
00032 
00033   gSystem->Load("StHiMicroAnalysis");
00034   Cut::SetCut(cut);
00035   Cut::ShowCuts();
00036 
00037   if(!inRoot){
00038     cout << "cannot find the infile" << endl;
00039     return;
00040   }
00041 
00042   // deduce the output root file directory from the input file
00043   TString sOutDir=inName;
00044   sOutDir.Remove(sOutDir.Last('/'),sOutDir.Length());
00045   cout << sOutDir << endl;
00046 
00047   TText* t = new TText;
00048   float ptAry[][100] = 
00049     {
00050       {1.5,2,2.5,3,4,5,6},
00051       {1.7, 1.8, 1.9, 2.0, 2.2, 2.45, 2.7, 3.0, 3.35, 3.8, 4.4, 5.1, 6.0}
00052     };
00053   char* ptType[] = { "PtPr","PtGl"};
00054   char* chargeType[]= { "","Pos","Neg"};
00055   char* dcaType[] = { "DcaXYGl","DcaGl", "SDcaGl"};
00056   char* ptRebin[] = { "", "Rebin" };
00057   int nDcaRebin[] = {2,2,2};
00058   float bkgMin[] = {-3, 2, -3 };
00059   float bkgMax[] = {-2, 3, -2 };
00060   float dcaCut = fabs(Cut::mSDcaGl[0]);
00061   float sigMin[] = {-dcaCut,0,-dcaCut};
00062   float sigMax[] = {dcaCut,dcaCut,dcaCut};
00063 
00064   TCanvas c1("c1","c1",400,500);
00065   TCanvas c2("c2","c2",400,500);
00066 
00067   TH2* h2[2], *h2b[2], *h2temp[2]; TH1* h1a[2], *h1b[2];
00068   TH2* h2Real; TH2* h2Mc, *h2Contam;
00069   TH1* h1Real; TH1* h1Mc, *h1Contam;
00070 
00071   int nDcaType=3;
00072   int nptType=1;
00073   int npt[] = {6,12};
00074   int nc=3;
00075   int nx[]={2,3}; int ny[]={3,4}; // per ptRebin
00076 
00077   gStyle->SetOptStat(0); gStyle->SetOptLogy(1);
00078 
00079   
00080   //--------------------------------
00081   for(int iPtRebin=0; iPtRebin<2; iPtRebin++){
00082 
00083     for(int iDcaType=0; iDcaType<nDcaType; iDcaType++){
00084       for(int ic=0;ic<nc;ic++){
00085         if(ic<2){
00086           sprintf(name,"%s.h2%sPtPrRebin",sPM[ic].Data(),dcaType[iDcaType]);
00087 
00088           // real
00089           h2Real=(TH2*)inRoot[0]->Get(name);
00090 
00091           // mc
00092           h2Mc=(TH2*)inRoot[1]->Get(name);
00093 
00094           sprintf(name,"%s.h2Contam%sPtPrRebin",
00095                   sPM[ic].Data(),dcaType[iDcaType]); 
00096           // contam
00097           h2Contam=(TH2*)inRoot[1]->Get(name);
00098 
00099         }
00100         else{ // sum up the charges
00101           // real
00102           sprintf(name,"Plus.h2%sPtPrRebin",
00103                   dcaType[iDcaType]);
00104           h2Real=(TH2*)inRoot[0]->Get(name);
00105           sprintf(name,"Minus.h2%sPtPrRebin",
00106                   dcaType[iDcaType],ptRebin[iPtRebin]);
00107           h2temp=(TH2*)inRoot[0]->Get(name);
00108           h2Real->Add(h2temp);
00109 
00110           // mc
00111           sprintf(name,"Plus.h2%sPtPrRebin",
00112                   dcaType[iDcaType]);
00113           h2Mc=(TH2*)inRoot[1]->Get(name);
00114           sprintf(name,"Minus.h2%sPtPrRebin",
00115                   dcaType[iDcaType]);
00116           h2temp=(TH2*)inRoot[1]->Get(name);
00117           h2Mc->Add(h2temp);
00118 
00119           // contam
00120           sprintf(name,"Plus.h2Contam%sPtPrRebin",
00121                   dcaType[iDcaType]);
00122           h2Contam=(TH2*)inRoot[1]->Get(name);
00123           sprintf(name,"Minus.h2Contam%sPtPrRebin",
00124                   dcaType[iDcaType]);
00125           h2temp=(TH2*)inRoot[1]->Get(name);
00126           h2Contam->Add(h2temp);
00127         }
00128         cout <<"### " << h2Real->GetName() << endl;
00129 
00130         // background and signal
00131         sprintf(title,"bkg + signal%s %s (cut %d)",
00132                 dcaType[iDcaType],chargeType[ic],cut);
00133         Divide(&c1,nx[iPtRebin],ny[iPtRebin],title,inName);
00134 
00135         c2.Clear();
00136         
00137         
00138         sprintf(title,"h1Background%s%s%s",
00139                 dcaType[iDcaType],chargeType[ic],ptRebin[iPtRebin]);
00140 
00141         TH1D* hb=new TH1D(title,title,npt[iPtRebin],ptAry[iPtRebin]);
00142 
00143         for(int ipt=0;ipt<npt[iPtRebin];ipt++){
00144           c1.cd(ipt+1);
00145           h1Real=HistSlice(h2Real,"","",0,
00146                            ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],
00147                            "x","e");
00148           h1Mc=HistSlice(h2Mc,"","",0,
00149                          ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],"x");
00150           h1Contam=HistSlice(h2Contam,"","",0,
00151                              ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],"x");
00152 
00153           
00154           if(nDcaRebin[iDcaType]>1){
00155             h1Real->Rebin(nDcaRebin[iDcaType]);
00156             h1Mc->Rebin(nDcaRebin[iDcaType]);
00157             h1Contam->Rebin(nDcaRebin[iDcaType]);
00158           }
00159           // check bin sizes
00160           if(h1Real->GetNbinsX() != h1Mc->GetNbinsX()){
00161             cout << "different bin widths? " << endl; return;
00162           }
00163 
00164           //      if(ipt==0) { dump(h1Mc); dump(h1Contam); }
00165           /*
00166           if(iDcaType==1){
00167             int bin = h1b->GetXaxis()->FindBin(.99999);
00168             cout << ptAry[iPtRebin][ipt] << "-" << ptAry[iPtRebin][ipt+1]
00169                  << ":\t" << h1a->Integral(1,bin)
00170                  << "\t" << h1b->Integral(1,bin)
00171                  << endl;
00172           }
00173           */
00174           if(show)cout << h1Real->GetTitle() << endl;
00175           
00176           // find the background
00177           int lowBin = h1Contam->FindBin(bkgMin[iDcaType]);
00178           int upBin  = h1Contam->FindBin(bkgMax[iDcaType]-.0000001);
00179                
00180           TAxis* axis=h1Contam->GetXaxis();
00181           if(show)
00182             cout << "bkg region : " << axis->GetBinLowEdge(lowBin) << " -- " 
00183                  << axis->GetBinUpEdge(upBin) << endl;
00184 
00185           // background yield in bad region
00186           float badBkgCount = h1Contam->Integral(lowBin,upBin);
00187           if(show)cout << "\tbad bkg count=" << badBkgCount << endl;
00188 
00189           if(badBkgCount<1) continue;
00190 
00191           // all counts in bad region
00192           float allBadCount = h1Mc->Integral(lowBin,upBin);
00193           
00194           if(show)cout << "\tall bad count=" << allBadCount << endl;
00195 
00196           // fraction of bkg in bad region
00197           float fractionBkg = badBkgCount/allBadCount;
00198           
00199           if(show)cout << "\tfraction bkg=" << fractionBkg << endl;
00200 
00201           // find background in signal region
00202           lowBin = h1Contam->FindBin(sigMin[iDcaType]);
00203           upBin  = h1Contam->FindBin(sigMax[iDcaType]-.0000001);
00204           
00205           if(show)cout << "signal region : " 
00206                        << axis->GetBinLowEdge(lowBin) << " -- " 
00207                        << axis->GetBinUpEdge(upBin) << endl;
00208 
00209           // bkg in signal region
00210           
00211           float signalBkgCount = h1Contam->Integral(lowBin,upBin);
00212           if(show)cout << "\tsignal bkg count=" << signalBkgCount << endl;
00213           
00214           float signalBkgOverBkg = signalBkgCount/badBkgCount;
00215           
00216           if(show)cout << "\tsignal background/background=" 
00217                        << signalBkgOverBkg << endl;
00218 
00219           // now look at real data
00220           lowBin = h1Real->FindBin(bkgMin[iDcaType]);
00221           upBin  = h1Real->FindBin(bkgMax[iDcaType]-.0000001);
00222 
00223           axis = h1Real->GetXaxis();
00224           if(show)cout << "real bkg region : " 
00225                        << axis->GetBinLowEdge(lowBin) << " -- " 
00226                        << axis->GetBinUpEdge(upBin) << endl;
00227           
00228           // count in bad region
00229           float realBadCount =h1Real->Integral(lowBin,upBin);
00230           if(show)cout << "\treal bad count=" << realBadCount << endl;
00231 
00232           if(realBadCount<1) continue;
00233 
00234           // estimate background in signal region
00235           float realSignalBkgCount = realBadCount*fractionBkg*signalBkgOverBkg;
00236           if(show)cout << "\treal signal bkg=" << realSignalBkgCount << endl;
00237           
00238 
00239           lowBin = h1Real->FindBin(sigMin[iDcaType]);
00240           upBin  = h1Real->FindBin(sigMax[iDcaType]-.0000001);
00241 
00242           float bkgPlusSignal = h1Real->Integral(lowBin,upBin);
00243           if(show)cout << "real bkg+sig=" << bkgPlusSignal << endl;
00244 
00245           float bkgOverData=realSignalBkgCount/bkgPlusSignal;
00246 
00247           if(show)cout << "fraction of bkg under signal region real="
00248                        << bkgOverData << endl;
00249 
00250           hb->SetBinContent(ipt+1,bkgOverData);
00251           
00252 
00253           h1Contam->SetFillColor(17);
00254 
00255           if(h1Mc->Integral())h1Mc->Scale(1./h1Mc->Integral());
00256           if(h1Contam->Integral())h1Contam->Scale(1./h1Contam->Integral());
00257           if(h1Real->Integral())h1Real->Scale(1./h1Real->Integral());
00258 
00259           h1Mc->Draw();
00260           h1Contam->Draw("same");
00261           h1Real->SetMarkerStyle(8); h1Real->SetMarkerSize(0.6); 
00262           h1Real->Draw("esame");
00263 
00264           c2.cd();
00265           hb->Draw();
00266           
00267             
00268         }
00269         Print(&c2,psDir,hb->GetName());
00270         sprintf(name,"contam%s%s%s",
00271                 dcaType[iDcaType],chargeType[ic],ptRebin[iPtRebin]);
00272         Print(&c1,psDir,name);
00273       } // ic
00274     } // iDcaType
00275     
00276   }
00277 }