StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
proLcpTree.C
1 // projects TTree in to histos using cuts (if any)
2 // former rdLcpTree()
3 
4 proLcpTree( int cut=0, TString run="R3010006", int maxEve=10000,
5  TString wrkDir="/star/data04/sim/balewski/LcpRun2/maxEta1.0/",
6  TString outPath="fixMe/"
7  ){
8  printf("#cut %d -->%s \n",cut,outPath.Data(),run.Data());
9 
10  gStyle->SetPalette(1,0);
11  gStyle->SetOptStat(1111111);
12 
13  printf("#run %s\n",run.Data());
14 
15 
16  TChain *chain = new TChain("T-LCP","dummName");
17 
18  TString item=wrkDir+"/tree/"+run;
19  item=item+".tree.root";
20  printf("#tree %s\n", item.Data());
21  int ret=chain->AddFile(item,-1);
22  printf("chain ret %d \n",ret);
23 
24 
25  // ============================ DEFINE OUTPUT HISTO
26  const float Pi=3.1416;
27  TString outF=wrkDir+outPath+run+".hist.root";
28 
29  hFile=new TFile(outF,"RECREATE");
30  assert(hFile->IsOpen());
31 
32  int i;
33  TH1F *h1[8];
34  TH2F *h2L[100]; int nh2L=0;
35 
36  //........................... QA histo
37  h1[0]=new TH1F("bx","rates vs. bXing",129,-0.5,128.5);
38  h1[1]=new TH1F("pt","Leading Particle pT (GeV/c)",100,0.,10.);
39  //h1[1]=new TH1F("pt","Leading Particle pT (GeV/c)",135,0.025,6.775); // Olga's definition
40  h1[2]=new TH1F("Zv","Z vertex/cm",100,-250.,250.);
41  h1[3]=new TH1F("nPrim","No.of valid prim TPC tracks ",50,-0.5,49.5);
42  h1[4]=new TH1F("cosm", "Mike's Cosmics Rejector return value",5,-0.5,4.5);
43  h1[5]=new TH1F("phi","LCP phi/rad",100,-3.15.,3.15);
44  h1[6]=new TH1F("eta","LCP Eta",100,2.,2.);
45  h1[7]=new TH1F("q","LCP charge",10,-2.,2.);
46  float bxF=-0.5, bxL=119.5;
47  int nbx=120;
48  int nphi=16;
49  TH2F *h2= new TH2F("PhiBxAll","LCP Phi/rad vs. bXing[0-119]",nbx,bxF,bxL,nphi,-Pi,Pi);
50 
51  /* nomanclautre for LCP cuts:
52  pTL,M,H for pt in[0.4,0.7,1.0,2.0]
53  pT1,2,3,4,5,6 for pt in[0.4,1,2,3,4,5,6]
54  etaBB,Bc,Fc,FF for eta in [low,-0.5,0.,0.5,high]
55  q+,- for charge + or -
56  */
57 
58  //.................. histo for PT 1,2,3,4,
59  const int mxPt=6;
60  TH2F *h2pT[mxPt];
61  for(i=0;i<mxPt;i++) {
62  char tit1[100], tit2[1000];
63  sprintf(tit1,"PhiBxPt%d",i+1);
64  sprintf(tit2,"LCP Phi/rad vs. bXing[0-119], cut=Pt%d",i+1);
65  h2pT[i]= new TH2F(tit1,tit2,nbx,bxF,bxL,nphi,-Pi,Pi);
66  h2L[nh2L++]=h2pT[i];
67  //printf("-%s-%s-\n",tit1,tit2);
68  }
69 
70  //.................. histo for PT L,M,H
71  const int mxRt=3;
72  TH2F *h2rT[mxRt];
73  char coreRt[mxRt]={'L','M','H'};
74  for(i=0;i<mxRt;i++) {
75  char tit1[100], tit2[1000];
76  sprintf(tit1,"PhiBxPt%c",coreRt[i]);
77  sprintf(tit2,"LCP Phi/rad vs. bXing[0-119], cut=Pt%c",coreRt[i]);
78  h2rT[i]= new TH2F(tit1,tit2,nbx,bxF,bxL,nphi,-Pi,Pi);
79  h2L[nh2L++]=h2rT[i];
80  }
81 
82  //.................. histo for Eta selection
83  const int mxEta=4;
84  TH2F *h2eta[mxEta];
85  char *coreEta[mxEta]={"BB","Bc","Fc","FF"};
86  for(i=0;i<mxEta;i++) {
87  char tit1[100], tit2[1000];
88  sprintf(tit1,"PhiBxEta%s",coreEta[i]);
89  sprintf(tit2,"LCP Phi/rad vs. bXing[0-119], cut=Eta%s",coreEta[i]);
90  h2eta[i]= new TH2F(tit1,tit2,nbx,bxF,bxL,nphi,-Pi,Pi);
91  h2L[nh2L++]=h2eta[i];
92  }
93 
94  //.................. histo for charge selection
95  const int mxq=2;
96  TH2F *h2q[mxq];
97  char *coreq[mxq]={"neg","pos"};
98  for(i=0;i<mxq;i++) {
99  char tit1[100], tit2[1000];
100  sprintf(tit1,"PhiBxQ%s",coreq[i]);
101  sprintf(tit2,"LCP Phi/rad vs. bXing[0-119], cut=Q%s",coreq[i]);
102  h2q[i]= new TH2F(tit1,tit2,nbx,bxF,bxL,nphi,-Pi,Pi);
103  h2L[nh2L++]=h2q[i];
104  }
105 
106 
107  //========================== HISTO CREATED
108 
109  int N=chain->GetEntries();
110  if(maxEve<=0) maxEve=N;
111 
112  printf("#sort %d of %d\n",maxEve,N);
113  // hFile->ls();
114 
115  // associate varibales with the tree
116  float vz,pt,phi,eta;
117  int nPrim,bx,cosm,id,nFit,charge;
118  chain->SetBranchAddress("bx120",&bx);
119  chain->SetBranchAddress("vz",&vz);
120  chain->SetBranchAddress("id",&id);
121  chain->SetBranchAddress("pt",&pt);
122  chain->SetBranchAddress("phi",&phi);
123  chain->SetBranchAddress("eta",&eta);
124  chain->SetBranchAddress("q",&charge);
125  chain->SetBranchAddress("nPrim",&nPrim);
126  chain->SetBranchAddress("nFit",&nFit);
127  chain->SetBranchAddress("cosm",&cosm);
128 
129  system("date");
130  printf("\n::::::::::::::: L O O P O V E R E V E N T S :::::\n");
131  int k;
132  int nLcp=0;
133  for(k=0;k<N;k++) {
134  if(k>=maxEve) break;
135  int ret=chain->GetEntry(k);
136  assert(ret);
137  if(k%1000==0)printf("%d %d %d %f id=%d nFit=%d\n",k,ret,nPrim ,vz,id,nFit);
138  // new final cuts
139  if(nPrim<3 ) continue;
140  if( pt>5. ) continue;
141 
142  if(cut==1 && (nPrim<5 || nPrim>20) ) continue;
143  if(cut==2 && fabs(vz) >50) continue;
144  if(cut==3 && (pt<1. || pt>3.) ) continue;
145  if(cut==4 && charge<0. ) continue;
146 
147  h1[3]->Fill(nPrim);
148  h1[0]->Fill(bx);
149  h1[2]->Fill(vz);
150  h1[4]->Fill(cosm);
151  if(nFit<=0) continue; // LCP must exist
152  nLcp++;
153  h1[1]->Fill(pt);
154  h1[5]->Fill(phi);
155  h1[6]->Fill(eta);
156  h1[7]->Fill(charge);
157  h2->Fill(bx,phi);
158 
159  // .............. cuts on LCP
160  //... Pt case
161  int i=(int)pt;
162  if(i>=0 && i<mxPt) h2pT[i]->Fill(bx,phi);
163  i=-1;
164  if(pt<0.7) i=0;
165  else if(pt<1.0) i=1;
166  else if(pt<3.0) i=2;
167  if(i>=0) h2rT[i]->Fill(bx,phi);
168 
169  //... Eta case
170  if(eta<-0.5) i=0;
171  else if(eta<0.) i=1;
172  else if(eta<0.5) i=2;
173  else i=3;
174  h2eta[i]->Fill(bx,phi);
175 
176  //... charge case
177  i=0;
178  if(charge>0) i=1;
179  h2q[i]->Fill(bx,phi);
180 
181  }
182  system("date");
183  printf("nLcp=%d\n",nLcp);
184 
185  hFile->Write();// Save all objects in this file
186 
187  //................... print out 2D-histo
188 
189  TAxis *X=h2->GetXaxis();
190  TAxis *Y=h2->GetYaxis();
191  int nx=X->GetNbins();
192  int ny=Y->GetNbins();
193  printf("#title %s\n", h2->GetTitle());
194  printf("#X-axis %f, %f %d %s \n", X->GetXmin(), X->GetXmax(),nx);
195  printf("#Y-axis %f, %f %d %s \n", Y->GetXmin(), Y->GetXmax(),ny);
196  printf("#Integral(valid)= %.1f \n",h2->Integral());
197 
198  for(i=0;i<nh2L;i++) {
199  TH2F *hx= h2L[i];
200  printf("#Integral(%s)= %.1f \n",hx->GetName(),hx->Integral());
201  }
202 
203  c=new TCanvas();
204  c->Divide(2,3);
205  c->cd(1); h1[1]->Draw(); gPad->SetLogy();
206  c->cd(2); h1[5]->Draw();
207  c->cd(3); h1[6]->Draw();
208  c->cd(4); h2->Draw("colz");
209  c->cd(5); h1[2]->Draw();
210  c->cd(6); h1[7]->Draw();
211 
212  TString outPS=outF.ReplaceAll("hist.root","pro.ps");
213 
214  c->Print(outPS);
215 
216 
217 }
Generated by   doxygen 1.8.5