StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ParityDevCorrFctn.cxx
1 #include "StHbtMaker/CorrFctn/ParityDevCorrFctn.h"
2 #include "StHbtMaker/Infrastructure/StHbtAnalysis.h"
3 #include "StHbtMaker/Cut/mikesEventCut.h"
4 
5 #include <cstdio>
6 
7 #ifdef __ROOT__
8 ClassImp(ParityDevCorrFctn)
9 #endif
10 
11 //____________________________
12 ParityDevCorrFctn::ParityDevCorrFctn(){
13  //mTagWriter = StHbtTagWriter::Instance(); // get the singleton
14 
15  // histograms
16  mKMiSamex = new StHbt1DHisto("SameKMinusX","Parity Same KMinus x",100,-0.4,0.4);
17  mKMiSamey = new StHbt1DHisto("SameKMinusY","Parity Same KMinus y",100,-0.4,0.4);
18  mKMiSamez = new StHbt1DHisto("SameKMinusZ","Parity Same KMinus z",100,-0.4,0.4);
19  mJcSamex = new StHbt1DHisto("SameJcX","Parity Same Jc X",100,-0.2,0.2);
20  mJcSamey = new StHbt1DHisto("SameJcY","Parity Same Jc Y",100,-0.2,0.2);
21  mJcSamez = new StHbt1DHisto("SameJcZ","Parity Same Jc Z",100,-0.2,0.2);
22  mJcKtSame = new StHbt1DHisto("SameJcKt","Parity Same JcKt",100,-0.005,0.005);
23  mKtwistKtSame = new StHbt1DHisto("SameKtwistKt","Parity Same KtwistKt",100,-0.005,0.005);
24 
25  mKMiMixedx = new StHbt1DHisto("MixedKMinusX","Parity Mixed KMinus x",100,-0.4,0.4);
26  mKMiMixedy = new StHbt1DHisto("MixedKMinusY","Parity Mixed KMinus y",100,-0.4,0.4);
27  mKMiMixedz = new StHbt1DHisto("MixedKMinusZ","Parity Mixed KMinus z",100,-0.4,0.4);
28  mJcMixedx = new StHbt1DHisto("MixedJcX","Parity Mixed Jc X",100,-0.2,0.2);
29  mJcMixedy = new StHbt1DHisto("MixedJcY","Parity Mixed Jc Y",100,-0.2,0.2);
30  mJcMixedz = new StHbt1DHisto("MixedJcZ","Parity Mixed Jc Z",100,-0.2,0.2);
31  mJcKtMixed = new StHbt1DHisto("MixedJcKt","Parity Mixed JcKt",100,-0.005,0.005);
32  mKtwistKtMixed = new StHbt1DHisto("MixedKtwistKt","Parity Mixed KtwistKt",100,-0.005,0.005);
33 
34  mJcKtBinomial = new StHbt1DHisto("JcKtBin","Number JcKt positive",100,-0.10,0.10);
35  mNumPairsBin = new StHbt1DHisto("Numpairs","Number of Pairs",2000,000.,2000.);
36 
37 
38  // end histograms
39 
40  // this next bit is unfortunately needed so that we can have many histos of same "title"
41  // it is neccessary if we typedef StHbt1DHisto to TH1d (which we do)
42  mKMiSamex->SetDirectory(0);
43  mKMiSamey->SetDirectory(0);
44  mKMiSamez->SetDirectory(0);
45  mJcSamex->SetDirectory(0);
46  mJcSamey->SetDirectory(0);
47  mJcSamez->SetDirectory(0);
48  mJcKtSame->SetDirectory(0);
49  mKtwistKtSame->SetDirectory(0);
50  mKMiMixedx->SetDirectory(0);
51  mKMiMixedy->SetDirectory(0);
52  mKMiMixedz->SetDirectory(0);
53  mJcMixedx->SetDirectory(0);
54  mJcMixedy->SetDirectory(0);
55  mJcMixedz->SetDirectory(0);
56  mJcKtMixed->SetDirectory(0);
57  mKtwistKtMixed->SetDirectory(0);
58  mJcKtBinomial->SetDirectory(0);
59  mNumPairsBin->SetDirectory(0);
60 
61  mKMiMixedx->Sumw2();
62  mKMiMixedy->Sumw2();
63  mKMiMixedz->Sumw2();
64  mJcMixedx->Sumw2();
65  mJcMixedy->Sumw2();
66  mJcMixedz->Sumw2();
67  mJcKtMixed->Sumw2();
68  mKtwistKtMixed->Sumw2();
69  mKMiSamex->Sumw2();
70  mKMiSamey->Sumw2();
71  mKMiSamez->Sumw2();
72  mJcSamex->Sumw2();
73  mJcSamey->Sumw2();
74  mJcSamez->Sumw2();
75  mJcKtSame->Sumw2();
76  mKtwistKtSame->Sumw2();
77  mJcKtBinomial->Sumw2();
78  mNumPairsBin->Sumw2();
79 
80 }
81 
82 //____________________________
83 ParityDevCorrFctn::~ParityDevCorrFctn(){
84  // histograms
85  delete mKMiSamex;
86  delete mKMiSamey;
87  delete mKMiSamez;
88  delete mJcSamex ;
89  delete mJcSamey ;
90  delete mJcSamez ;
91  delete mJcKtSame ;
92  delete mKtwistKtSame ;
93 
94  delete mKMiMixedx;
95  delete mKMiMixedy;
96  delete mKMiMixedz;
97  delete mJcMixedx ;
98  delete mJcMixedy ;
99  delete mJcMixedz ;
100  delete mJcKtMixed ;
101  delete mKtwistKtMixed ;
102  // end histograms;
103 }
104 //_________________________
105 void ParityDevCorrFctn::Finish(){
106  cout << " alive in finish " << endl;
107 
108 }
109 //____________________________
110 StHbtString ParityDevCorrFctn::Report(){
111  string stemp = "Parity Correlation Function Report:\n";
112  StHbtString returnThis = stemp;
113  return returnThis;
114 }
115 //____________________________
116 inline void ParityDevCorrFctn::AddRealPair(const StHbtPair* pair){
117  cout << "WARNING ParityDevCorrFctn::AddRealPair shouldn't be called" << endl;
118 }
119 //_________________________
120 
121 //_________________________
122 inline void ParityDevCorrFctn::AddMixedPair(const StHbtPair* pair){
123  cout << "WARNING ParityDevCorrFctn::AddMixedPair shouldn't be called" << endl;
124 }
125 
126 //_________________________
127 inline void ParityDevCorrFctn::ParityCompute(ParityBuff *Plus, ParityBuff *Minus, int mxd){
128 
129  StHbtLorentzVector PlusTrack;
130  StHbtLorentzVector MinusTrack;
131  StHbtThreeVector KMinusPos;
132  StHbtThreeVector KMinusNeg;
133  StHbtThreeVector KMinus;
134  StHbtThreeVector Kt;
135 
136  double JcKt =0;
137  double PairJcKt = 0;
138  double JcKtBinom = 0;
139  double KtwistKt = 0;
140 
141  int plusSize = Plus->size();
142  int minusSize = Minus->size();
143 
144  if (mxd == SAME) {
145  cout << "******** we got to the ParityCompute for SAME event" << endl;
146  cout << " size of same event plus is " << Plus->size() << " ";
147  cout << " size of same event minus is " << Minus->size() << endl;
148  }
149  if (mxd == MIXED) {
150  cout << "******** we got to the ParityCompute for MIXED event" << endl;
151  cout << " size of mixed event plus is " << Plus->size() << " ";
152  cout << " size of mixed event minus is " << Minus->size() << endl;
153  }
154  // ****** Compute K- **********
155 
156  {for (int jjj = 0; jjj <plusSize; jjj++){
157  KMinusPos += (*Plus)[jjj].vect().unit();
158  }}
159  {for (int jjj = 0; jjj < minusSize; jjj++){
160  KMinusNeg += (*Minus)[jjj].vect().unit();
161  }}
162  KMinus = (KMinusPos/plusSize) - (KMinusNeg/minusSize);
163  Kt = KMinus;
164  Kt.setZ(0.);
165  cout << "Kt is " << Kt.x() << "," << Kt.y() << endl ;
166 
167  // now calculate Ktwist*Kt
168 
169  double plusSum = 0.;
170  {for (int jjj = 0; jjj < plusSize; jjj++){
171  StHbtThreeVector TempV = (*Plus)[jjj].vect();
172  double step1 = (TempV.y() * Kt.x()-TempV.x() * Kt.y() ); // this is y-component rotated to Kt coordinates
173  double step2 = step1*TempV.z(); // multiply by z-component
174  double step3 = step2/(TempV.mag2() * Kt.mag()); // normalize
175  plusSum += step3;
176  }}
177  double minusSum = 0.;
178  {for (int jjj = 0; jjj < minusSize; jjj++){
179  StHbtThreeVector TempV = (*Minus)[jjj].vect();
180  double step1 = (TempV.y() * Kt.x()-TempV.x() * Kt.y() ); // this is y-component rotated to Kt coordinates
181  double step2 = step1*TempV.z(); // multiply by z-component
182  double step3 = step2/(TempV.mag2() * Kt.mag()); // normalize
183  minusSum += step3;
184  }}
185 
186  double Ktwist = (plusSum/ plusSize) - (minusSum/ minusSize);
187  KtwistKt = Ktwist*Kt.mag();
188 
189  cout << "KtwistKt is " << KtwistKt << endl;
190  // end Ktwist
191 
192  // now JcKt -here we use unique pairs since it's approx. as sensitive as using all pairs (?)
193 
194  int smallSize = plusSize;
195  int numPairs = 0;
196  int numJcKtPlus = 0;
197  if (minusSize < plusSize) smallSize = minusSize;
198  StHbtThreeVector Jc;
199  {for (int jjj = 0; jjj < smallSize; jjj++){
200  PlusTrack = (*Plus)[jjj];
201  MinusTrack = (*Minus)[jjj];
202  StHbtThreeVector vOne = PlusTrack.vect();
203  StHbtThreeVector vTwo = MinusTrack.vect();
204  StHbtThreeVector PairJc = vOne.unit().cross(vTwo.unit());
205  if ( (vOne.z()*vTwo.z()) < 0.){
206  PairJc.setX(-PairJc.x());
207  }
208  if ( (vOne.x()*vTwo.x()) < 0.){
209  PairJc.setY(-PairJc.y());
210  }
211  if ( (vOne.y()*vTwo.y()) < 0.){
212  PairJc.setZ(-PairJc.z());
213  }
214  Jc += PairJc;
215  //this is for binomial only
216  numPairs++;
217  PairJcKt = PairJc.dot(Kt);
218  if (PairJcKt > 0){
219  numJcKtPlus++;
220  }
221  // end binomial
222  }}
223  Jc = Jc/(smallSize); // remember to change normalization if we go back to all pairs
224  JcKt = Jc.dot(Kt);
225  cout << "JcKt is " << JcKt << endl;
226  JcKtBinom = ( double(numJcKtPlus)/double(numPairs) ) -1./2. ;
227  cout << "JcKt bin = " << numJcKtPlus << "/" <<numPairs << "-1/2 ="<< JcKtBinom << endl;
228  // end JcKt
229 
230  if (mxd == SAME) {
231  mKMiSamex->Fill(KMinus.x());
232  mKMiSamey->Fill(KMinus.y());
233  mKMiSamez->Fill(KMinus.z());
234  mJcSamex->Fill(Jc.x());
235  mJcSamey->Fill(Jc.y());
236  mJcSamez->Fill(Jc.z());
237  mJcKtSame->Fill(JcKt);
238  mKtwistKtSame->Fill(KtwistKt);
239 
240  mJcKtBinomial->Fill(JcKtBinom); // binomial
241  mNumPairsBin->Fill(numPairs); // binomial
242  }
243 
244  if (mxd == MIXED) {
245  mKMiMixedx->Fill(KMinus.x());
246  mKMiMixedy->Fill(KMinus.y());
247  mKMiMixedz->Fill(KMinus.z());
248  mJcMixedx->Fill(Jc.x());
249  mJcMixedy->Fill(Jc.y());
250  mJcMixedz->Fill(Jc.z());
251  mJcKtMixed->Fill(JcKt);
252  mKtwistKtMixed->Fill(KtwistKt);
253  }
254 
255 
256 }
257 //_________________________
Generated by   doxygen 1.8.5