StRoot: StEEmcPool/macros/LCP/addRuns.C Source File

00001 const int mxBx=120, mxPol=4;
00002 
00003 char *cpolBY[mxPol]={"B+Y+","B+Y-","B-Y+","B-Y-"};
00004 char bXselection[mxBx];
00005 int polPattBY[mxBx];
00006 double Lum[mxPol],Lum2[mxPol];
00007 TFile  *outH=0;
00008 TString cutL=" ";
00009 
00010 
00011 //-------------------------------------
00012 //-------------------------------------
00013 void  updatePattern(char *run) {
00014   int irun=atoi(run+1);
00015   assert(irun>0);
00016   printf("Update pattern for run=%d\n",irun);
00017 
00018   char * polPattY,* polPattB,*dum;
00019   
00020   char *  //yellow
00021     selPatt="a.........k....BBBBBa.........k....pppppa.........k....YYYYY";
00022   dum      ="|         |         |         |         |         |        |";
00023   dum      ="1         11        21        31        41        51       60";   
00024 
00025   if(irun<3017036) { // (STAR) Pattern1  
00026     polPattY= "--++--++--++--++--++--++--++--++--++--++--++--++--++--+aaaaa";
00027     polPattB= "-+-+-+-+-+-+-+-aaaaa-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+";
00028     dum     = "|         |         |         |         |         |        |";
00029     dum     = "1         11        21        31        41        51       60";   
00030   } else { // (STAR) Pattern2  with 3 unpolarized bXings
00031     polPattY= "0-+-+-+-+-+-+-+-+-+-0+-+-+-+-+-+-+-+-+-+0-+-+-+-+-+-+-+aaaaa";
00032     polPattB= "0++--++--++--++aaaaa0--++--++--++--++--+0+--++--++--++--++--";
00033     dum     = "|         |         |         |         |         |        |";
00034     dum     = "1         11        21        31        41        51       60";   
00035   }
00036 
00037   int i;
00038   for(i=0;i<60;i++) {
00039     bXselection[2*i]=selPatt[i];
00040     bXselection[2*i+1]='i'; // ignore
00041   }
00042   
00043   for(i=0;i<60;i++) {
00044     int k=2*i;
00045     polPattBY[k]=-2; //default
00046     if     (polPattB[i]=='+' && polPattY[i]=='+' ) polPattBY[k]=0;
00047     else if(polPattB[i]=='+' && polPattY[i]=='-' ) polPattBY[k]=1;
00048     else if(polPattB[i]=='-' && polPattY[i]=='+' ) polPattBY[k]=2;
00049     else if(polPattB[i]=='-' && polPattY[i]=='-' ) polPattBY[k]=3;
00050   }
00051 
00052 }
00053 
00054 
00055 
00056 //-------------------------------------
00057 //        M A I N 
00058 //-------------------------------------
00059 addRuns(TString fill=0, TString runL=0, TString wrkDir=0) {
00060   gStyle->SetPalette(1,0);
00061   gStyle->SetOptStat(1111111);
00062   
00063   //fill="F2201x";  runL="R2362128 R777";
00064   //wrkDir="/star/data04/sim/balewski/LcpRun2/setEta10/";
00065 
00066   // create output histo
00067   TString outF=wrkDir+fill+".hist.root";
00068   //TString outF=fill+".hist.root";
00069   outH=new TFile(outF,"RECREATE");
00070   assert(outH->IsOpen());
00071   printf("save outH -->%s\n", outF.Data());
00072 
00073 
00074   // ----- access scaler data
00075   TString inpScal="JoannaScal/scalVsRun.hist.root";
00076   scal=new TFile(inpScal);
00077   assert(scal->IsOpen());  
00078 
00079   char *run=strtok(runL.Data()," "); // init 'strtok'
00080 
00081   float sumIn=0;
00082   int nRun=0, nLum=0; // counts runs & scaler data for them
00083   do { // loop over runs
00084      printf("\n\n process run %02d '%s' \n",nRun,run);
00085     updatePattern(run);
00086 
00087     //.................. access LCP histos
00088     TString lcpHist=wrkDir+run+".hist.root";
00089     TFile * lcp=new TFile(lcpHist);
00090      assert(lcp->IsOpen());
00091     if(!lcp->IsOpen()) {
00092       printf("#fail %s does not open, skip\n",lcpHist.Data());
00093       continue;
00094     }
00095     nRun++;
00096 
00097 
00098     // ................  access scale data (if exist)
00099     TH1F *hSc=(TH1F*)scal->Get(run);
00100     if(hSc) { 
00101       printf(" use absolute normalization\n");
00102       sumScal(hSc); 
00103       calcNorm(run);
00104       nLum++;
00105     }  else { // do not scal by lumin
00106       printf(" NO Scaler data, --> no renormalization\n");
00107       for(int i=0;i<mxPol;i++) {Lum[i]=Lum2[i]=1.;}
00108     }
00109     
00110     //lcp->ls();
00111 
00112     TList *top = lcp->GetListOfKeys();
00113     TIter iter(top);
00114     TObject*ob=0;
00115     // ........... loop over histo with LCP
00116     while( ob=iter())     {
00117       char *name=ob->GetName();
00118       if(strstr(name,"PhiBx")==0) continue;
00119       // if(strstr(name,"PhiBxPt1")==0 && strstr(name,"PhiBxAll")==0) continue;
00120 
00121       TH2F *hIn=(TH2F*)lcp->Get(name);
00122       printf("  %10s integral=%.1f\n",name, hIn->Integral());
00123       if(strstr(name,"PhiBxAll"))sumIn+=hIn->Integral();
00124       
00125       if(nRun==1) createHisto1D(hIn); 
00126       TH1F * h[mxPol];
00127       TString cut=name+5;
00128       fetch1D(cut,h);
00129       accumulateLcp(hIn,h); 
00130 
00131     } // end of loop over histos in one run
00132     lcp->Close();
00133     delete lcp;
00134     
00135     
00136  } while(run=strtok(0," "));  // advance by one nam
00137 
00138   printf("acumulation done nRun=%d nLum=%d\n", nRun, nLum);
00139   assert(nLum==0 || nRun==nLum);
00140 
00141   if(nLum==0) autoLum("All");
00142   printf("left2=%s=\n",cutL.Data());
00143   printf("#fill %s , inp= %d ,  ",fill.Data(),sumIn);
00144   float sumAll=summaryLcp();
00145   printf(" AllEff= %.2f\n",sumAll/sumIn );
00146   outH->Write();
00147 
00148   //printLcp("PhiBxAll");
00149 
00150 
00151   //  printLcp("PhiBxPtM");
00152 
00153   // outH->ls();
00154   (outH->Get("PhiB+Y-All"))->Draw();
00155   TH1F*h1=(TH1F*) outH->Get("PhiB+Y+All");
00156   h1->Draw("same"); h1->SetLineColor(kRed);
00157   
00158 
00159   // hsc->Draw();
00160 
00161 }
00162 
00163 
00164 
00165 //-------------------------------------
00166 //-------------------------------------
00167 int bx2pol(int bx1) {
00168   if(bXselection[bx1-1]!='.') return -3;
00169   int pol= polPattBY[bx1-1];
00170   assert(pol<mxPol);
00171   return pol;
00172 }
00173 
00174 //-------------------------------------
00175 //-------------------------------------
00176 void   createHisto1D(TH2F *hIn){
00177   TString name=hIn->GetName();
00178   TString cut=name.Data()+5;
00179   cutL+=" "+cut;
00180   printf("create outHist for %s\n",name.Data());
00181   TAxis *phiAx=hIn->GetYaxis();
00182 
00183   outH->cd();
00184   int k;
00185   for(k=0;k<mxPol;k++) {
00186     TString name1="Phi";
00187     name1+=cpolBY[k];
00188     name1+=cut;
00189     char tit2[100]={"aaa bbb"};
00190     sprintf(tit2,"LCP vs. Phi/rad, pol=%s, cut=%s",cpolBY[k],cut.Data());
00191     printf("k=%d %s '%s'\n",k,name1.Data(), tit2);
00192     TH1F *h=new TH1F(name1,tit2,phiAx->GetNbins(),phiAx->GetXmin(), phiAx->GetXmax());
00193     h->Sumw2(); // be carefull, not use Fill(x,w) ! 
00194   }
00195 }
00196 
00197 
00198 //-------------------------------------
00199 //-------------------------------------
00200 void   accumulateLcp(TH2F *hIn, TH1F**h){
00201   TString name=hIn->GetName();
00202   //  printf("accu histo %s\n",name.Data());
00203 
00204   // fetch proper 1D histo
00205   int k;
00206   for(k=0;k<mxPol;k++) {
00207     TString name1=hIn->GetName();
00208     name1.ReplaceAll("Bx",cpolBY[k]);    
00209     h[k]=(TH1F *)outH->Get(name1);
00210     assert(h[k]);
00211   }
00212 
00213   assert(hIn);
00214   TAxis *X=hIn->GetXaxis();  
00215   int nx=X->GetNbins();
00216   TAxis *Yax=hIn->GetYaxis();
00217   assert(Yax);
00218   int ny=Yax->GetNbins();
00219   int ix,iy;
00220   
00221   double nSel=0;
00222   for(iy=1; iy<=ny; iy++) {
00223     float phi=Yax->GetBinCenter(iy);
00224     //intf("iy=%d phi/deg=%f\n", iy,phi/3.1416*180);
00225     for(ix=1;ix<=nx;ix++) {
00226       int pol=bx2pol(ix);
00227       if(pol<0) continue;
00228       double val=hIn->GetBinContent(ix,iy);
00229       if (val==0) continue;
00230       nSel+=val;
00231       TH1F * h1=h[pol];
00232       h1->AddBinContent(iy,val/Lum[pol]);
00233       // do error propagation by hand
00234       double err=h1->GetBinError(iy);
00235       h1->SetBinError(iy,sqrt(err*err+val/Lum2[pol]));
00236       // printf("add %d %f %d\n",pol,val, h[pol]->Integral());
00237     }
00238   } 
00239 
00240   // cross checks
00241 
00242   return;
00243   printf("%s --> %.1f nSel=%.1f\n",hIn->GetName(), hIn->Integral(),nSel);
00244   for(k=0;k<mxPol;k++) {
00245     TH1F *h1=h[k];
00246     printf("%s --> %.0f\n",h1->GetName(), h1->Integral());
00247   }
00248  
00249 }
00250 
00251 
00252 //-------------------------------------
00253 //-------------------------------------
00254 void   printLcp(char *cut){
00255   printf("print histo %s, format: kPhi  polBY ++ +-  -+ --\n#tot: ",cut);
00256   
00257   TH1F * h[mxPol];
00258 
00259   // fetch proper 1D histo
00260   int k;
00261   double sum=0;
00262   for(k=0;k<mxPol;k++) {
00263     TString name1=cut;
00264     name1.ReplaceAll("Bx",cpolBY[k]);    
00265     h[k]=(TH1F *)outH->Get(name1);
00266     assert(h[k]);
00267     sum+=h[k]->Integral();
00268     printf(" %.1f ", h[k]->Integral());
00269   }
00270 
00271   printf("\n#rel: ",cut);
00272   for(k=0;k<mxPol;k++) {
00273     printf(" %.3f ", h[k]->Integral()*4/sum);
00274   }
00275   printf("\n");
00276   TAxis *X=h[0]->GetXaxis();  
00277   int nx=X->GetNbins();
00278   int ix;
00279   
00280   for(ix=1; ix<=nx; ix++) {
00281     printf("%2d ",ix);
00282     for(k=0;k<mxPol;k++) {
00283       float val=h[k]->GetBinContent(ix);
00284       float err=h[k]->GetBinError(ix);
00285       printf("%8.1f +/- %5.1f,  ",val,err);
00286     }
00287     printf("\n");
00288   } 
00289 }
00290 
00291 //-------------------------------------
00292 //-------------------------------------
00293 void   fetch1D(TString cut,TH1F **h) {// fetch proper 1D histo
00294   int k;
00295   for(k=0;k<mxPol;k++) {
00296     TString name1="Phi";
00297     name1+=cpolBY[k];
00298     name1+=cut;
00299     h[k]=(TH1F *)outH->Get(name1);
00300     assert(h[k]);
00301   }
00302 }
00303 
00304 //-------------------------------------
00305 //-------------------------------------
00306 void sumScal(TH1F *h) {
00307   memset(Lum,0,sizeof(Lum)); // clear normalization
00308   memset(Lum2,0,sizeof(Lum2)); 
00309 
00310   assert(h);
00311   TAxis *X=h->GetXaxis();  
00312   int nx=X->GetNbins();
00313   int ix;
00314   int nBx=0;
00315   for(ix=1;ix<=nx;ix++) {
00316     int pol=bx2pol(ix);
00317     if(pol<0) continue;
00318     assert(pol<mxPol);
00319     Lum[pol]+=h->GetBinContent(ix);
00320     nBx++;
00321     // printf("%3d %c %2d  %8.1f %d \n",ix,bXselection[ix-1],polPattBY[ix-1], h->GetBinContent(ix),nBx);
00322   }
00323   
00324   int i;
00325   printf("Lum: ");
00326   for(i=0;i<mxPol;i++)
00327     printf("%d ",Lum[i]);
00328  
00329   printf(" usedBx=%d\n",nBx);
00330 
00331 }
00332   
00333 //-------------------------------------
00334 //-------------------------------------
00335 void calcNorm(char *name){
00336   int i;
00337   float sum=0;
00338   printf("\nyield ");
00339   for(i=0;i<mxPol;i++){
00340     sum+=Lum[i];
00341     printf("%d ",Lum[i]);
00342   }
00343   
00344   int id=atoi(name+1);
00345   //float day=(id-3000000.)/1000.;
00346 
00347   printf("\n#norm %s %d  ",name , id);
00348   sum/=4.;
00349   // renormalize
00350   for(i=0;i<mxPol;i++){
00351     Lum[i]/=sum;
00352     Lum2[i]=Lum[i]*Lum[i];
00353     printf("%.4f ",Lum[i]);
00354   }
00355   printf("\n");
00356 }
00357 
00358 
00359 
00360 //-------------------------------------
00361 //-------------------------------------
00362 void autoLum( char *cut0) {
00363   printf("autoLuminosiy using cut=%s\n",cut0);
00364   
00365   TH1F * h[mxPol];
00366   fetch1D(cut0,h);
00367   int k;
00368   for(k=0;k<mxPol;k++) {
00369     Lum[k]=h[k]->Integral();
00370   }// end of loop over pol  
00371   
00372   calcNorm("AutoLum7777");
00373 
00374   // scale all histos for all cuts 
00375   TString cut1=" "+cutL;
00376   char *cut=strtok(cut1.Data()," "); // init 'strtok'
00377   do {
00378     TH1F * h[mxPol];
00379     fetch1D(cut,h);
00380     int k;
00381     for(k=0;k<mxPol;k++) {
00382       h[k]->Scale(1./Lum[k]);
00383     }// end of loop over pol  
00384   } while(cut=strtok(0," "));  // advance by one nam
00385 }  
00386 
00387 
00388 //-------------------------------------
00389 //-------------------------------------
00390 float summaryLcp(){
00391   TString cut1=cutL;
00392   //printf("summary for cutL=%s\n",cut1.Data());
00393   char *cut=strtok(cut1.Data()," "); // init 'strtok'
00394   int sumAll=-1;
00395   do {
00396     // printf("cut:  %s\n",cut);
00397     float sum[mxPol];
00398     float sum0=0;
00399     TH1F * h[mxPol];
00400     fetch1D(cut,h);
00401     int k;
00402     for(k=0;k<mxPol;k++) {
00403       // printf("  %s %.1f -->",name1.Data(),h->Integral());
00404       sum[k]=h[k]->Integral();
00405       sum0+=sum[k];
00406     }// end of loop over pol  
00407     
00408     printf("%s= %.1f , ",cut,sum0);
00409     if(strstr(cut,"All"))sumAll=sum0;
00410 
00411   } while(cut=strtok(0," "));  // advance by one nam
00412   
00413   return sumAll;
00414 }
00415 
00416 
00417 
00418 
00419 
00420 
00421 
00422