StRoot: StHighptPool/macros/yieldVvtxZ.C Source File

00001 #include "commonmacro/common.h"
00002 #include "common/Name.cc"
00003 #include "commonmacro/histutil.h"
00004 
00005 void yieldVvtxZ(const char* inName=
00006                 "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
00007                 const char* psDir="ps",
00008                 int cut = 1,
00009                 const char* outDir="./",
00010                 const char* more = "embeddingfile",
00011                    float extraValue = 0) 
00012 {
00013  
00014   gSystem->Load("StHiMicroAnalysis");
00015  
00016   Cut::SetCut(cut);
00017   Cut::ShowCuts();
00018 
00019   TFile* inRoot = new TFile(inName);
00020   TFile* embedRoot;
00021   if(more){
00022     strcpy(name,more);
00023   }
00024   else if(extraValue>0){
00025     // attempt to find the embed file
00026     char* trigger = (strstr(inName,"central")) ? "central" :"minbias";
00027     sprintf(name,"~/wrk/assoc/links/P01hi.%s.HighpT_piminus_%d.hist/efficiency_cut%d.hist.root",trigger,(int)extraValue,cut);
00028   }
00029   else{
00030     cout << "Unknown embed file" << endl; return;
00031   }
00032 
00033   TString half = "";
00034   if(strstr(inName,"east")){ half="east"; }
00035   if(strstr(inName,"west")){ half="west"; }
00036 
00037   embedRoot = new TFile(name);
00038   cout << "--------------------------" << endl;
00039   cout << "in name=" << inName << endl
00040        << "ps dir=" << psDir << endl
00041        << "cut=" << cut << endl
00042        << "embed name=" << name << endl
00043        << "half = " << half << endl;
00044   cout << "--------------------------" << endl;
00045 
00046   if(!inRoot){
00047     cout << "cannot find the infile" << endl;
00048     return;
00049   }
00050   if(!embedRoot){
00051     cout << "cannot find " << name << endl;
00052   }
00053   float minRealPt=2,maxRealPt=3;
00054   float minMcPt=2,maxMcPt=8;
00055 
00056   float minVtxZ=-160,maxVtxZ=160;
00057   float etaAry[][2] = { { .0,.1}, {.2,.3}, {.4,.5}, {.6,.7},
00058     {-.1,0}, {-.3,-.2},{-.5,-.4},{-.7,-.6} };
00059   const int nEta = 8; int nx=2; int ny=4;
00060   int nRebin=2;
00061   int nMcRebin=4;
00062   
00063 
00064   TH3* h3;
00065   TH3* h3Real; // 10 cm bins, 
00066   TH3* h3Mc[2]; // 5 cm bins
00067 
00068   TH2* h2Real,*h2Mc[2];
00069   TH1* h1Real[2]; // raw and corrected
00070   TH1* h1Mc[2][nEta]; // raw,matched
00071   TH1* h1VtxZ,*h1;
00072 
00073   // real stuff
00074   h3Real=(TH3*)inRoot.Get("Minus.h3VtxZEtaPrPtPr");
00075   h3 = (TH3*)inRoot.Get("Plus.h3VtxZEtaPrPtPr");
00076   h3Real->Add(h3);
00077   //slice off pt
00078   h2Real=HistSlice(h3Real,"","",0,minRealPt,maxRealPt,"yx","e");
00079 
00080   // mc stuff
00081   h3Mc[0]=(TH3*)embedRoot.Get("h3McRawVtxZEtaPt");
00082   h3Mc[1]=(TH3*)embedRoot.Get("h3MatchedVtxZEtaPt");
00083 
00084   h2Mc[0]=HistSlice(h3Mc[0],"","",0,minMcPt,maxMcPt,"yx","e");
00085   h2Mc[1]=HistSlice(h3Mc[1],"","",0,minMcPt,maxMcPt,"yx","e");
00086 
00087   h3=(TH3*)inRoot.Get("h3VtxXYZ");
00088   h1VtxZ=(TH3*)h3->Project3D("ze"); h1VtxZ->Rebin(nMcRebin);
00089 
00090   char* type[] = {"real","mc"}; 
00091   int marker[] = {4,8};
00092   char* opt[] = { "p","psame"};
00093   char* opte[] = {"e","esame"};
00094   float markerSize = 0.4;
00095   
00096 
00097   TCanvas c1("c1","c1",400,500);
00098   TText* t=new TText;
00099   gStyle->SetOptStat(0);
00100 
00101   Cut::SetCut(cut);
00102   
00103   //-------------------------------------------------------------
00104   // see what the eta slices look like
00105   TH1* hSlice = new TH1D("hSlice","hSlice",20,-200,200);
00106   SetRange(hSlice->GetXaxis(),minVtxZ,maxVtxZ);
00107   hSlice->SetMinimum(0); hSlice->SetMaximum(200);
00108 
00109   sprintf(title,"eta range(cut %d)",cut);
00110   TLine line;
00111   Divide(&c1,nx,ny,title,inName);
00112   for(int iEta=0;iEta<nEta;iEta++){
00113     c1.cd(iEta+1);
00114 
00115     sprintf(title,"%.1f<eta<%.1f",etaAry[iEta][0],etaAry[iEta][1]);
00116     h1=(TH1*)hSlice->Clone();
00117     h1->SetTitle(title);
00118     h1->Draw(); TAxis* axis=h1->GetXaxis();    
00119 
00120     for(int i=axis->GetFirst();i<=axis->GetLast();i++){
00121       float dip1 = 3.14159/2.-2.*atan(TMath::Exp(-etaAry[iEta][0]));
00122       float dip2 = 3.14159/2.-2.*atan(TMath::Exp(-etaAry[iEta][1]));
00123       
00124       //cout << "dip1=" << dip1*180./3.14159 
00125       //           << ",dip2="<<dip2*180./3.14159 << endl;
00126 
00127       cout << "eta="<< etaAry[iEta][1] << ",z=" << 192*tan(dip2) << endl;
00128 
00129       float xmax1 = 200*tan(dip2) + axis->GetBinLowEdge(i);
00130       float xmax2 = 200*tan(dip1) + axis->GetBinLowEdge(i);
00131 
00132       float xmin = axis->GetBinLowEdge(i);
00133       line.SetLineStyle(2);
00134       line.SetLineColor(kBlue);
00135       line.DrawLine(xmin,0,xmax1,200);
00136       line.SetLineStyle(1);
00137       line.SetLineColor(kRed);
00138       line.DrawLine(xmin,0,xmax2,200);
00139     }
00140     line.SetLineStyle(1);
00141     line.SetLineColor(kBlack);
00142     line.DrawLine(0,0,0,200);
00143 
00144     line.DrawLine(minVtxZ,60,maxVtxZ,60);
00145 
00146   }
00147   sprintf(name,"etaRange");
00148   Print(&c1,psDir,name);
00149 
00150   //  return;
00151   // look at raw and matched counts
00152 
00153   sprintf(title,"raw and matched yield V vtxZ (eta slices) (cut %d)",cut);
00154   Divide(&c1,nx,ny,title,inName);
00155   for(int iEta=0;iEta<nEta;iEta++){
00156     c1.cd(iEta+1);
00157     for(int i=0; i<2; i++){
00158       h1Mc[i][iEta]=HistSlice(h2Mc[i],"",h2Mc[i]->GetTitle(),0,
00159                               etaAry[iEta][0],etaAry[iEta][1],"x","e");
00160       h1Mc[i][iEta]->SetMarkerStyle(marker[i]);
00161       if(nMcRebin>1)h1Mc[i][iEta]->Rebin(nMcRebin);
00162       if(i==0)SetMinMax(h1Mc[i][iEta],h1Mc[i][iEta]->GetMinimum()*0.2,
00163                         h1Mc[i][iEta]->GetMaximum()*1.2);
00164       SetRange(h1Mc[i][iEta]->GetXaxis(),minVtxZ,maxVtxZ);
00165       h1Mc[i][iEta]->SetMarkerSize(markerSize);
00166       h1Mc[i][iEta]->Draw(opt[i]);
00167     }
00168   }
00169   sprintf(name,"mcYieldVvtxZEtaSlices");
00170   Print(&c1,psDir,name);
00171 
00172   // look at efficiency*acceptance
00173   sprintf(title,"efficiency*acceptance V vtxZ (eta slices) (cut %d)",cut);
00174   Divide(&c1,nx,ny,title,inName);
00175   for(int iEta=0;iEta<nEta;iEta++){
00176     c1.cd(iEta+1); gPad->SetGridx(); gPad->SetGridy();
00177     h1Mc[1][iEta]->Divide(h1Mc[0][iEta]); SetMinMax(h1Mc[1][iEta],0,1);
00178     h1Mc[1][iEta]->SetXTitle("vtxZ");
00179     h1Mc[1][iEta]->Draw("e"); 
00180    
00181   }
00182   sprintf(name,"efficiencyVvtxZEtaSlices");
00183   Print(&c1,psDir,name);
00184 
00185   // look at raw and uncorrected real yields V vtx Z
00186   sprintf(title,"raw and corrected yield/Nevent V vtxZ (eta slices) (cut %d)",cut);
00187   Divide(&c1,nx,ny,title,inName);
00188   for(int iEta=0;iEta<nEta;iEta++){
00189     c1.cd(iEta+1); gPad->SetGridx(); gPad->SetGridy();
00190     h1Real[0]=HistSlice(h2Real,"",h2Real->GetTitle(),0,
00191                         etaAry[iEta][0],etaAry[iEta][1],"x","e");
00192     h1Real[0]->Rebin(nRebin);
00193     h1Real[0]->Divide(h1VtxZ); // uncorrected
00194     h1Real[1]=(TH1*)h1Real[0]->Clone(); // corrected
00195     h1Real[1]->Divide(h1Mc[1][iEta]);
00196     
00197 
00198     for(int i=0; i<2; i++){
00199       h1Real[i]->SetMarkerStyle(marker[i]);
00200       if(i==0)SetMinMax(h1Real[i],0,.5);
00201       SetRange(h1Mc[i][iEta]->GetXaxis(),minVtxZ,maxVtxZ);
00202       h1Real[i]->SetMarkerSize(markerSize);
00203       h1Real[i]->SetXTitle("vtxZ");
00204       h1Real[i]->Draw(opte[i]);
00205     }
00206   }
00207   Print(&c1,psDir,"yieldVvtxZEtaSlices");
00208   
00209 }
00210 
00211 
00212 
00213 
00214   
00215   
00216