fitPts.C

#include "commonmacro/common.h"
#include "common/Name.cc"
#include "commonmacro/histutil.h"

void fitPts(const char* inName=
            "links/P01hi.minbias.2000.hist/hianalysis_1000.hist.root",
            const char* psDir="ps",
            int cut = 1,
            const char* outDir="./",
            const char* more = "west",
            float extraValue = 0)
{
  cout << "--------------------------" << endl;
  cout << "in name=" << inName << endl
       << "ps dir=" << psDir << endl
       << "cut=" << cut << endl;
  cout << "--------------------------" << endl;
 
  inRoot = new TFile(inName);

  gSystem->Load("StHiMicroEvent");
  gSystem->Load("StHiMicroAnalysis");
  Cut::SetCut(cut);
  Cut::ShowCuts();

  if(!inRoot){
    cout << "cannot find the infile" << endl;
    return;
  }
  TText* t = new TText;
  float ptAry[] = {2,3,4,5,6};
  char* ptType[] = { "PtPr","PtGl"};
  char* chargeType[]= { "Pos","Neg","PosPlusNeg"};
  TH1 *ha, *hb;
  TCanvas c1("c1","c1",400,500);
  TH3* h3[2][3]; // ptType; plus,minus,plus+minus
  TH2* h2a, *h2b; TH1* h1a, *h1b;
  float ptmin=1.5,ptmax=4.0;
  char cname[100],ctitle[100];
  int nc=2;
  int nptType=1;

  int fitPtsCut=Cut::mFitPts[0];
  float fitptsary[] = { 3,fitPtsCut};
  char* fitptstitle[] = { "Wide","Cut"};

  gStyle->SetOptStat(0);
  //-------------------------------------------------------------------
  // 
  // plus,minus,plus+minus
  
  for(int iptType=0;iptType<nptType;iptType++){
    //plus
    sprintf(name,"Minus.h3PhiPrFitPts%s",ptType[iptType]);
    cout << name << endl;
    h3[iptType][0]=(TH3*)inRoot.Get(name);
    // minus
    sprintf(name,"Plus.h3PhiPrFitPts%s",ptType[iptType]);
    h3[iptType][1]=(TH3*)inRoot.Get(name);
    // plus+minus
    h3[iptType][2]=(TH3*)h3[iptType][0]->Clone();
    sprintf(name,"PlusPlusMinus.h3PhiPrFitPts");
    h3[iptType][2]->SetName(name); h3[iptType][2]->SetTitle(name);
    h3[iptType][2]->Add(h3[iptType][1]);
  }
  //---------------------------------------------------------------
  for(int if=0;if<2; if++){
    for(int iptType=0;iptType<nptType;iptType++){

      //**** fit pts, Pos and Neg on same plot
      sprintf(cname,"fitPts%s%s",ptType[iptType],
              fitptstitle[if]);
      sprintf(ctitle,"fit pts, %s slices  %d-45 (cut %d)",
              ptType[iptType],fitptsary[if],cut);
      Divide(&c1,2,2,ctitle,inName);
    
      for(int j=0;j<4;j++){
        c1.cd(j+1);
        if(j==0){
          TLegend l(0.1,0.4,0.3,0.5); l.SetBorderSize(0);
          l.SetFillColor(0);
        }
        for(int ic=0;ic<nc;ic++){
          h2a=(TH2*)HistSlice(h3[iptType][ic],"","",0,
                              ptAry[j],ptAry[j+1],"xy");
          h1a=HistSlice(h2a,"",h2a->GetTitle(),0,0,0,"x");
          SetRange(h1a->GetXaxis(),fitptsary[if],46);
          h1a->SetMarkerSize(0.6);
          if(ic<1){
            SetMinMax(h1a,0,h1a->GetMaximum()*2.5);
            h1a->SetMarkerStyle(4); 
            h1a->Draw("p");
            PrintMeanRms(h1a,0.1,.8);
            
          }
          else if (ic==1){
            h1a->SetMarkerStyle(8);
            h1a->Draw("psame");
            PrintMeanRms(h1a,0.1,.7);
          }     
          else if (ic==2){
            h1a->SetMarkerStyle(12);
            h1a->Draw("psame");
            PrintMeanRms(h1a,0.1,.6);
          }
          if(j==0) l.AddEntry(h1a,chargeType[ic],"p");
        }
        if(j==0) l.Draw();
      } // pt loop
      Print(&c1,psDir,cname);

      //*** fit pts distribution, phi slices
      // combine pt 1.5-4
      // one plot for each charge type
      
      for(int ic=0;ic<nc;ic++){
        sprintf(cname,"fitPts%s%sPhiSlices%s",
                ptType[iptType],chargeType[ic],fitptstitle[if]);
        sprintf(ctitle,"fit pts %d-45, phi slices (%.1f<%s<%.1f) %s (cut %d)",
                fitptsary[if],ptmin,ptType[iptType],ptmax,chargeType[ic],cut);
        Divide(&c1,3,4,ctitle,inName);
        
        h2a=(TH2*)HistSlice(h3[iptType][ic],"","",0,
                            ptmin,ptmax,"xy");
        for(int iPhi=1; iPhi<=h2a->GetNbinsY(); iPhi++){
          c1.cd(iPhi);
          TAxis *axis = h2a->GetYaxis();
          sprintf(title,"%d<phi<%d",
                  axis->GetBinLowEdge(iPhi),axis->GetBinUpEdge(iPhi));
          sprintf(name,"%s%s%s%s",
                  chargeType[ic],ptType[iptType],fitptstitle[if],title);
          
          h1a = h2a->ProjectionX(name,iPhi,iPhi);
          h1a->SetMarkerStyle(8); h1a->SetMarkerSize(0.6);
          h1a->SetTitle(title); SetRange(h1a->GetXaxis(),fitptsary[if],45.5);
          h1a->Draw("p");
          PrintMeanRms(h1a,0.1,0.6,0.08);
        }
        Print(&c1,psDir,cname);

      } // ic
      
      //** mean fit, profile vs phi
      for(int ic=0;ic<nc;ic++){
        sprintf(cname,"fitPtsMeanVPhi%s%s%s",
                ptType[iptType],chargeType[ic],fitptstitle[if]);
        sprintf(ctitle,"mean fit pts V phi %d-45 (%.1f<%s<%.1f) %s (cut %d)",
                fitptsary[if],ptmin,ptType[iptType],ptmax,chargeType[ic],cut);
        Divide(&c1,1,1,ctitle,inName);
        c1.cd(1); gPad->SetGridx(); gPad->SetGridy();
        h2a=(TH2*)HistSlice(h3[iptType][ic],"","",0,
                            ptmin,ptmax,"xy");
        int low=h2a->GetXaxis()->FindBin(fitptsary[if]);
        sprintf(name,"meanVphi%s%s%s",
                chargeType[ic],ptType[iptType],fitptstitle[if]);
        TProfile *p=h2a->ProfileY(name,low,9999); p->SetTitle("");
        SetMinMax(p,30,40);
        p->SetMarkerStyle(8);   
        p->Draw();     

        Print(&c1,psDir,cname);
        
      } // ic
      
    }
  }
}