StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
genWalgoYields.C
1 // simulates population of spin dependent W yields, including background
2 // version 2.0
3 /* mapping of spin4-index to helicity at STAR */
4 enum { ka=10, /* STAR pol B+ Y + */
5  kb=9, /* STAR pol B+ Y - */
6  kc=6, /* STAR pol B- Y + */
7  kd=5, /* STAR pol B- Y - */
8  mxSS=4,
9  mxEta=8, /* # of used STAR eta bins */
10  mxQ=2 /*0=W+, 1=W- */
11 };
12 
13 int kA[mxSS]={ka,kb,kc,kd}; // mapping of of valid spin4 values to blue & yellow polarization
14 
15 TRandom3 *rnd=0;
16 
17 //=================================
18 void genWalgoYields(int nEve=6e2, double pol1=0.5,double pol2=0.99 , double lumiSpread=0.1){
19  gROOT->LoadMacro("setWALmodel.C");
20  TString outF="run12toySetB";
21 
22  rnd=new TRandom3();
23  //rnd->SetSeed(1234); // fix random sequence
24  rnd->SetSeed(0); // job-unique random sequence
25 
26  // create output histos
27  TH1F * hLumi=newSpin4Histo("AEta8spinY1","relative lums from toy MC");
28 
29  TH1F * hWyield[mxQ][mxEta];
30  for (int iq=0;iq<mxQ;iq++) {
31  TString Q="P"; if (iq) Q="N";
32  for(int k=0;k<mxEta;k++) {
33  int etaBin=k+1;
34  TString cEta=Form("Eta%d",k+1);
35  hWyield[iq][k]=newSpin4Histo("A"+cEta+"spinY2_"+Q,"charge="+Q+", thrown W+BG events from toy MC, starEtaBin="+cEta);
36  //printf("k=%d iq=%d\n",k,iq);
37  }
38  }
39 
40  //------- generate lumi model -------
41  double lumi4[mxSS], prob4[mxSS];
42  for(int is=0; is<mxSS; is++) {
43  float val=1.+rnd->Uniform(lumiSpread);
44  lumi4[is]=val;
45  hLumi->SetBinContent(kA[is]+1,val);
46  }
47  TString core="2012";
48  TString bgName="8.15.12/bkgdBetaModel_"+core+".hist.root";
49  TFile *fdBG=new TFile(bgName); assert(fdBG->IsOpen());
50  TFile *fdWAL=new TFile("WALModel.hist.root"); assert(fdWAL->IsOpen());
51 
52  printf(" generate W algo yields, nEve=%d year=%s pol1=%.2f, pol2=%.2f lumiSpread=%.2f\n",nEve,core.Data(), pol1, pol2,lumiSpread);
53 
54  //fdWAL->ls(); // fdBG->ls();
55 
56  for( int iq=0; iq<mxQ; iq++) {
57  TString Q="P"; if (iq) Q="N";
58  // if(iq==1) continue; // tmp work only with + charge
59 
60  // Ws
61  TH1F *hWAL=fdWAL->Get("modW_AL_"+Q); assert(hWAL);
62  TH1F *hWALL=fdWAL->Get("modW_ALL_"+Q); assert(hWALL);
63  // background
64  // TH1F *hAlpha=fdBG->Get("alpha"+Q+"_"+core); assert(hAlpha);
65  TH1F *hBeta=fdBG->Get("beta"+Q+"_"+core); assert(hBeta);
66 
67  // loopp over eta bins mapped to different eta bins
68  for( int starPhysEtaBin=1; starPhysEtaBin<=6; starPhysEtaBin++) {
69  int nThrow=nEve/6 *(1.-iq*0.3);
70  if (starPhysEtaBin>=5) nThrow/=5;
71  printf("\n\n........ processing Q=%s ..... starEtaBin=%d ....nThrow=%d\n",Q.Data(),starPhysEtaBin,nThrow);
72 
73  //map starEta on to polBeam eta
74  int polBeam1EtaBin=12+starPhysEtaBin;
75  int polBeam2EtaBin=17-starPhysEtaBin;
76  printf("starPhysEtaBin=%d --> polBeam1=%d polBeam2=%d\n",starPhysEtaBin,polBeam1EtaBin,polBeam2EtaBin);
77 
78  // compute probabilities ....
79  double WAL1=hWAL->GetBinContent(polBeam1EtaBin);
80  double WAL2=hWAL->GetBinContent(polBeam2EtaBin);
81  double WALL=hWALL->GetBinContent(polBeam1EtaBin); // does not matter which etaBin - it is even-function
82 
83  // double alpha1=hAlpha->GetBinContent(polBeam1EtaBin);
84  // double alpha2=hAlpha->GetBinContent(polBeam2EtaBin);
85  double beta=hBeta->GetBinContent(starPhysEtaBin);
86 
87  // tmp
88  double alpha1=0, alpha2=0;
89  // beta=starPhysEtaBin/10.;
90 
91  //------- generate pol x-section model ------
92  getProbOneEta( pol1, pol2,WAL1, WAL2, WALL, alpha1,alpha2, beta, lumi4, prob4);
93 
94  //---------- generate events ---------
95  for( int ieve=0; ieve<nThrow; ieve++) {
96  int spin4=throwSpin4( prob4 );
97  int kBin=starPhysEtaBin-1;
98  //printf("eve=%d got spin4=%d kBin=%d\n",ieve,spin4,kBin);
99  hWyield[iq][kBin]->Fill(spin4);
100 
101  }
102  }//..... end of loop over eta bins
103 
104  // compute yields for B-only and E-only bins, no overlap assumed
105  // Endcap : starPhysEtaBin7=5+6
106  hWyield[iq][7-1]->Add( hWyield[iq][6-1]);
107  hWyield[iq][7-1]->Add( hWyield[iq][5-1]);
108  for(int jj=0;jj<4;jj++)
109  hWyield[iq][8-1]->Add( hWyield[iq][jj]);
110 
111  }//....... end of loop over + - charge
112 
113  TFile *fdOut=new TFile(outF+".wana.hist.root","recreate"); assert(fdOut->IsOpen());
114  printf("saving histo w/ events to %s\n",fdOut->GetName());
115  // loopp over all eta bins
116  for(int k=0;k<mxEta;k++) {
117  TString cEta=Form("Eta%d",k+1);
118  fdOut->mkdir(cEta);
119  fdOut->cd(cEta);
120  for (int iq=0;iq<mxQ;iq++) hWyield[iq][k]->Write();
121  //fdOut->ls();
122  fdOut->cd();
123  }
124 
125  fdOut->cd("Eta8");
126  hLumi->Write();
127 
128  //fdOut->ls();
129  fdOut->Close();
130 
131 }
132 
133 
134 //=====================
135 int throwSpin4(double *prob4) { // assigne event to one of 4 spin states
136  double x=rnd->Uniform();
137  for(int is=0; is<mxSS; is++)
138  if(x<prob4[is]) break;
139  assert(is<mxSS); // probability should not exeed 2.
140  int spin4=kA[is];
141  return spin4;
142 }
143 
144 //=====================
145 void getProbOneEta(double P1, double P2, double WAL1, double WAL2, double WALL,
146  double alpha1, double alpha2, double beta, double *lumi,
147  double *prob ) { // output
148  printf("getProbOneEta INPUT: pol1=%.2f, pol2=%.2f, W: AL1=%g AL2=%g ALL=%g \n",P1,P2,WAL1,WAL2,WALL);
149  printf("input lumis: "); for(int is=0;is<mxSS; is++) printf("%8.3f ", lumi[is]);printf("\n");
150  if(beta<0.1) {printf("changed beta WARN\n"); beta=0.9;}
151  printf("input BG: beta=%.3f alpha1=%.3f alpha2=%.3f\n",beta,alpha1, alpha2);
152 
153  double ppW , pnW, npW, nnW ; // spin dependent x-sections, relative
154 
155  ppW=lumi[0]*(1. +(WAL1*beta +alpha1)*P1 +( WAL2*beta +alpha2)*P2 +WALL*beta*P1*P2 );
156  pnW=lumi[1]*(1. +(WAL1*beta +alpha1)*P1 -( WAL2*beta +alpha2)*P2 -WALL*beta*P1*P2 );
157  npW=lumi[2]*(1. -(WAL1*beta +alpha1)*P1 +( WAL2*beta +alpha2)*P2 -WALL*beta*P1*P2 );
158  nnW=lumi[3]*(1. -(WAL1*beta +alpha1)*P1 -( WAL2*beta +alpha2)*P2 +WALL*beta*P1*P2 );
159 
160  // get cumulative probabilities after folding in lumi asymmetries
161  prob[0]= 0. +ppW;
162  prob[1]=prob[0]+pnW;
163  prob[2]=prob[1]+npW;
164  prob[3]=prob[2]+nnW;
165 
166  for(int is=0;is<mxSS; is++) prob[is]/=prob[3]; // renormalize probabilities
167 
168  printf("computed probs:%8.3f ", prob[0]); for(int is=1;is<mxSS; is++) printf("%8.3f ", prob[is]-prob[is-1]); printf("\n");
169 
170  // Returns relative probabilities for 4 spin states for 'signal' yields
171  // (including pol & unpol backgrounds)
172 }
Generated by   doxygen 1.8.5