StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
backgroundDca.C
1 #include "commonmacro/common.h"
2 #include "common/Name.cc"
3 #include "commonmacro/histutil.h"
4 
5 int show=0;
6 
7 void backgroundDca(const char* inName=
8  "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
9  const char* psDir="ps",
10  int cut = 1,
11  const char* outDir="./",
12  const char* more = "west",
13  float extraValue = 0)
14 {
15  cout << "--------------------------" << endl;
16  cout << "in name=" << inName << endl
17  << "ps dir=" << psDir << endl
18  << "cut=" << cut << endl;
19  cout << "--------------------------" << endl;
20 
21  TFile* inRoot[2]={0};
22 
23  // real
24  inRoot[0] = new TFile(inName);
25 
26  // mc
27  inRoot[1] = new TFile(more);
28 
29  if(!inRoot[0] || !inRoot[1]){
30  cout << "cannot find input file" << endl; return;
31  }
32 
33  gSystem->Load("StHiMicroAnalysis");
34  Cut::SetCut(cut);
35  Cut::ShowCuts();
36 
37  if(!inRoot){
38  cout << "cannot find the infile" << endl;
39  return;
40  }
41 
42  // deduce the output root file directory from the input file
43  TString sOutDir=inName;
44  sOutDir.Remove(sOutDir.Last('/'),sOutDir.Length());
45  cout << sOutDir << endl;
46 
47  TText* t = new TText;
48  float ptAry[][100] =
49  {
50  {1.5,2,2.5,3,4,5,6},
51  {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}
52  };
53  char* ptType[] = { "PtPr","PtGl"};
54  char* chargeType[]= { "","Pos","Neg"};
55  char* dcaType[] = { "DcaXYGl","DcaGl", "SDcaGl"};
56  char* ptRebin[] = { "", "Rebin" };
57  int nDcaRebin[] = {2,2,2};
58  float bkgMin[] = {-3, 2, -3 };
59  float bkgMax[] = {-2, 3, -2 };
60  float dcaCut = fabs(Cut::mSDcaGl[0]);
61  float sigMin[] = {-dcaCut,0,-dcaCut};
62  float sigMax[] = {dcaCut,dcaCut,dcaCut};
63 
64  TCanvas c1("c1","c1",400,500);
65  TCanvas c2("c2","c2",400,500);
66 
67  TH2* h2[2], *h2b[2], *h2temp[2]; TH1* h1a[2], *h1b[2];
68  TH2* h2Real; TH2* h2Mc, *h2Contam;
69  TH1* h1Real; TH1* h1Mc, *h1Contam;
70 
71  int nDcaType=3;
72  int nptType=1;
73  int npt[] = {6,12};
74  int nc=3;
75  int nx[]={2,3}; int ny[]={3,4}; // per ptRebin
76 
77  gStyle->SetOptStat(0); gStyle->SetOptLogy(1);
78 
79 
80  //--------------------------------
81  for(int iPtRebin=0; iPtRebin<2; iPtRebin++){
82 
83  for(int iDcaType=0; iDcaType<nDcaType; iDcaType++){
84  for(int ic=0;ic<nc;ic++){
85  if(ic<2){
86  sprintf(name,"%s.h2%sPtPrRebin",sPM[ic].Data(),dcaType[iDcaType]);
87 
88  // real
89  h2Real=(TH2*)inRoot[0]->Get(name);
90 
91  // mc
92  h2Mc=(TH2*)inRoot[1]->Get(name);
93 
94  sprintf(name,"%s.h2Contam%sPtPrRebin",
95  sPM[ic].Data(),dcaType[iDcaType]);
96  // contam
97  h2Contam=(TH2*)inRoot[1]->Get(name);
98 
99  }
100  else{ // sum up the charges
101  // real
102  sprintf(name,"Plus.h2%sPtPrRebin",
103  dcaType[iDcaType]);
104  h2Real=(TH2*)inRoot[0]->Get(name);
105  sprintf(name,"Minus.h2%sPtPrRebin",
106  dcaType[iDcaType],ptRebin[iPtRebin]);
107  h2temp=(TH2*)inRoot[0]->Get(name);
108  h2Real->Add(h2temp);
109 
110  // mc
111  sprintf(name,"Plus.h2%sPtPrRebin",
112  dcaType[iDcaType]);
113  h2Mc=(TH2*)inRoot[1]->Get(name);
114  sprintf(name,"Minus.h2%sPtPrRebin",
115  dcaType[iDcaType]);
116  h2temp=(TH2*)inRoot[1]->Get(name);
117  h2Mc->Add(h2temp);
118 
119  // contam
120  sprintf(name,"Plus.h2Contam%sPtPrRebin",
121  dcaType[iDcaType]);
122  h2Contam=(TH2*)inRoot[1]->Get(name);
123  sprintf(name,"Minus.h2Contam%sPtPrRebin",
124  dcaType[iDcaType]);
125  h2temp=(TH2*)inRoot[1]->Get(name);
126  h2Contam->Add(h2temp);
127  }
128  cout <<"### " << h2Real->GetName() << endl;
129 
130  // background and signal
131  sprintf(title,"bkg + signal%s %s (cut %d)",
132  dcaType[iDcaType],chargeType[ic],cut);
133  Divide(&c1,nx[iPtRebin],ny[iPtRebin],title,inName);
134 
135  c2.Clear();
136 
137 
138  sprintf(title,"h1Background%s%s%s",
139  dcaType[iDcaType],chargeType[ic],ptRebin[iPtRebin]);
140 
141  TH1D* hb=new TH1D(title,title,npt[iPtRebin],ptAry[iPtRebin]);
142 
143  for(int ipt=0;ipt<npt[iPtRebin];ipt++){
144  c1.cd(ipt+1);
145  h1Real=HistSlice(h2Real,"","",0,
146  ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],
147  "x","e");
148  h1Mc=HistSlice(h2Mc,"","",0,
149  ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],"x");
150  h1Contam=HistSlice(h2Contam,"","",0,
151  ptAry[iPtRebin][ipt],ptAry[iPtRebin][ipt+1],"x");
152 
153 
154  if(nDcaRebin[iDcaType]>1){
155  h1Real->Rebin(nDcaRebin[iDcaType]);
156  h1Mc->Rebin(nDcaRebin[iDcaType]);
157  h1Contam->Rebin(nDcaRebin[iDcaType]);
158  }
159  // check bin sizes
160  if(h1Real->GetNbinsX() != h1Mc->GetNbinsX()){
161  cout << "different bin widths? " << endl; return;
162  }
163 
164  // if(ipt==0) { dump(h1Mc); dump(h1Contam); }
165  /*
166  if(iDcaType==1){
167  int bin = h1b->GetXaxis()->FindBin(.99999);
168  cout << ptAry[iPtRebin][ipt] << "-" << ptAry[iPtRebin][ipt+1]
169  << ":\t" << h1a->Integral(1,bin)
170  << "\t" << h1b->Integral(1,bin)
171  << endl;
172  }
173  */
174  if(show)cout << h1Real->GetTitle() << endl;
175 
176  // find the background
177  int lowBin = h1Contam->FindBin(bkgMin[iDcaType]);
178  int upBin = h1Contam->FindBin(bkgMax[iDcaType]-.0000001);
179 
180  TAxis* axis=h1Contam->GetXaxis();
181  if(show)
182  cout << "bkg region : " << axis->GetBinLowEdge(lowBin) << " -- "
183  << axis->GetBinUpEdge(upBin) << endl;
184 
185  // background yield in bad region
186  float badBkgCount = h1Contam->Integral(lowBin,upBin);
187  if(show)cout << "\tbad bkg count=" << badBkgCount << endl;
188 
189  if(badBkgCount<1) continue;
190 
191  // all counts in bad region
192  float allBadCount = h1Mc->Integral(lowBin,upBin);
193 
194  if(show)cout << "\tall bad count=" << allBadCount << endl;
195 
196  // fraction of bkg in bad region
197  float fractionBkg = badBkgCount/allBadCount;
198 
199  if(show)cout << "\tfraction bkg=" << fractionBkg << endl;
200 
201  // find background in signal region
202  lowBin = h1Contam->FindBin(sigMin[iDcaType]);
203  upBin = h1Contam->FindBin(sigMax[iDcaType]-.0000001);
204 
205  if(show)cout << "signal region : "
206  << axis->GetBinLowEdge(lowBin) << " -- "
207  << axis->GetBinUpEdge(upBin) << endl;
208 
209  // bkg in signal region
210 
211  float signalBkgCount = h1Contam->Integral(lowBin,upBin);
212  if(show)cout << "\tsignal bkg count=" << signalBkgCount << endl;
213 
214  float signalBkgOverBkg = signalBkgCount/badBkgCount;
215 
216  if(show)cout << "\tsignal background/background="
217  << signalBkgOverBkg << endl;
218 
219  // now look at real data
220  lowBin = h1Real->FindBin(bkgMin[iDcaType]);
221  upBin = h1Real->FindBin(bkgMax[iDcaType]-.0000001);
222 
223  axis = h1Real->GetXaxis();
224  if(show)cout << "real bkg region : "
225  << axis->GetBinLowEdge(lowBin) << " -- "
226  << axis->GetBinUpEdge(upBin) << endl;
227 
228  // count in bad region
229  float realBadCount =h1Real->Integral(lowBin,upBin);
230  if(show)cout << "\treal bad count=" << realBadCount << endl;
231 
232  if(realBadCount<1) continue;
233 
234  // estimate background in signal region
235  float realSignalBkgCount = realBadCount*fractionBkg*signalBkgOverBkg;
236  if(show)cout << "\treal signal bkg=" << realSignalBkgCount << endl;
237 
238 
239  lowBin = h1Real->FindBin(sigMin[iDcaType]);
240  upBin = h1Real->FindBin(sigMax[iDcaType]-.0000001);
241 
242  float bkgPlusSignal = h1Real->Integral(lowBin,upBin);
243  if(show)cout << "real bkg+sig=" << bkgPlusSignal << endl;
244 
245  float bkgOverData=realSignalBkgCount/bkgPlusSignal;
246 
247  if(show)cout << "fraction of bkg under signal region real="
248  << bkgOverData << endl;
249 
250  hb->SetBinContent(ipt+1,bkgOverData);
251 
252 
253  h1Contam->SetFillColor(17);
254 
255  if(h1Mc->Integral())h1Mc->Scale(1./h1Mc->Integral());
256  if(h1Contam->Integral())h1Contam->Scale(1./h1Contam->Integral());
257  if(h1Real->Integral())h1Real->Scale(1./h1Real->Integral());
258 
259  h1Mc->Draw();
260  h1Contam->Draw("same");
261  h1Real->SetMarkerStyle(8); h1Real->SetMarkerSize(0.6);
262  h1Real->Draw("esame");
263 
264  c2.cd();
265  hb->Draw();
266 
267 
268  }
269  Print(&c2,psDir,hb->GetName());
270  sprintf(name,"contam%s%s%s",
271  dcaType[iDcaType],chargeType[ic],ptRebin[iPtRebin]);
272  Print(&c1,psDir,name);
273  } // ic
274  } // iDcaType
275 
276  }
277 }
Generated by   doxygen 1.8.5