StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plTpcEffCorrect.C
1 
2 void plTpcEffCorrect(bool print=true){
3  gStyle->SetPalette(1);
4  gStyle->SetOptStat(100000000);
5  const float PI=TMath::Pi();
6  TList *Lx; TLine *ln;
7 
8  //input files for TPC effic correction
9  TString core="/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/";
10  TFile *fdQcd = TFile::Open("outMC/tpcEffic/uncorrected/rcf10016.wana.hist.root");
11  //TFile *fdQcd = TFile::Open("outMC/corrected/rcf10016.wana.hist.root");
12  TFile *fdData = TFile::Open(core+"forJoe/data/run9setABCD.wana.hist.root");
13  TFile *fdCorrected = TFile::Open("/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/feb01/data/run9setABCD.wana.hist.root");
14 
15  //#if 0
16  //plots for pt>1
17  //hQcd=(TH2F*)fdQcd->Get("muTr2D1"); assert(hQcd);
18  //hData=(TH2F*)fdData->Get("muTr2D1"); assert(hData);
19 
20  //plots for pt>5
21  hQcd=(TH2F*)fdQcd->Get("muTr2D1pt5"); assert(hQcd);
22  hData=(TH2F*)fdData->Get("muTr2D1pt5"); assert(hData);
23  hCorr=(TH2F*)fdCorrected->Get("muTr2D1pt5weight"); assert(hCorr);
24 
25  //get plots rebinned in my (eta,phi) mapping
26  hQcdRebin=(TH2F*) hQcd->Clone();
27  hDataRebin=(TH2F*) hData->Clone();
28  hCorrRebin=(TH2F*) hCorr->Clone();
29  hQcdRebin->Rebin2D(10,10); hDataRebin->Rebin2D(10,10);
30  hCorrRebin->Rebin2D(10,10);
31  //hDataRebin->Draw("colz"); hDataRebin->SetMinimum(0);
32 
33 
34  //add etabin lines to original plot
35  Lx=hQcd->GetListOfFunctions();
36  //Lx=hData->GetListOfFunctions();
37  for(int i=0; i<9; i++){
38  float binlim=(i*.22)-.88;
39  ln=new TLine(binlim,-PI,binlim,PI);
40  ln->SetLineColor(kRed);
41  ln->SetLineWidth(2);
42  Lx->Add(ln);
43  }
44  //hQcd->Rebin2D(3,2); hQcd->Draw("colz");
45  //hData->Draw("colz");
46 
47  //ratio of QCD MC to Data
48  hRatio=(TH2F*) hQcdRebin->Clone();
49  hRatio->Divide(hDataRebin);
50  //hRatio->Draw("colz"); hRatio->SetMaximum(.13);
51 
52 #if 0
53  //final slices in |eta| to compare ratio
54  gStyle->SetOptStat(1011);
55  hRatio->ProjectionY("hRatioFin1",1,2);
56  hRatio->ProjectionY("hRatioFin2",3,4);
57  hRatio->ProjectionY("hRatioFin3",5,6);
58  hRatio->ProjectionY("hRatioFin4",7,8);
59  hRatio->ProjectionY("hRatioFin5",9,10);
60  hRatioFin1->SetMinimum(0);hRatioFin1->SetMaximum(0.25); hRatioFin1->SetTitle("Ratio MC/Data tracks eta<-0.66");
61  hRatioFin2->SetMinimum(0);hRatioFin2->SetMaximum(0.25); hRatioFin2->SetTitle("Ratio MC/Data tracks -0.66<eta<-0.22");
62  hRatioFin3->SetMinimum(0); hRatioFin3->SetMaximum(0.25); hRatioFin3->SetTitle("Ratio MC/Data tracks |eta|<0.22");
63  hRatioFin4->SetMinimum(0);hRatioFin4->SetMaximum(0.25); hRatioFin4->SetTitle("Ratio MC/Data tracks 0.22<eta<0.66");
64  hRatioFin5->SetMinimum(0);hRatioFin5->SetMaximum(0.25); hRatioFin5->SetTitle("Ratio MC/Data tracks eta>0.66");
65  c3=new TCanvas("cc","cc",1500,500);
66  c3->Divide(5,1);
67  c3->cd(1);
68  hRatioFin1->Draw();
69  c3->cd(2);
70  hRatioFin2->Draw();
71  c3->cd(3);
72  hRatioFin3->Draw();
73  c3->cd(4);
74  hRatioFin4->Draw();
75  c3->cd(5);
76  hRatioFin5->Draw();
77 #endif
78 
79  //make projections to get ratio of MC/Data for each eta bin
80  TH1 *hRatioAll[10];
81  hRatio->ProjectionY("hRatioA",1,1); hRatioAll[0]=hRatioA;
82  hRatio->ProjectionY("hRatioB",2,2); hRatioAll[1]=hRatioB;
83  hRatio->ProjectionY("hRatioC",3,3); hRatioAll[2]=hRatioC;
84  hRatio->ProjectionY("hRatioH",8,8); hRatioAll[7]=hRatioH;
85  hRatio->ProjectionY("hRatioI",9,9); hRatioAll[8]=hRatioI;
86  hRatio->ProjectionY("hRatioJ",10,10); hRatioAll[9]=hRatioJ;
87 
88  float normal[10]={0.07,0.075,0.08,0,0,0,0,0.08,0.075,0.07};//determined from looking at data where yield was high
89  //float normal[10]={0.049,0.065,0.069,0,0,0,0,0.068,0.062,0.055};//for syst error chosen for lowest phi bin in each eta slice
90  if(print) c=new TCanvas("aa","aa",800,600);
91  for(int j=0; j<10; j++){
92  if(j>2 && j<7) continue;
93  hRatioAll[j]->SetTitle(Form("Eta bin %c: Ratio of tracks QCD MC to run 9 data (track pt > 5)",'A'+j));
94  hRatioAll[j]->SetMinimum(0); hRatioAll[j]->SetLineWidth(2);
95  Lx=hRatioAll[j]->GetListOfFunctions(); hRatioAll[j]->SetStats(false);
96  ln=new TLine(-PI,normal[j],PI,normal[j]);
97  ln->SetLineColor(kRed); ln->SetLineWidth(2);
98  Lx->Add(ln);
99  if(print){
100  hRatioAll[j]->Draw();
101  c->Print(Form("etabin%c.png",'A'+j));
102  }
103  }
104 
105 
106  float phiRad=999;
107  float etaDet=999;
108 
109  //initialize and clear eff values
110  float effic[10][24];
111  float efficSec[10][24];
112  for(int i=0;i<10;i++){
113  for(int l=0;l<24;l++){
114  effic[i][l]=0;
115  efficSec[i][l]=0;
116  }
117  }
118 
119  //calculate inefficiency
120  for(int ieta=0; ieta<10; ieta++){
121  if(ieta>2 && ieta<7) continue;
122  for(int iphi=0; iphi<24; iphi++){
123  float binVal=hRatioAll[ieta]->GetBinContent(iphi+1);
124  effic[ieta][iphi]=binVal/normal[ieta];
125  //cout<<"Eta "<<Form("%c",ieta+'A')<<" phi "<<iphi<<" effic "<<binVal<<endl;
126  phiRad=hRatioAll[ieta]->GetXaxis()->GetBinCenter(iphi+1);
127 
128  //find tpc sector
129  if(ieta<3) etaDet=-1;
130  else if(ieta>6) etaDet=1;
131  const float PI=TMath::Pi();
132  int sec=0;
133  float phi=phiRad/PI*180; // now in degrees
134  if (etaDet>0) { // West TPC
135  float x=75-phi;
136  while(x<0) x+=360;
137  sec=1+(int)( x/30.);
138  } else {
139  float x=phi-105;
140  while(x<0) x+=360;
141  sec=13+(int)( x/30.);
142  }
143 
144  efficSec[ieta][sec-1]+=effic[ieta][iphi]/2; //average 2 bins in sector
145  }
146 #if 0 //print weight arrays to use on data in W algo
147  cout<<"mWeight"<<ieta+1<<"={";
148  for(int isec=0; isec<24; isec++){
149  //cout<<"sector "<<isec+1<<" eta "<<ieta<<" effic "<<efficSec[ieta][isec]<<endl;
150  cout<<efficSec[ieta][isec]<<",";
151  }
152  cout<<"};"<<endl;
153 #endif
154  //set weight arrays to use on MC in W algo
155  for(int isec=0; isec<24; isec++){
156  if(efficSec[ieta][isec]!=0) efficSec[ieta][isec]=1/efficSec[ieta][isec];
157  // cout<<"sector "<<isec+1<<" eta "<<ieta<<" effic "<<efficSec[ieta][isec]<<endl;
158  }
159 
160 
161  }
162 
163  //#endif
164 
165 
166  //Data corrected for efficiencies in outer eta bins
167  //hCorr->Draw("colz");
168  hCorr->Rebin2D(10,10);
169  hRatioCorr=(TH2F*) hQcdRebin->Clone();
170  hRatioCorr->Divide(hCorrRebin);
171  //hRatioCorr->Draw("colz");
172 
173  //add sector across pi boundry separately
174  hCorr->ProjectionX("h1",1,1);
175  hCorr->ProjectionX("h24",24,24);
176  TH1 * hMidCorr[12];
177  hMidCorr[11]=h1;hMidCorr[11]->Add(h24);
178 
179  bool print2=false;
180  if(print2) c2=new TCanvas("bb","bb",800,600);
181  for(int iphi=1; iphi<22; iphi+=2){ // get sector slices
182  hCorr->ProjectionX(Form("h",iphi+1),iphi+1,iphi+2);
183  hMidCorr[(iphi-1)/2]=h;
184 
185  //find tpc sector
186  float phiRad=hCorr->GetYaxis()->GetBinCenter(iphi+1);
187  float etaDet=1;
188  const float PI=TMath::Pi();
189  int sec1=0; int sec2=0;
190  float phi=phiRad/PI*180; // now in degrees
191  float x=75-phi;
192  while(x<0) x+=360;
193  sec1=1+(int)( x/30.);
194  float x=phi-105;
195  while(x<0) x+=360;
196  sec2=13+(int)( x/30.);
197 
198 
199  //#if 0
200  //calculate efficiency at midrapidity
201  float etaVal[10]; float etaEffic[10];
202  for(int ieta=0; ieta<10; ieta++)
203  etaVal[ieta]=hMidCorr[(iphi-1)/2]->GetBinContent(ieta+1);
204  float midRapVal=(etaVal[2]+etaVal[7])/2;
205  if(sec2==13) float midRapVal=etaVal[2];
206  for(int ieta2=3; ieta2<7; ieta2++){
207  etaEffic[ieta2]=etaVal[ieta2]/midRapVal;
208  if(ieta2<5){
209  sec=sec2;
210  //cout<<"sec "<<sec2<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
211  }
212  else{
213  sec=sec1;
214  //cout<<"sec "<<sec1<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
215  }
216  efficSec[ieta2][sec-1]=etaEffic[ieta2];
217  }
218  //#endif
219 
220  if(print2){
221  Lx=h->GetListOfFunctions();
222  h->SetStats(false);
223  ln=new TLine(-0.44,midRapVal,0.44,midRapVal);
224  ln->SetLineColor(kRed); ln->SetLineWidth(2);
225  Lx->Add(ln);
226  h->SetMinimum(0); h->SetFillColor(kBlue);
227  h->Draw();
228  h->SetTitle(Form("Sectors %d and %d Tracks vs Eta",sec1,sec2));
229  //c2->Print(Form("sec%d_%d.png",sec1,sec2));
230  }
231  }
232 
233  //do sector 9,15 separately
234  for(int ieta=0; ieta<10; ieta++)
235  etaVal[ieta]=hMidCorr[11]->GetBinContent(ieta+1);
236  float midRapVal=(etaVal[2]+etaVal[7])/2;
237  if(print2) {
238  hMidCorr[11]->SetTitle("Sectors 9 and 15 Tracks vs Eta");
239  hMidCorr[11]->Draw(); hMidCorr[11]->SetMinimum(0); c2->Print("sec9_15.ps");
240  Lx=h->GetListOfFunctions();
241  ln=new TLine(-PI,midRapVal,PI,midRapVal);
242  ln->SetLineColor(kRed);
243  Lx->Add(ln);
244  }
245  for(int ieta2=3; ieta2<7; ieta2++){
246  etaEffic[ieta2]=etaVal[ieta2]/midRapVal;
247  if(ieta2<5){
248  sec=15;
249  //cout<<"sec "<<sec2<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
250  }
251  else{
252  sec=9;
253  //cout<<"sec "<<sec1<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
254  }
255  efficSec[ieta2][sec-1]=etaEffic[ieta2];
256  }
257 
258 #if 0
259  //print effic numbers
260  for(int jeta=0;jeta<10;jeta++){
261  cout<<"float mWeight"<<jeta+1<<"={";
262  for(int jsec=0;jsec<24;jsec++){
263  //cout<<"sec "<<jsec+1<<" eta "<<jeta<<" effic "<<efficSec[jeta][jsec]<<endl;
264  cout<<efficSec[jeta][jsec]<<",";
265  }
266  cout<<"};"<<endl;
267  }
268 #endif
269 
270 
271 }
Generated by   doxygen 1.8.5