StRoot: StSpinPool/StWalgoB2009/macros/plTpcEffCorrect.C Source File

00001 
00002 void plTpcEffCorrect(bool print=true){
00003   gStyle->SetPalette(1);
00004   gStyle->SetOptStat(100000000);
00005   const float PI=TMath::Pi();
00006   TList *Lx;  TLine *ln;
00007 
00008   //input files for TPC effic correction
00009   TString core="/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/";
00010   TFile *fdQcd = TFile::Open("outMC/tpcEffic/uncorrected/rcf10016.wana.hist.root");
00011   //TFile *fdQcd = TFile::Open("outMC/corrected/rcf10016.wana.hist.root");
00012   TFile *fdData = TFile::Open(core+"forJoe/data/run9setABCD.wana.hist.root");
00013   TFile *fdCorrected = TFile::Open("/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/feb01/data/run9setABCD.wana.hist.root");
00014 
00015   //#if 0
00016   //plots for pt>1
00017   //hQcd=(TH2F*)fdQcd->Get("muTr2D1"); assert(hQcd);
00018   //hData=(TH2F*)fdData->Get("muTr2D1"); assert(hData);
00019 
00020   //plots for pt>5
00021   hQcd=(TH2F*)fdQcd->Get("muTr2D1pt5"); assert(hQcd);
00022   hData=(TH2F*)fdData->Get("muTr2D1pt5"); assert(hData);
00023   hCorr=(TH2F*)fdCorrected->Get("muTr2D1pt5weight"); assert(hCorr);
00024 
00025   //get plots rebinned in my (eta,phi) mapping
00026   hQcdRebin=(TH2F*) hQcd->Clone();
00027   hDataRebin=(TH2F*) hData->Clone();
00028   hCorrRebin=(TH2F*) hCorr->Clone();
00029   hQcdRebin->Rebin2D(10,10); hDataRebin->Rebin2D(10,10);
00030   hCorrRebin->Rebin2D(10,10);
00031   //hDataRebin->Draw("colz"); hDataRebin->SetMinimum(0);
00032   
00033 
00034   //add etabin lines to original plot
00035   Lx=hQcd->GetListOfFunctions();
00036   //Lx=hData->GetListOfFunctions();
00037   for(int i=0; i<9; i++){
00038     float binlim=(i*.22)-.88;
00039     ln=new TLine(binlim,-PI,binlim,PI);
00040     ln->SetLineColor(kRed);
00041     ln->SetLineWidth(2);
00042     Lx->Add(ln);
00043   }
00044   //hQcd->Rebin2D(3,2); hQcd->Draw("colz");
00045   //hData->Draw("colz");
00046 
00047   //ratio of QCD MC to Data
00048   hRatio=(TH2F*) hQcdRebin->Clone();
00049   hRatio->Divide(hDataRebin);
00050   //hRatio->Draw("colz"); hRatio->SetMaximum(.13);
00051   
00052 #if 0
00053   //final slices in |eta| to compare ratio
00054   gStyle->SetOptStat(1011);
00055   hRatio->ProjectionY("hRatioFin1",1,2);
00056   hRatio->ProjectionY("hRatioFin2",3,4);
00057   hRatio->ProjectionY("hRatioFin3",5,6);
00058   hRatio->ProjectionY("hRatioFin4",7,8);
00059   hRatio->ProjectionY("hRatioFin5",9,10);
00060   hRatioFin1->SetMinimum(0);hRatioFin1->SetMaximum(0.25); hRatioFin1->SetTitle("Ratio MC/Data tracks eta<-0.66");
00061   hRatioFin2->SetMinimum(0);hRatioFin2->SetMaximum(0.25); hRatioFin2->SetTitle("Ratio MC/Data tracks -0.66<eta<-0.22");
00062   hRatioFin3->SetMinimum(0); hRatioFin3->SetMaximum(0.25); hRatioFin3->SetTitle("Ratio MC/Data tracks |eta|<0.22");
00063   hRatioFin4->SetMinimum(0);hRatioFin4->SetMaximum(0.25); hRatioFin4->SetTitle("Ratio MC/Data tracks 0.22<eta<0.66");
00064   hRatioFin5->SetMinimum(0);hRatioFin5->SetMaximum(0.25); hRatioFin5->SetTitle("Ratio MC/Data tracks eta>0.66");
00065   c3=new TCanvas("cc","cc",1500,500);
00066   c3->Divide(5,1);
00067   c3->cd(1);
00068   hRatioFin1->Draw();
00069   c3->cd(2);
00070   hRatioFin2->Draw();
00071   c3->cd(3);
00072   hRatioFin3->Draw();
00073   c3->cd(4);
00074   hRatioFin4->Draw();
00075   c3->cd(5);
00076   hRatioFin5->Draw();
00077 #endif
00078 
00079   //make projections to get ratio of MC/Data for each eta bin
00080   TH1 *hRatioAll[10];
00081   hRatio->ProjectionY("hRatioA",1,1);   hRatioAll[0]=hRatioA;
00082   hRatio->ProjectionY("hRatioB",2,2);   hRatioAll[1]=hRatioB;
00083   hRatio->ProjectionY("hRatioC",3,3);   hRatioAll[2]=hRatioC;
00084   hRatio->ProjectionY("hRatioH",8,8);   hRatioAll[7]=hRatioH;
00085   hRatio->ProjectionY("hRatioI",9,9);   hRatioAll[8]=hRatioI;
00086   hRatio->ProjectionY("hRatioJ",10,10); hRatioAll[9]=hRatioJ;
00087 
00088   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
00089   //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
00090   if(print) c=new TCanvas("aa","aa",800,600);
00091   for(int j=0; j<10; j++){
00092     if(j>2 && j<7) continue;
00093     hRatioAll[j]->SetTitle(Form("Eta bin %c: Ratio of tracks QCD MC to run 9 data (track pt > 5)",'A'+j));
00094     hRatioAll[j]->SetMinimum(0); hRatioAll[j]->SetLineWidth(2);
00095     Lx=hRatioAll[j]->GetListOfFunctions(); hRatioAll[j]->SetStats(false);
00096     ln=new TLine(-PI,normal[j],PI,normal[j]);
00097     ln->SetLineColor(kRed); ln->SetLineWidth(2);  
00098     Lx->Add(ln);
00099     if(print){
00100       hRatioAll[j]->Draw();
00101       c->Print(Form("etabin%c.png",'A'+j));
00102     }
00103   }
00104   
00105   
00106   float phiRad=999;
00107   float etaDet=999;
00108   
00109   //initialize and clear eff values
00110   float effic[10][24];
00111   float efficSec[10][24];
00112   for(int i=0;i<10;i++){ 
00113     for(int l=0;l<24;l++){
00114       effic[i][l]=0;
00115       efficSec[i][l]=0;
00116     }
00117   }
00118   
00119   //calculate inefficiency
00120   for(int ieta=0; ieta<10; ieta++){
00121     if(ieta>2 && ieta<7) continue;
00122     for(int iphi=0; iphi<24; iphi++){
00123       float binVal=hRatioAll[ieta]->GetBinContent(iphi+1);
00124       effic[ieta][iphi]=binVal/normal[ieta];
00125       //cout<<"Eta "<<Form("%c",ieta+'A')<<" phi "<<iphi<<" effic "<<binVal<<endl;
00126       phiRad=hRatioAll[ieta]->GetXaxis()->GetBinCenter(iphi+1);
00127       
00128       //find tpc sector
00129       if(ieta<3) etaDet=-1;
00130       else if(ieta>6) etaDet=1;
00131       const float PI=TMath::Pi();
00132       int sec=0;
00133       float phi=phiRad/PI*180; // now in degrees
00134       if (etaDet>0) { // West TPC
00135         float x=75-phi;
00136         while(x<0) x+=360;
00137         sec=1+(int)( x/30.);
00138       } else {
00139         float x=phi-105;
00140         while(x<0) x+=360;
00141         sec=13+(int)( x/30.);
00142       }
00143       
00144       efficSec[ieta][sec-1]+=effic[ieta][iphi]/2; //average 2 bins in sector
00145     }
00146 #if 0 //print weight arrays to use on data in W algo
00147     cout<<"mWeight"<<ieta+1<<"={";
00148     for(int isec=0; isec<24; isec++){
00149       //cout<<"sector "<<isec+1<<" eta "<<ieta<<" effic "<<efficSec[ieta][isec]<<endl;
00150       cout<<efficSec[ieta][isec]<<",";
00151     }
00152     cout<<"};"<<endl;
00153 #endif
00154     //set weight arrays to use on MC in W algo
00155     for(int isec=0; isec<24; isec++){
00156       if(efficSec[ieta][isec]!=0) efficSec[ieta][isec]=1/efficSec[ieta][isec];
00157     //  cout<<"sector "<<isec+1<<" eta "<<ieta<<" effic "<<efficSec[ieta][isec]<<endl;
00158     }
00159 
00160 
00161   }
00162 
00163   //#endif
00164   
00165   
00166   //Data corrected for efficiencies in outer eta bins
00167   //hCorr->Draw("colz");
00168   hCorr->Rebin2D(10,10);
00169   hRatioCorr=(TH2F*) hQcdRebin->Clone();
00170   hRatioCorr->Divide(hCorrRebin);
00171   //hRatioCorr->Draw("colz");
00172 
00173   //add sector across pi boundry separately
00174   hCorr->ProjectionX("h1",1,1);
00175   hCorr->ProjectionX("h24",24,24);
00176   TH1 * hMidCorr[12];
00177   hMidCorr[11]=h1;hMidCorr[11]->Add(h24);
00178   
00179   bool print2=false;
00180   if(print2) c2=new TCanvas("bb","bb",800,600);
00181   for(int iphi=1; iphi<22; iphi+=2){ // get sector slices
00182     hCorr->ProjectionX(Form("h",iphi+1),iphi+1,iphi+2);   
00183     hMidCorr[(iphi-1)/2]=h;
00184    
00185     //find tpc sector
00186     float phiRad=hCorr->GetYaxis()->GetBinCenter(iphi+1);
00187     float etaDet=1;
00188     const float PI=TMath::Pi();
00189     int sec1=0; int sec2=0;
00190     float phi=phiRad/PI*180; // now in degrees
00191     float x=75-phi;
00192     while(x<0) x+=360;
00193     sec1=1+(int)( x/30.);
00194     float x=phi-105;
00195     while(x<0) x+=360;
00196     sec2=13+(int)( x/30.);
00197 
00198     
00199     //#if 0
00200     //calculate efficiency at midrapidity
00201     float etaVal[10]; float etaEffic[10];
00202     for(int ieta=0; ieta<10; ieta++)
00203       etaVal[ieta]=hMidCorr[(iphi-1)/2]->GetBinContent(ieta+1);
00204     float midRapVal=(etaVal[2]+etaVal[7])/2;
00205     if(sec2==13) float midRapVal=etaVal[2];
00206     for(int ieta2=3; ieta2<7; ieta2++){
00207       etaEffic[ieta2]=etaVal[ieta2]/midRapVal;
00208       if(ieta2<5){
00209         sec=sec2;
00210         //cout<<"sec "<<sec2<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
00211       }
00212       else{
00213         sec=sec1;
00214         //cout<<"sec "<<sec1<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
00215       }
00216       efficSec[ieta2][sec-1]=etaEffic[ieta2];
00217     }
00218     //#endif
00219 
00220     if(print2){
00221       Lx=h->GetListOfFunctions();
00222       h->SetStats(false);
00223       ln=new TLine(-0.44,midRapVal,0.44,midRapVal);
00224       ln->SetLineColor(kRed);  ln->SetLineWidth(2);
00225       Lx->Add(ln);
00226       h->SetMinimum(0); h->SetFillColor(kBlue);
00227       h->Draw();
00228       h->SetTitle(Form("Sectors %d and %d Tracks vs Eta",sec1,sec2)); 
00229       //c2->Print(Form("sec%d_%d.png",sec1,sec2));
00230     }
00231   }
00232 
00233   //do sector 9,15 separately
00234   for(int ieta=0; ieta<10; ieta++)
00235     etaVal[ieta]=hMidCorr[11]->GetBinContent(ieta+1);
00236   float midRapVal=(etaVal[2]+etaVal[7])/2;
00237   if(print2) {
00238     hMidCorr[11]->SetTitle("Sectors 9 and 15 Tracks vs Eta"); 
00239     hMidCorr[11]->Draw(); hMidCorr[11]->SetMinimum(0); c2->Print("sec9_15.ps");
00240     Lx=h->GetListOfFunctions();
00241     ln=new TLine(-PI,midRapVal,PI,midRapVal);
00242     ln->SetLineColor(kRed);  
00243     Lx->Add(ln);
00244   }
00245   for(int ieta2=3; ieta2<7; ieta2++){
00246     etaEffic[ieta2]=etaVal[ieta2]/midRapVal;
00247     if(ieta2<5){
00248       sec=15;
00249       //cout<<"sec "<<sec2<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
00250     }
00251     else{
00252       sec=9;
00253       //cout<<"sec "<<sec1<<" iphi "<<iphi<<" ieta "<<ieta2<<" effic "<<etaEffic[ieta2]<<endl;
00254     }
00255     efficSec[ieta2][sec-1]=etaEffic[ieta2];
00256   }
00257 
00258 #if 0  
00259   //print effic numbers
00260   for(int jeta=0;jeta<10;jeta++){
00261     cout<<"float mWeight"<<jeta+1<<"={";
00262     for(int jsec=0;jsec<24;jsec++){
00263       //cout<<"sec "<<jsec+1<<" eta "<<jeta<<" effic "<<efficSec[jeta][jsec]<<endl;
00264       cout<<efficSec[jeta][jsec]<<",";
00265     }
00266     cout<<"};"<<endl;
00267   }
00268 #endif
00269   
00270 
00271 }