StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStruct2ptCorrelations.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStruct2ptCorrelations.cxx,v 1.31 2013/02/08 19:32:43 prindle Exp $
4  *
5  * Author: Jeff Porter adaptation of Aya's 2pt-analysis
6  *
7  **********************************************************************
8  *
9  * Description: Analysis code for 2pt-analysis.
10  * The analysis runs as follows:
11  * 1D and 2D arrays (in yt,eta,phi) are setup
12  * and filled for each of the 8 pair types:
13  * Sibling (++, +-, -+, --)
14  * Mixed (++, +-, -+, --)
15  * Particle id is done via dEdx and introduced into the analysis via cut-bins.
16  * Order particles so pi always come before K before p and plus comes before
17  * minus. 2D histograms will not be guaranteed to be symmetric.
18  * The 2D versions are additionally divided into z-vertex (via the StEStructBuffer)
19  * After arrays are filled (looped over all events/job), Histograms are
20  * created, filled, and written out to the data file for further
21  * processing.
22  *
23  * Note that if we find charge symmetry we can form LS, US, CD and CI
24  * combinations using these histograms.
25  *
26  *
27  ***********************************************************************/
28 #include "StEStruct2ptCorrelations.h"
29 
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH1D.h"
33 #include "TH2D.h"
34 #include "TFile.h"
35 
36 #include "StEStructPool/EventMaker/StEStructEvent.h"
37 #include "StEStructPool/EventMaker/StEStructTrack.h"
38 #include "StEStructPool/Correlations/StEStructMaxB.h"
39 #include "StEStructPool/Correlations/StEStructBuffer.h" // to read some constants
40 
41 #include "StTimer.hh"
42 #include "StEStructCutBin.h"
43 
44 #include "Stiostream.h"
45 
46 #include "StMessMgr.h"
47 
48 #include <iostream>
49 using namespace std;
50 
51 ClassImp(StEStruct2ptCorrelations);
52 
53 //--------------------------------------------------------------------------
54 StEStruct2ptCorrelations::StEStruct2ptCorrelations(int mode) {
55  initInternalData();
56  manalysisMode=mode;
57  mPairCuts = new StEStructPairCuts;
58  mHcb = NULL;
59 }
60 
61 
62 //--------------------------------------------------------------------------
63 StEStruct2ptCorrelations::StEStruct2ptCorrelations(StEStructPairCuts* pcuts, int mode) {
64  initInternalData();
65  manalysisMode=mode;
66  mPairCuts=pcuts;
67  mHcb = NULL;
68 }
69 
70 //----------------------------------------------------------
71 void StEStruct2ptCorrelations::initInternalData(){
72 
73  mPairCuts = NULL;
74  mQAHists = NULL;
75  mCurrentEvent = NULL;
76  mMixingEvent = NULL;
77  moutFileName = NULL;
78  mqaoutFileName = NULL;
79  mOneZBuffer = NULL;
80 
81  mskipPairCuts = false;
82  mdoPairCutHistograms = false;
83  mdoPairDensityHistograms = false;
84  mskipEtaDeltaWeight = false;
85  mdoInvariantMassHistograms = false;
86  mdoFillEtaEta = false;
87  mdoFillSumHistograms = false;
88  mdontFillMeanPt = false;
89  mdontFillYtYt = false;
90  mFillQInv = false;
91  mFillASSS = false;
92 
93  mInit = false;
94  mDeleted = false;
95  mHistosWritten = false;
96  mlocalQAHists = false;
97  manalysisIndex = -1;
98 
99  kZBuffMin = -75.0;
100  kZBuffMax = +75.0;
101  kBuffWidth = 5.0;
102  kNumBuffers = int( (kZBuffMax-kZBuffMin)/kBuffWidth ); // Better be 30.
103  // for(int i=0;i<kNumBuffers;i++)mbuffCounter[i]=0;
104  mZBufferCutBinning = 0;
105 }
106 
107 
108 
109 //--------------------------------------------------------------------------
110 StEStruct2ptCorrelations::~StEStruct2ptCorrelations(){ cleanUp(); };
111 
112 
113 void StEStruct2ptCorrelations::init() {
114 
115  cout << "Initializing with analysis mode " << manalysisMode << endl;
116  cout << "Use Z Buffer cut binning = " << mZBufferCutBinning << endl;
117 
118  mCurrentEvent=NULL;
119  mtimer=NULL;
120 
121  //--> code to simplify hist-creation via class held name defs
122  const char* _tmpName[]={"Sibpp","Sibpm","Sibmp","Sibmm","Mixpp","Mixpm","Mixmp","Mixmm"};
123  const char* _tmpTitle[]={"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-",
124  "Mixed : +.+","Mixed : +.-","Mixed : -.+","Mixed : -.-"};
125 
126  for(int i=0;i<8;i++){
127  bName[i]=new char[6];
128  strcpy(bName[i],_tmpName[i]);
129  bTitle[i]=new char[20];
130  strcpy(bTitle[i],_tmpTitle[i]);
131  }
132 
133  if (manalysisMode & 0x1) {
134  mskipPairCuts=true;
135  }
136  if (manalysisMode & 0x2) {
137  mdoPairCutHistograms=true;
138  }
139  if (manalysisMode & 0x4) {
140  mdoPairDensityHistograms=true;
141  }
142  if (manalysisMode & 0x8) {
143  mskipEtaDeltaWeight = true;
144  }
145  if (manalysisMode & 0x10) {
146  mdoInvariantMassHistograms = true;
147  }
148  if (manalysisMode & 0x20) {
149  mdoFillEtaEta = true;
150  }
151  if (manalysisMode & 0x40) {
152  mdoFillSumHistograms = true;
153  }
154  if (manalysisMode & 0x80) {
155  mdontFillMeanPt = true;
156  }
157  if (manalysisMode & 0x100) {
158  mdontFillYtYt = true;
159  }
160  if (manalysisMode & 0x200) {
161  mFillQInv = true;
162  }
163  if (manalysisMode & 0x400) {
164  mFillASSS = true;
165  }
166 
167  cout << " Skip Pair Cuts = " << mskipPairCuts << endl;
168  cout << " Do Pair Cut Hists = " << mdoPairCutHistograms << endl;
169  cout << " Do Pair Density Hists = " << mdoPairDensityHistograms << endl;
170  cout << " Skip EtaDelta weights = " << mskipEtaDeltaWeight << endl;
171  cout << " Do Invariant mass histograms = " << mdoInvariantMassHistograms << endl;
172  cout << " Fill EtaEta (and PhiPhi) = " << mdoFillEtaEta << endl;
173  cout << " Fill Sum Histograms (SYt, SEta) = " << mdoFillSumHistograms << endl;
174  cout << " Don't fill mean pt histograms = " << mdontFillMeanPt << endl;
175  cout << " Don't fill YtYt histograms = " << mdontFillYtYt << endl;
176  cout << " Fill QInv histograms = " << mFillQInv << endl;
177  cout << " Fill AS,SS (eta,eta) histograms = " << mFillASSS << endl;
178  if (mdoPairCutHistograms) {
179  cout << " >>>>> You have tried to turn on Pair Cut Histograms. Those have not been re-implemented yet. " << endl;
180  cout << " >>>>> If you really want them we should do something about it. " << endl;
181  mdoPairCutHistograms = false;
182  }
183 
184  for(int i=0;i<8;i++)numPairs[i]=numPairsProcessed[i]=mpossiblePairs[i]=0;
185 
186  initArrays();
187 // Try allocating histograms at start of job.
188 // If we don't add histograms to directory we don't get the obnoxious
189 // Potential memory leak error.
190  TH1::AddDirectory(kFALSE);
191  //initHistograms();
192 
193  /* Event count via Nch distribution */
194  char name[2048];
195  int nZBins = 1;
196  if (mZBufferCutBinning) {
197  nZBins = kNumBuffers;
198  }
199  mHNEventsSib = new TH1D*[nZBins];
200  mHNEventsMix = new TH1D*[nZBins];
201  mHNEventsPosSib = new TH1D*[nZBins];
202  mHNEventsPosMix = new TH1D*[nZBins];
203  mHNEventsNegSib = new TH1D*[nZBins];
204  mHNEventsNegMix = new TH1D*[nZBins];
205  for (int iz=0;iz<nZBins;iz++) {
206  sprintf(name,"NEventsSib_zBuf_%i",iz);
207  mHNEventsSib[iz]=new TH1D(name,name,1000,0.,2000.);
208  sprintf(name,"NEventsMix_zBuf_%i",iz);
209  mHNEventsMix[iz]=new TH1D(name,name,1000,0.,2000.);
210  sprintf(name,"NEventsPosSib_zBuf_%i",iz);
211  mHNEventsPosSib[iz]=new TH1D(name,name,1000,0.,2000.);
212  sprintf(name,"NEventsPosMix_zBuf_%i",iz);
213  mHNEventsPosMix[iz]=new TH1D(name,name,1000,0.,2000.);
214  sprintf(name,"NEventsNegSib_zBuf_%i",iz);
215  mHNEventsNegSib[iz]=new TH1D(name,name,1000,0.,2000.);
216  sprintf(name,"NEventsNegMix_zBuf_%i",iz);
217  mHNEventsNegMix[iz]=new TH1D(name,name,1000,0.,2000.);
218  }
219 
220  StEStructCutBin* cb = StEStructCutBin::Instance();
221  if (mReader->mTCuts) {
222  cb->setMaxDEta(mReader->mTCuts->maxVal("Eta")-mReader->mTCuts->minVal("Eta"));
223  }
224  int ncutbins = cb->getNumBins();
225  int nQAbins = cb->getNumQABins();
226 
227  // Do we want to make invariant mass histograms?
228  // (Only works in cutbin mode 5, which uses PID assignments.)
229  if (mdoInvariantMassHistograms) {
230  cb->setCutBinHistMode(1);
231  cb->initCutBinHists();
232  } else {
233  cb->setCutBinHistMode(0);
234  }
235 
236  /* QA histograms */
237  if(!mQAHists) {
238  mlocalQAHists = true;
239  cout<<"creating QA hists"<<endl;
240  mQAHists = new StEStructQAHists(); // if not set .. assume data
241  mQAHists->initTrackHistograms(nQAbins,analysisIndex());
242  } else {
243  cout<<"init QA Hists"<<endl;
244  mQAHists->initTrackHistograms(nQAbins); // will do it for only 1 analysis
245  }
246 
247  /* Inclusive pt distribution */
248  mHptAll = new TH1D("ptAll","ptAll",500,0,10.);
249 
250  /*
251  * pt distributions for parent distributions.
252  * Use to get estimate of mean pt of parent distributionb but also
253  * mean number of tracks in the sample.
254  */
255  int numParentBins=cb->getNumParentBins();
256  int nzb = 1;
257  if (mZBufferCutBinning) {
258  nzb = kNumBuffers;
259  }
260 
261  mHMeanPtTot = new TH1D*[numParentBins*nzb];
262  mHMeanPtP = new TH1D*[numParentBins*nzb];
263  mHMeanPtM = new TH1D*[numParentBins*nzb];
264  mHMeanYtTot = new TH1D*[numParentBins*nzb];
265  mHMeanYtP = new TH1D*[numParentBins*nzb];
266  mHMeanYtM = new TH1D*[numParentBins*nzb];
267  mHEtaTot = new TH1D*[numParentBins*nzb];
268  mHEtaP = new TH1D*[numParentBins*nzb];
269  mHEtaM = new TH1D*[numParentBins*nzb];
270  TString hname;
271  StEStructBinning* b = StEStructBinning::Instance();
272  if (mReader->mTCuts) {
273  b->setEtaRange(mReader->mTCuts->minVal("Eta"),mReader->mTCuts->maxVal("Eta"));
274  }
275  for(int p=0;p<numParentBins;p++){
276  for (int z=0;z<nzb;z++) {
277  int pz = p*nzb + z;
278  hname = "meanPtTot_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
279  mHMeanPtTot[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ptMin(),b->ptMax());
280  hname = "meanPtP_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
281  mHMeanPtP[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ptMin(),b->ptMax());
282  hname = "meanPtM_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
283  mHMeanPtM[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ptMin(),b->ptMax());
284  hname = "meanYtTot_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
285  mHMeanYtTot[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ytMin(),b->ytMax());
286  hname = "meanYtP_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
287  mHMeanYtP[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ytMin(),b->ytMax());
288  hname = "meanYtM_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
289  mHMeanYtM[pz] = new TH1D(hname.Data(),hname.Data(),200,b->ytMin(),b->ytMax());
290  hname = "etaTot_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
291  mHEtaTot[pz] = new TH1D(hname.Data(),hname.Data(),100,b->etaMin(),b->etaMax());
292  hname = "etaP_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
293  mHEtaP[pz] = new TH1D(hname.Data(),hname.Data(),100,b->etaMin(),b->etaMax());
294  hname = "etaM_parentBin"; hname += p; hname += "_zBuf_"; hname += z;
295  mHEtaM[pz] = new TH1D(hname.Data(),hname.Data(),100,b->etaMin(),b->etaMax());
296  }
297  }
298  hname = "ptTot_toward";
299  mHPtTot[0] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
300  hname = "ptP_toward";
301  mHPtP[0] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
302  hname = "ptM_toward";
303  mHPtM[0] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
304  hname = "ptTot_transverse";
305  mHPtTot[1] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
306  hname = "ptP_transverse";
307  mHPtP[1] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
308  hname = "ptM_transverse";
309  mHPtM[1] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
310  hname = "ptTot_away";
311  mHPtTot[2] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
312  hname = "ptP_away";
313  mHPtP[2] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
314  hname = "ptM_awayAway";
315  mHPtM[2] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
316  hname = "ptTot_awayAway";
317  mHPtTot[3] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
318  hname = "ptP_awayAway";
319  mHPtP[3] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
320  hname = "ptM_away";
321  mHPtM[3] = new TH2D(hname.Data(),hname.Data(),200,0,10,20,0,5);
322  hname = "ytTot_toward";
323  mHYtTot[0] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
324  hname = "ytP_toward";
325  mHYtP[0] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
326  hname = "ytM_toward";
327  mHYtM[0] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
328  hname = "ytTot_transverse";
329  mHYtTot[1] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
330  hname = "ytP_transverse";
331  mHYtP[1] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
332  hname = "ytM_transverse";
333  mHYtM[1] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
334  hname = "ytTot_away";
335  mHYtTot[2] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
336  hname = "ytP_away";
337  mHYtP[2] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
338  hname = "ytM_away";
339  mHYtM[2] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
340  hname = "ytTot_awayAway";
341  mHYtTot[3] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
342  hname = "ytP_awayAway";
343  mHYtP[3] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
344  hname = "ytM_awayAway";
345  mHYtM[3] = new TH2D(hname.Data(),hname.Data(),25,1,4.5,25,1,4.5);
346  mHPhiAssocTot = new TH2D("phiAssocTot","phiAssocTot",90,0,3.1415926,16,1,4.5);
347  mHPhiAssocP = new TH2D("phiAssocP","phiAssocP",90,0,3.1415926,16,1,4.5);
348  mHPhiAssocM = new TH2D("phiAssocM","phiAssocM",90,0,3.1415926,16,1,4.5);
349  mHPhiAssocPtTot = new TH2D("phiAssocPtTot","phiAssocPtTot",90,0,3.1415926,40,0,5);
350  mHPhiAssocPtP = new TH2D("phiAssocPtP","phiAssocPtP",90,0,3.1415926,40,0,5);
351  mHPhiAssocPtM = new TH2D("phiAssocPtM","phiAssocPtM",90,0,3.1415926,40,0,5);
352  mHPtTrigTot = new TH1D("PtTrigTot","PtTrigTot",40,0,5);
353  mHPtTrigP = new TH1D("PtTrigP","PtTrigP",40,0,5);
354  mHPtTrigM = new TH1D("PtTrigM","PtTrigM",40,0,5);
355  mHYtTrigTot = new TH1D("YtTrigTot","YtTrigTot",25,1,4.5);
356  mHYtTrigP = new TH1D("YtTrigP","YtTrigP",25,1,4.5);
357  mHYtTrigM = new TH1D("YtTrigM","YtTrigM",25,1,4.5);
358 
359 
360  /* Event Mixing Parameters */
361  mHMixZdN = new TH2D("Mixed_Z_dN","Event Mixing: average-Z vs delta-N",50,-25,25, 50,-25,25);
362  mHMixZN = new TH2D("Mixed_Z_N","Event Mixing: average-Z vs average-N",50,-25,25, 50,0,1500);
363  mHMixZdC = new TH2D("Mixed_Z_dC","Event Mixing: average-Z vs delta-C",50,-25,25, 50,-100,100);
364  mHMixZC = new TH2D("Mixed_Z_C","Event Mixing: average-Z vs average-C",50,-25,25,50,0,10000);
365  mHMixZdZ = new TH2D("Mixed_Z_dZ","Event Mixing: average-Z vs delta-Z",50,-25,25,50,-5,5);
366  mHMixdZdN = new TH2D("Mixed_dZ_dN","Event Mixing: delta-Z vs delta-N",50,-5,5, 50,-25,25);
367  mHMixdZN = new TH2D("Mixed_dZ_N","Event Mixing: delta-Z vs average-N",50,-5,5, 50,0,1500);
368  mHMixdZdC = new TH2D("Mixed_dZ_dC","Event Mixing: delta-Z vs delta-C",50,-5,5, 50,-100,100);
369  mHMixdZC = new TH2D("Mixed_dZ_C","Event Mixing: delta-Z vs average-C",50,-5,5,50,0,10000);
370  mHMixNdC = new TH2D("Mixed_N_dC","Event Mixing: average-N vs delta-C",50,0,1500, 50,-100,100);
371  mHMixNC = new TH2D("Mixed_N_C","Event Mixing: average-N vs average-C",50,0,1500,50,0,10000);
372  mHMixNdN = new TH2D("Mixed_N_dN","Event Mixing: average-N vs delta-N",50,0,1500,50,-25,25);
373  mHMixdNdC = new TH2D("Mixed_dN_dC","Event Mixing: delta-N vs delta-C",50,-50,50, 50,-100,100);
374  mHMixdNC = new TH2D("Mixed_dN_C","Event Mixing: delta-N vs average-C",50,-50,50,50,0,10000);
375  mHMixCdC = new TH2D("Mixed_C_dC","Event Mixing: average-C vs delta-C",50,0,10000,50,-100,100);
376  mHcb = new TH2D("hcb","Cutbin usage",ncutbins,-0.5,ncutbins - 0.5,8,-0.5,7.5); // local
377  TH1::AddDirectory(kTRUE);
378 
379  mInit=true;
380 }
381 
382 //--------------------------------------------------------------------------
383 void StEStruct2ptCorrelations::finish(){
384 
385  if(!moutFileName){
386  cout<<" NO OUTPUTFILE TO WRITE TO ..... giving up ...."<<endl;
387  return;
388  }
389 
390  if(!mInit){
391  cout<<" WARNING: init=false"<<endl;
392  cout<<"No events were processed, either there was a problem reading input files or cuts are too restrictive"<<endl;
393  return;
394  }
395 
396  // We call finish explicitly in doEStruct, but it also gets called when
397  // we tell root to quit. I don't know the details of why. If we have
398  // already written histograms don't bother a second time.
399  // (Also I think finish() will invoke cleanUp() so second time here
400  // all the arrays will be gone and we get a core dump if we try
401  // touching them.)
402  if (mHistosWritten) {
403  return;
404  }
405  TH1::AddDirectory(kFALSE);
406  initHistograms();
407  fillHistograms();
408  TH1::AddDirectory(kTRUE);
409  TFile * tf=new TFile(moutFileName,"RECREATE");
410  tf->cd();
411  writeHistograms();
412  tf->Close();
413  // When we don't delete tf it gets reported as a memory leak
414  // (which it is, but we are quitting anyway so why bother with it?)
415  mHistosWritten = true;
416 }
417 
418 //--------------------------------------------------------------------------
419 void StEStruct2ptCorrelations::cleanUp(){
420  if(mDeleted) return;
421  deleteHistograms();
422  deleteArrays();
423 
424  int nZBins = 1;
425  if (mZBufferCutBinning) {
426  nZBins = kNumBuffers;
427  }
428  for (int iz=0;iz<nZBins;iz++) {
429  delete mHNEventsSib[iz];
430  delete mHNEventsMix[iz];
431  delete mHNEventsPosSib[iz];
432  delete mHNEventsPosMix[iz];
433  delete mHNEventsNegSib[iz];
434  delete mHNEventsNegMix[iz];
435  }
436  delete [] mHNEventsSib;
437  delete [] mHNEventsMix;
438  delete [] mHNEventsPosSib;
439  delete [] mHNEventsPosMix;
440  delete [] mHNEventsNegSib;
441  delete [] mHNEventsNegMix;
442  delete mHptAll;
443  mDeleted=true;
444 }
445 
446 //--------------------------------------------------------------------------
447 // Parse cuts file for limits on vertex position.
448 // Want 5cm wide buffers.
449 void StEStruct2ptCorrelations::setZBuffLimits(StEStructCuts* ecut) {
450 
451  kZBuffMin = ecut->minVal("primaryVertexZ");
452  kZBuffMax = ecut->maxVal("primaryVertexZ");
453  kNumBuffers = int( (kZBuffMax-kZBuffMin)/kBuffWidth );
454 
455  if(kZBuffMin==kZBuffMax){ // no z vertex cut ...
456  kZBuffMin=-100.;
457  kZBuffMax=100.;
458  kNumBuffers = 1;
459  kBuffWidth= (kZBuffMax - kZBuffMin) / kNumBuffers; // adjust widths
460  }
461 
462  cout<<"Setting ZBuffers: Max="<<kZBuffMax<<" Min="<<kZBuffMin<<" NumBuff="<<kNumBuffers<<endl;
463  if(kNumBuffers<=_MAX_ZVBINS_) return; // we're good to go
464 
465  kNumBuffers=_MAX_ZVBINS_;
466  kBuffWidth= (kZBuffMax - kZBuffMin) / kNumBuffers; // adjust widths
467  if(kBuffWidth>6.5) {
468  gMessMgr->Warning()<<" Zvertex Width ="<<kBuffWidth<<" gt 6.5 cm"<<endm;
469  }
470 }
471 //---------------------------------------------------------------------
472 void StEStruct2ptCorrelations::writeDiagnostics(){
473 
474  if (mOneZBuffer) {
475  return;
476  }
477  int nIn=0;
478  int nOut=0;
479  for(int i=0;i<kNumBuffers;i++){
480  nIn += mbuffer[i].numEventsIn();
481  nOut += mbuffer[i].numEventsDeleted();
482  }
483  cout<<" Analysis["<<analysisIndex()<<"] recieved "<<nIn<<" useful events, deleted "<<nOut<<" of them"<<endl;
484 
485 }
486 
487 //
488 //------- Analysis Function ------//
489 //
490 //--------------------------------------------------------------------------
491 bool StEStruct2ptCorrelations::doEvent(StEStructEvent* event){
492 
493  if(!event) {
494  mCurrentEvent=NULL;
495  return false;
496  }
497 
498  if(mInit == false) init();
499 
500  if(2>event->Ntrack()){
501  delete event;
502  return true;
503  }
504 
505  //cout << "Doing event " << event->EventID() << ": " << event->Ntrack() << " tracks." << endl;
506 
507  moveEvents();
508  mCurrentEvent=event;
509  int nZBin = 1;
510  int iZBin = 0;
511  if (mZBufferCutBinning) {
512  nZBin = kNumBuffers;
513  iZBin = bufferIndex();
514  }
515  mHNEventsSib[iZBin]->Fill(event->Ntrack());
516  mHNEventsPosSib[iZBin]->Fill(event->Npos());
517  mHNEventsNegSib[iZBin]->Fill(event->Nneg());
518 
519  // inclusive pt distribution
520  // Note that cuts are done in yt-yt space and bins are
521  // set up before pid is done (so assume pion mass.)
522  // I think these are just diagnostic histograms. Even so
523  // I am uncomfortable with this pt histogramming.
524  // (dEdx histogramming is ok though.)
525  // Assign mass according to dEdx pid here.
526  // djp 11/28/2005
527 
528  float Mass[] = {0.1396, 0.1396, 0.497, 0.9383};
529 
530  StEStructCutBin* cb = StEStructCutBin::Instance();
531 
532  StEStructTrackIterator iTotTrigger = 0;
533  double pTotTrigger = 0;
534  double yTotTrigger = 0;
535  double phiTotTrigger = 0;
536 
537  StEStructTrackCollection *tc;
538  for(int ich=0;ich<2;ich++){
539  if(ich==0){
540  tc=mCurrentEvent->TrackCollectionP();
541  } else {
542  tc=mCurrentEvent->TrackCollectionM();
543  }
544 
545  // Define event class by highest pt track.
546  // For each event class increment away, transverse and toward
547  // Within the ich loop we only have positive or negative. I had been adding plus and minus
548  // together thinking I was getting CD. With triggered stuff it doesn't work that way.
549  // Need one loop through to get trigger, another to get associated.
550  StEStructTrackIterator iTrigger = 0;
551  double pTrigger = 0;
552  double yTrigger = 0;
553  double phiTrigger = 0;
554  for(StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
555  if ((*iter)->Pt()>pTrigger) {
556  pTrigger = (*iter)->Pt();
557  yTrigger = (*iter)->Yt();
558  phiTrigger = (*iter)->Phi();
559  iTrigger = iter;
560  }
561  if ((*iter)->Pt()>pTotTrigger) {
562  pTotTrigger = (*iter)->Pt();
563  yTotTrigger = (*iter)->Yt();
564  phiTotTrigger = (*iter)->Phi();
565  iTotTrigger = iter;
566  }
567  }
568  if (pTrigger>5.0) {
569  pTrigger = 5.0;
570  }
571  if (yTrigger>4.5) {
572  yTrigger = 4.5;
573  }
574  double pi = acos(-1);
575  for(StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
576  double phi = (*iter)->Phi();
577  double dphi = phi - phiTrigger;
578  double pt = (*iter)->Pt();
579  double yt = (*iter)->Yt();
580  if (dphi<-pi) {
581  dphi += 2*pi;
582  } else if (dphi>pi) {
583  dphi -= 2*pi;
584  }
585  dphi = fabs(dphi);
586  if (iter == iTrigger) {
587  if (ich==0) {
588  mHPtTrigP->Fill(pTrigger);
589  mHYtTrigP->Fill(yTrigger);
590  } else if (ich==1) {
591  mHPtTrigM->Fill(pTrigger);
592  mHYtTrigM->Fill(yTrigger);
593  }
594  } else {
595  if (ich==0) {
596  mHPhiAssocP->Fill(dphi,yTrigger);
597  mHPhiAssocPtP->Fill(dphi,pTrigger);
598  if (dphi<pi/3) {
599  mHPtP[0]->Fill(pt,pTrigger);
600  mHYtP[0]->Fill(yt,yTrigger);
601  } else if (dphi<2*pi/3) {
602  mHPtP[1]->Fill(pt,pTrigger);
603  mHYtP[1]->Fill(yt,yTrigger);
604  } else if (dphi<5*pi/6) {
605  mHPtP[2]->Fill(pt,pTrigger);
606  mHYtP[2]->Fill(yt,yTrigger);
607  } else {
608  mHPtP[3]->Fill(pt,pTrigger);
609  mHYtP[3]->Fill(yt,yTrigger);
610  }
611  } else if (ich==1) {
612  mHPhiAssocM->Fill(dphi,yTrigger);
613  mHPhiAssocPtM->Fill(dphi,pTrigger);
614  if (dphi<pi/3) {
615  mHPtM[0]->Fill(pt,pTrigger);
616  mHYtM[0]->Fill(yt,yTrigger);
617  } else if (dphi<2*pi/3) {
618  mHPtM[1]->Fill(pt,pTrigger);
619  mHYtM[1]->Fill(yt,yTrigger);
620  } else if (dphi<5*pi/6) {
621  mHPtM[2]->Fill(pt,pTrigger);
622  mHYtM[2]->Fill(yt,yTrigger);
623  } else {
624  mHPtM[3]->Fill(pt,pTrigger);
625  mHYtM[3]->Fill(yt,yTrigger);
626  }
627  }
628  }
629  }
630 
631  for(StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
632  int parentClass = cb->getParentBin(mPairCuts,(*iter));
633  int ipb = parentClass*nZBin + iZBin;
634  mHMeanPtTot[ipb]->Fill((*iter)->Pt());
635  mHMeanYtTot[ipb]->Fill((*iter)->Yt());
636  mHEtaTot[ipb]->Fill((*iter)->Eta());
637  if (0 == ich) {
638  mHMeanPtP[ipb]->Fill((*iter)->Pt());
639  mHMeanYtP[ipb]->Fill((*iter)->Yt());
640  mHEtaP[ipb]->Fill((*iter)->Eta());
641  } else {
642  mHMeanPtM[ipb]->Fill((*iter)->Pt());
643  mHMeanYtM[ipb]->Fill((*iter)->Yt());
644  mHEtaM[ipb]->Fill((*iter)->Eta());
645  }
646  mQAHists->fillTrackHistograms(*iter,parentClass);
647 
648 if (0) {
649  // Choose mass according to dEdx (in which case transverse and longitudinal
650  // rapidities will be calculated as actual rapidities) or set mass to 0
651  // in which case we will use the quantity Jeff was using for transverse rapidity
652  // (a mid-rapidity approximation for pions) and eta for longitudinal rapidity.
653  // Sould really have a switch accessible from macro level instead of
654  // using cutMode.
655 
656  if (cb->getMode() == 5) {
657  // Try using real rapidity calculation.
658  (*iter)->SetMassAssignment(Mass[parentClass]);
659  // (*iter)->SetMassAssignment(0);
660  } else {
661  // 0 should be default, but just to be explicit here/
662  (*iter)->SetMassAssignment(0);
663  }
664  // Always use psuedo-rapidity for now.
665  // We are used to looking at the eta-phi plots (more or less.)
666  // (*iter)->SetMassAssignment(0);
667 } else {
668  (*iter)->SetMassAssignment(0);
669 }
670  }
671  }
672 
673  // Fill uncharged triggered histograms
674  for (int ich=0;ich<2;ich++) {
675  if (ich==0) {
676  tc=mCurrentEvent->TrackCollectionP();
677  } else {
678  tc=mCurrentEvent->TrackCollectionM();
679  }
680  if (pTotTrigger>5.0) {
681  pTotTrigger = 5.0;
682  }
683  if (yTotTrigger>4.5) {
684  yTotTrigger = 4.5;
685  }
686  double pi = acos(-1);
687  for (StEStructTrackIterator iter = tc->begin(); iter != tc->end(); iter++) {
688  double phi = (*iter)->Phi();
689  double dphi = phi - phiTotTrigger;
690  double pt = (*iter)->Pt();
691  double yt = (*iter)->Yt();
692  if (dphi<-pi) {
693  dphi += 2*pi;
694  } else if (dphi>pi) {
695  dphi -= 2*pi;
696  }
697  dphi = fabs(dphi);
698 
699  if (iter == iTotTrigger) {
700  mHPtTrigTot->Fill(pTotTrigger);
701  mHYtTrigTot->Fill(yTotTrigger);
702  } else {
703  mHPhiAssocTot->Fill(dphi,yTotTrigger);
704  mHPhiAssocPtTot->Fill(dphi,pTotTrigger);
705  if (dphi<pi/3) {
706  mHPtTot[0]->Fill(pt,pTotTrigger);
707  mHYtTot[0]->Fill(yt,yTotTrigger);
708  } else if (dphi<2*pi/3) {
709  mHPtTot[1]->Fill(pt,pTotTrigger);
710  mHYtTot[1]->Fill(yt,yTotTrigger);
711  } else if (dphi<5*pi/6) {
712  mHPtTot[2]->Fill(pt,pTotTrigger);
713  mHYtTot[2]->Fill(yt,yTotTrigger);
714  } else {
715  mHPtTot[3]->Fill(pt,pTotTrigger);
716  mHYtTot[3]->Fill(yt,yTotTrigger);
717  }
718  }
719  }
720  }
721 
722 
723  // Set magnetic field for pair cuts
724  mPairCuts->SetBField( mCurrentEvent->BField() );
725  return makeSiblingAndMixedPairs();
726 }
727 
728 
729 //-----------------------------------------------------------------------
730 int StEStruct2ptCorrelations::bufferIndex(){
731  if(!mCurrentEvent || kBuffWidth==0.) return -1;
732 
733  // this should only be in 1 place
734  return (int) floor((mCurrentEvent->VertexZ()-kZBuffMin)/kBuffWidth);
735 }
736 
737 //--------------------------------------------------------------------------
738 void StEStruct2ptCorrelations::moveEvents(){
739 
740  if(!mCurrentEvent) return;
741  if (mCurrentEvent->VertexZ() > kZBuffMax) {
742  return;
743  }
744 
745  if (mOneZBuffer) {
746  mOneZBuffer->addEvent(mCurrentEvent);
747  } else {
748  int i=bufferIndex();
749 
750  if(i<0 || i>kNumBuffers-1) return;
751  mbuffer[i].addEvent(mCurrentEvent);
752  // mbuffCounter[i]++;
753  }
754 
755 }
756 
757 //--------------------------------------------------------------------------
758 bool StEStruct2ptCorrelations::makeSiblingAndMixedPairs() {
759 
760  if(!mCurrentEvent) return false; // logic problem!
761  if (mCurrentEvent->VertexZ() > kZBuffMax) {
762  return false;
763  }
764 
765  int i=bufferIndex();
766 
767  if(i<0 || i>kNumBuffers-1) return false;
768 
769  // I finally realized the cases 2 and 6 could be handled
770  // by 1 and 5 (loop over -+ includes everything +- does.)
771  // Those cases should be dropped for a small gain in speed,
772  // but I would need to make more changes in the code than I
773  // want to right now. djp June 21, 2006
774  mInterestingPair = 0;
775  makePairs(mCurrentEvent,mCurrentEvent,0);
776  makePairs(mCurrentEvent,mCurrentEvent,1);
777  makePairs(mCurrentEvent,mCurrentEvent,2);
778  makePairs(mCurrentEvent,mCurrentEvent,3);
779 
780  if (mOneZBuffer) {
781  mOneZBuffer->resetCounter();
782  } else {
783  mbuffer[i].resetCounter();
784  }
785  int mult = mCurrentEvent->Ntrack();
786  //cout << "Event " << mCurrentEvent->EventID() << ": " << mult << " tracks; Z = " << mCurrentEvent->VertexZ() << ", i = " << i << endl;
787  //cout << i <<"\t"; mbuffer[i].Print();
788  while (1) {
789  if (mOneZBuffer) {
790  mMixingEvent = mOneZBuffer->nextEvent(mult,mCurrentEvent->VertexZ(),mCurrentEvent->ZDCCoincidence());
791  } else {
792  mMixingEvent = mbuffer[i].nextEvent(mult);
793  }
794  if (!mMixingEvent) break;
795  // Require magnetic field to have same sign.
796  // Changing BField sign seems to interchage + and -.
797  // Do we want to require same magnitude?
798  if (mCurrentEvent->BField()*mMixingEvent->BField() < 0) {
799  continue;
800  }
801  //cout << "\t mixing " << mMixingEvent->EventID() << ": " << mMixingEvent->Ntrack() << " tracks." << endl;
802  int iZBin = 0;
803  if (mZBufferCutBinning) {
804  iZBin = bufferIndex();
805  }
806  mHNEventsMix[iZBin]->Fill(mMixingEvent->Ntrack());
807  mHNEventsPosMix[iZBin]->Fill(mMixingEvent->Npos());
808  mHNEventsNegMix[iZBin]->Fill(mMixingEvent->Nneg());
809  float deltaZ = mCurrentEvent->VertexZ() - mMixingEvent->VertexZ();
810  float aveZ = (mCurrentEvent->VertexZ() + mMixingEvent->VertexZ())/2;
811  float deltaN = mCurrentEvent->Ntrack() - mMixingEvent->Ntrack();
812  float aveN = (mCurrentEvent->Ntrack() + mMixingEvent->Ntrack()) / 2;
813  float deltaC = mCurrentEvent->ZDCCoincidence() - mMixingEvent->ZDCCoincidence();
814  float aveC = (mCurrentEvent->ZDCCoincidence() + mMixingEvent->ZDCCoincidence())/2;
815  mHMixZdN->Fill(aveZ,deltaN);
816  mHMixZN->Fill(aveZ,aveN);
817  mHMixZdC->Fill(aveZ,deltaC);
818  mHMixZC->Fill(aveZ,aveC);
819  mHMixZdZ->Fill(aveZ,deltaZ);
820  mHMixdZdN->Fill(deltaZ,deltaN);
821  mHMixdZN->Fill(deltaZ,aveN);
822  mHMixdZdC->Fill(deltaZ,deltaC);
823  mHMixdZC->Fill(deltaZ,aveC);
824  mHMixNdC->Fill(aveN,deltaC);
825  mHMixNC->Fill(aveN,aveC);
826  mHMixNdN->Fill(aveN,deltaN);
827  mHMixdNdC->Fill(deltaN,deltaC);
828  mHMixdNC->Fill(deltaN,aveC);
829  mHMixCdC->Fill(aveC,deltaC);
830 
831  makePairs(mCurrentEvent,mMixingEvent,4);
832  makePairs(mCurrentEvent,mMixingEvent,5);
833  makePairs(mCurrentEvent,mMixingEvent,6);
834  makePairs(mMixingEvent,mCurrentEvent,5);
835  makePairs(mMixingEvent,mCurrentEvent,6);
836  makePairs(mCurrentEvent,mMixingEvent,7);
837  }
838 
839  return true;
840 }
841 
842 //--------------------------------------------------------------------------
843 /* For non-pid case:
844  * ++ put pair into array pp at +/- \delta\eta, +/- \delta\phi.
845  * -- put pair into array mm at +/- \delta\eta, +/- \delta\phi.
846  * +- put pair into array pm at +/- \delta\eta, +/- \delta\phi.
847  * -+ ignore pair.
848  * This is because we have an "extra" loop and the same pair will
849  * turn up with opposite sign \delta\eta.
850  * Always put pair into y1-y2 and y2-y1.
851  *
852  * For pid case we don't want to symmetrize yt-yt and we want to
853  * distinguish +- from -+ in the case of non-identical particles.
854  * ++ put pair into array pp at +/- \delta\eta, +/- \delta\phi.
855  * -- put pair into array mm at +/- \delta\eta, +/- \delta\phi.
856  * If pid1 == pid2
857  * +- put pair into array pm at +/- \delta\eta, +/- \delta\phi.
858  * -+ ignore
859  * put pair into
860  * If pid1 != pid2
861  * put pi+K-, pi+\bar p, K+\bar p pair into array pm at +/- \delta\eta, +/- \delta\phi.
862  * put pi-K+, pi-p, K-p pair into array mp at +/- \delta\eta, +/- \delta\phi.
863  *
864  */
865 
866 void StEStruct2ptCorrelations::makePairs(StEStructEvent* e1, StEStructEvent* e2, int j){
867 
868  double pi = acos(-1);
869 
870  if(j>=8) return;
871  StEStructTrackCollection* t1;
872  StEStructTrackCollection* t2;
873 
874  StEStructPairCuts& mPair = *mPairCuts;
875 
876  StEStructBinning* b = StEStructBinning::Instance();
877  StEStructCutBin* cb = StEStructCutBin::Instance();
878 
879  ytBins** ytyt;
880  ytBins** nytyt;
881  if (!mdontFillYtYt) {
882  ytyt = mYtYt[j];
883  nytyt = mNYtYt[j];
884  }
885 
886  dphiBins** jtdetadphi = mJtDEtaDPhi[j];
887  dphiBins** prjtdetadphi;
888  dphiBins** pajtdetadphi;
889  dphiBins** pbjtdetadphi;
890  if (!mdontFillMeanPt) {
891  prjtdetadphi = mPrJtDEtaDPhi[j];
892  pajtdetadphi = mPaJtDEtaDPhi[j];
893  pbjtdetadphi = mPbJtDEtaDPhi[j];
894  }
895 
896  etaBins** etaeta;
897  phiBins** phiphi;
898  phiBins** nphiphi;
899  etaBins** pretaeta;
900  etaBins** paetaeta;
901  etaBins** pbetaeta;
902  phiBins** prphiphi;
903  phiBins** paphiphi;
904  phiBins** pbphiphi;
905 
906  etaBins** etaetaSS;
907  etaBins** pretaetaSS;
908  etaBins** paetaetaSS;
909  etaBins** pbetaetaSS;
910  etaBins** etaetaAS;
911  etaBins** pretaetaAS;
912  etaBins** paetaetaAS;
913  etaBins** pbetaetaAS;
914 
915  if (mdoFillEtaEta) {
916  etaeta = mEtaEta[j];
917  if (mFillASSS) {
918  etaetaSS = mEtaEtaSS[j];
919  etaetaAS = mEtaEtaAS[j];
920  }
921  phiphi = mPhiPhi[j];
922  nphiphi = mNPhiPhi[j];
923  if (!mdontFillMeanPt) {
924  pretaeta = mPrEtaEta[j];
925  paetaeta = mPaEtaEta[j];
926  pbetaeta = mPbEtaEta[j];
927  if (mFillASSS) {
928  pretaetaSS = mPrEtaEtaSS[j];
929  paetaetaSS = mPaEtaEtaSS[j];
930  pbetaetaSS = mPbEtaEtaSS[j];
931  pretaetaAS = mPrEtaEtaAS[j];
932  paetaetaAS = mPaEtaEtaAS[j];
933  pbetaetaAS = mPbEtaEtaAS[j];
934  }
935  prphiphi = mPrPhiPhi[j];
936  paphiphi = mPaPhiPhi[j];
937  pbphiphi = mPbPhiPhi[j];
938  }
939  }
940 
941  dytBins** atytyt;
942  dytBins** atnytyt;
943  dphiBins** jtsetadphi;
944  dphiBins** jtnsetadphi;
945  dphiBins** prjtsetadphi;
946  dphiBins** pajtsetadphi;
947  dphiBins** pbjtsetadphi;
948  if (mdoFillSumHistograms) {
949  atytyt = mAtSYtDYt[j];
950  atnytyt = mAtNSYtDYt[j];
951  jtsetadphi = mJtSEtaDPhi[j];
952  jtnsetadphi = mJtNSEtaDPhi[j];
953  if (!mdontFillMeanPt) {
954  prjtsetadphi = mPrJtSEtaDPhi[j];
955  pajtsetadphi = mPaJtSEtaDPhi[j];
956  pbjtsetadphi = mPbJtSEtaDPhi[j];
957  }
958  }
959 
960  qBins* qinv;
961  qBins* nqinv;
962  if (mFillQInv) {
963  qinv = mQinv[j];
964  nqinv = mNQinv[j];
965  }
966 
967  TPCSepBins* avgtsep = mTPCAvgTSep[j];
968  TPCSepBins* avgzsep = mTPCAvgZSep[j];
969  TPCSepBins* enttsep = mTPCEntTSep[j];
970  TPCSepBins* entzsep = mTPCEntZSep[j];
971  TPCSepBins* midtsep = mTPCMidTSep[j];
972  TPCSepBins* midzsep = mTPCMidZSep[j];
973  TPCSepBins* exittsep = mTPCExitTSep[j];
974  TPCSepBins* exitzsep = mTPCExitZSep[j];
975  TPCSepBins* midtp = mTPCMidTdptP[j];
976  TPCSepBins* midtn = mTPCMidTdptN[j];
977  TPCSepBins* midzp = mTPCMidZdptP[j];
978  TPCSepBins* midzn = mTPCMidZdptN[j];
979  TPCSepBins* pairqual = mTPCQuality[j];
980 
981  TPCSepBins** avgtz = mTPCAvgTZ[j];
982  TPCSepBins** enttz = mTPCEntTZ[j];
983  TPCSepBins** midtz = mTPCMidTZ[j];
984  TPCSepBins** midtzc = mTPCMidTZC[j];
985  TPCSepBins** midtznc = mTPCMidTZNC[j];
986  TPCSepBins** exittz = mTPCExitTZ[j];
987  TPCSepBins** entqz = mTPCEntQZ[j];
988  TPCSepBins** midqz = mTPCMidQZ[j];
989  TPCSepBins** entqt = mTPCEntQT[j];
990  TPCSepBins** midqt = mTPCMidQT[j];
991  TPCSepBins** entqzt = mTPCEntQZT[j];
992  TPCSepBins** midqzt = mTPCMidQZT[j];
993  dptBins** enttd = mTPCEntTdpt[j];
994  dptBins** midtd = mTPCMidTdpt[j];
995  dptBins** exittd = mTPCExitTdpt[j];
996 
997 
998  switch(j) {
999  case 0:
1000  {
1001  t1=e1->TrackCollectionP();
1002  t2=e2->TrackCollectionP();
1003  mPair.setPairType(0);
1004  mpossiblePairs[j]+=(int)floor(0.5*(t1->getEntries()*(t2->getEntries()-1)));
1005  break;
1006  }
1007 
1008  case 1:
1009  {
1010  t1=e1->TrackCollectionP();
1011  t2=e2->TrackCollectionM();
1012  mPair.setPairType(1);
1013  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1014  break;
1015  }
1016  case 2:
1017  {
1018  t1=e1->TrackCollectionM();
1019  t2=e2->TrackCollectionP();
1020  mPair.setPairType(1);
1021  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1022  break;
1023  }
1024  case 3:
1025  {
1026  t1=e1->TrackCollectionM();
1027  t2=e2->TrackCollectionM();
1028  mPair.setPairType(0);
1029  mpossiblePairs[j]+=(int)floor(0.5*(t1->getEntries()*(t2->getEntries()-1)));
1030  break;
1031  }
1032  case 4:
1033  {
1034  t1=e1->TrackCollectionP();
1035  t2=e2->TrackCollectionP();
1036  mPair.setPairType(2);
1037  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1038  break;
1039  }
1040  case 5:
1041  {
1042  t1=e1->TrackCollectionP();
1043  t2=e2->TrackCollectionM();
1044  mPair.setPairType(3);
1045  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1046  break;
1047  }
1048  case 6:
1049  {
1050  t1=e1->TrackCollectionM();
1051  t2=e2->TrackCollectionP();
1052  mPair.setPairType(3);
1053  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1054  break;
1055  }
1056  case 7:
1057  {
1058  t1=e1->TrackCollectionM();
1059  t2=e2->TrackCollectionM();
1060  mPair.setPairType(2);
1061  mpossiblePairs[j]+=(int)(t1->getEntries()*t2->getEntries());
1062  break;
1063  }
1064  }
1065 
1066  // for pair cut analysis, in mixed events you need to correct for the offset from different primary vertex Z values
1067 
1068 // Why did we require mdoPairDensityHistograms to set the ZOffset????
1069 // if(j>=4 && mdoPairDensityHistograms) mPair.SetZoffset(e2->VertexZ() - e1->VertexZ());
1070 // else mPair.SetZoffset(0);
1071  // Note that for j < 4 e1 = e2 (so Vz2-Vz1 = 0)
1072  mPair.SetZoffset(e2->VertexZ() - e1->VertexZ());
1073 
1074  StEStructTrackIterator Iter1;
1075  StEStructTrackIterator Iter2;
1076 
1077  int iyt1,iyt2,idyt,isyt;
1078  int ipt1,ipt2;
1079  int ieta1,ieta2,ideta,iseta;
1080  int iphi1,iphi2,idphi;
1081  int isavgt, isavgz, isentt, isentz;
1082  int ismidt, ismidz, isexitt, isexitz, iqual;
1083  float pt1, pt2;
1084 
1085  if(mtimer)mtimer->start();
1086 
1087  int nZBin = 1;
1088  int iZBin = 0;
1089  if (mZBufferCutBinning) {
1090  nZBin = kNumBuffers;
1091  iZBin = bufferIndex();
1092  }
1093  // Should be extracting mass from pid information.
1094  // For now use mass=0 so we can compare eta with non-pid analysis.
1095  //>>>>>>>>>>>>
1096  float mass1 = 0, mass2 = 0;
1097 
1098  int it1 = -1;
1099  for(Iter1=t1->begin();Iter1!=t1->end();++Iter1){
1100  it1++;
1101 
1102  mPair.SetTrack1(*Iter1);
1103  int it2;
1104  if (j==0 || j==3) {
1105  Iter2 = Iter1+1;
1106  it2 = it1;
1107  } else {
1108  Iter2 = t2->begin();
1109  it2 = -1;
1110  }
1111 
1112  if (j == 0 || j==3) {
1113  mHptAll->Fill(mPair.Track1()->Pt());
1114  }
1115 
1116  for(; Iter2!=t2->end(); ++Iter2){
1117  it2++;
1118  numPairs[j]++;
1119  // QA pair histograms have not been implemented yet.
1120  // Looks like the mdoPairDensityHistograms sections do that job.
1121  // I think the idea of those was that they were turned on during initial data
1122  // check, then turned off to save CPU time.
1123  // mQAHists->fillPairHistograms(*Iter1,*Iter2,i,0);
1124  //
1125  // Pair density histograms are only being filled for pairs that pass pair cuts.
1126  // If you want to use pair density histograms to tune pair cuts you should
1127  // turn pair cuts off.
1128  mPair.SetTrack2(*Iter2);
1129  if( mskipPairCuts || (mPair.cutPair(mdoPairCutHistograms)==0) ){
1130  numPairsProcessed[j]++;
1131 
1132  int icb, jcb;
1133  if (!cb->ignorePair(&mPair) && ((jcb=cb->getCutBin(&mPair,j)) >= 0)) {
1134  int ncutbins=cb->getNumBins();
1135  mHcb->Fill(jcb,j); // now filled for all cb modes
1136  if(jcb>=ncutbins) {
1137  cout << "ERROR, got cutbin " << jcb << " of " << ncutbins << " possible." << endl;
1138  return;
1139  }
1140  icb = jcb*nZBin + iZBin;
1141 // mQAHists->fillPairHistograms(*Iter1,*Iter2,jcb,1);
1142 
1143  // The following is special code to help save interesting p-p events.
1144  // Meant to be used in run where we only look at low multiplicity
1145  // events. Set a flag indicating this event has an interesting pair.
1146  // Routine calling makePairs may want to save event.
1147 // if (2*5 == (icb & 2*5)) {
1148 // mInterestingPair++;;
1149 // }
1150 
1151  pt1 = mPair.Track1()->Pt();
1152  pt2 = mPair.Track2()->Pt();
1153 
1154  iyt1 = b->iyt(mPair.Track1()->Yt(mass1));
1155  ipt1 = b->ipt(pt1);
1156  ieta1 = b->ieta(mPair.Track1()->Eta(mass1));
1157  iphi1 = b->iphi(mPair.Track1()->Phi());
1158  iyt2 = b->iyt(mPair.Track2()->Yt(mass2));
1159  ipt2 = b->ipt(pt2);
1160  ieta2 = b->ieta(mPair.Track2()->Eta(mass2));
1161  iphi2 = b->iphi(mPair.Track2()->Phi());
1162 
1163  float pwgt=pt1*pt2;
1164  float spta=pt1;
1165  float sptb=pt2;
1166  float nwgt = 1.0;
1167  float ytwgt = 1.0;
1168  if( !mskipEtaDeltaWeight ) {
1169  ytwgt = b->getDEtaWeight(mPair.DeltaEta());
1170 // pwgt*=nwgt;
1171 // spta*=nwgt;
1172 // sptb*=nwgt;
1173  }
1174 
1175  /*
1176  * Makes no sense for switchXX _and_ symmetrizeXX to be true at same time.
1177  */
1178  int symmetrizeXX = cb->symmetrizeXX(&mPair);
1179  int switchXX = cb->switchXX(&mPair);
1180  if (switchXX) {
1181  if (!mdontFillYtYt) {
1182  ytyt[icb][iyt2].yt[iyt1]+=ytwgt;
1183  nytyt[icb][iyt2].yt[iyt1]+=1;
1184  }
1185 
1186  if (mdoFillEtaEta) {
1187  bool SS = fabs(mPair.DeltaPhi()) < pi/4 || fabs(fabs(mPair.DeltaPhi())-2*pi) < pi/4;
1188  bool AS = fabs(fabs(mPair.DeltaPhi())-pi) < pi/4;
1189  etaeta[icb][ieta2].eta[ieta1]+=1; // nwgt;
1190  phiphi[icb][iphi2].phi[iphi1]+=nwgt;
1191  nphiphi[icb][iphi2].phi[iphi1]+=1;
1192  if (mFillASSS) {
1193  if (SS) {
1194  etaetaSS[icb][ieta2].eta[ieta1]+=1;
1195  } else if (AS) {
1196  etaetaAS[icb][ieta2].eta[ieta1]+=1;
1197  }
1198  }
1199  if (!mdontFillMeanPt) {
1200  pretaeta[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1201  paetaeta[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1202  pbetaeta[icb][ieta2].eta[ieta1]+=spta/nwgt;
1203  prphiphi[icb][iphi2].phi[iphi1]+=pwgt;
1204  paphiphi[icb][iphi2].phi[iphi1]+=sptb;
1205  pbphiphi[icb][iphi2].phi[iphi1]+=spta;
1206  if (mFillASSS) {
1207  if (SS) {
1208  pretaetaSS[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1209  paetaetaSS[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1210  pbetaetaSS[icb][ieta2].eta[ieta1]+=spta/nwgt;
1211  } else if (AS) {
1212  pretaetaAS[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1213  paetaetaAS[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1214  pbetaetaAS[icb][ieta2].eta[ieta1]+=spta/nwgt;
1215  }
1216  }
1217  }
1218  }
1219  } else {
1220  if (!mdontFillYtYt) {
1221  ytyt[icb][iyt1].yt[iyt2]+=ytwgt;
1222  nytyt[icb][iyt1].yt[iyt2]+=1;
1223  }
1224 
1225  if (mdoFillEtaEta) {
1226  bool SS = fabs(mPair.DeltaPhi()) < pi/4 || fabs(fabs(mPair.DeltaPhi())-2*pi) < pi/4;
1227  bool AS = fabs(fabs(mPair.DeltaPhi())-pi) < pi/4;
1228  etaeta[icb][ieta1].eta[ieta2]+=1; // nwgt;
1229  phiphi[icb][iphi1].phi[iphi2]+=nwgt;
1230  nphiphi[icb][iphi1].phi[iphi2]+=1;
1231  if (mFillASSS) {
1232  if (SS) {
1233  etaetaSS[icb][ieta1].eta[ieta2]+=1;
1234  } else if (AS) {
1235  etaetaAS[icb][ieta1].eta[ieta2]+=1;
1236  }
1237  }
1238  if (!mdontFillMeanPt) {
1239  pretaeta[icb][ieta1].eta[ieta2]+=pwgt/nwgt;
1240  paetaeta[icb][ieta1].eta[ieta2]+=spta/nwgt;
1241  pbetaeta[icb][ieta1].eta[ieta2]+=sptb/nwgt;
1242  prphiphi[icb][iphi1].phi[iphi2]+=pwgt;
1243  paphiphi[icb][iphi1].phi[iphi2]+=spta;
1244  pbphiphi[icb][iphi1].phi[iphi2]+=sptb;
1245  if (mFillASSS) {
1246  if (SS) {
1247  pretaetaSS[icb][ieta1].eta[ieta2]+=pwgt/nwgt;
1248  paetaetaSS[icb][ieta1].eta[ieta2]+=sptb/nwgt;
1249  pbetaetaSS[icb][ieta1].eta[ieta2]+=spta/nwgt;
1250  } else if (AS) {
1251  pretaetaAS[icb][ieta1].eta[ieta2]+=pwgt/nwgt;
1252  paetaetaAS[icb][ieta1].eta[ieta2]+=sptb/nwgt;
1253  pbetaetaAS[icb][ieta1].eta[ieta2]+=spta/nwgt;
1254  }
1255  }
1256  }
1257  }
1258  if (symmetrizeXX) {
1259  if (!mdontFillYtYt) {
1260  ytyt[icb][iyt2].yt[iyt1]+=ytwgt;
1261  nytyt[icb][iyt2].yt[iyt1]+=1;
1262  }
1263 
1264  //-> X vs X (symmetry)
1265  if (mdoFillEtaEta) {
1266  bool SS = fabs(mPair.DeltaPhi()) < pi/4 || fabs(fabs(mPair.DeltaPhi())-2*pi) < pi/4;
1267  bool AS = fabs(fabs(mPair.DeltaPhi())-pi) < pi/4;
1268  etaeta[icb][ieta2].eta[ieta1]+=1; // nwgt;
1269  phiphi[icb][iphi2].phi[iphi1]+=nwgt;
1270  nphiphi[icb][iphi2].phi[iphi1]+=1;
1271  if (mFillASSS) {
1272  if (SS) {
1273  etaetaSS[icb][ieta2].eta[ieta1]+=1;
1274  } else if (AS) {
1275  etaetaAS[icb][ieta2].eta[ieta1]+=1;
1276  }
1277  }
1278  if (!mdontFillMeanPt) {
1279  pretaeta[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1280  paetaeta[icb][ieta2].eta[ieta1]+=spta/nwgt;
1281  pbetaeta[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1282  prphiphi[icb][iphi2].phi[iphi1]+=pwgt;
1283  paphiphi[icb][iphi2].phi[iphi1]+=spta;
1284  pbphiphi[icb][iphi2].phi[iphi1]+=sptb;
1285  if (mFillASSS) {
1286  if (SS) {
1287  pretaetaSS[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1288  paetaetaSS[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1289  pbetaetaSS[icb][ieta2].eta[ieta1]+=spta/nwgt;
1290  } else if (AS) {
1291  pretaetaAS[icb][ieta2].eta[ieta1]+=pwgt/nwgt;
1292  paetaetaAS[icb][ieta2].eta[ieta1]+=sptb/nwgt;
1293  pbetaetaAS[icb][ieta2].eta[ieta1]+=spta/nwgt;
1294  }
1295  }
1296  }
1297  }
1298  }
1299  }
1300 
1301  //-> delta y vs delta x
1302  ideta = b->ideta(mPair.DeltaEta(mass1,mass2));
1303  idphi = b->idphi(mPair.DeltaPhi());
1304 
1305  jtdetadphi[icb][ideta].dphi[idphi] +=nwgt;
1306  if (!mdontFillMeanPt) {
1307  prjtdetadphi[icb][ideta].dphi[idphi] += pwgt;
1308  if (switchXX) {
1309  pajtdetadphi[icb][ideta].dphi[idphi] += sptb;
1310  pbjtdetadphi[icb][ideta].dphi[idphi] += spta;
1311  } else {
1312  pajtdetadphi[icb][ideta].dphi[idphi] += spta;
1313  pbjtdetadphi[icb][ideta].dphi[idphi] += sptb;
1314  }
1315  }
1316 
1317  //-> Sum y vs delta x
1318  // For symmetry only reflect around the delta axis.
1319  // Symmetrization done later.
1320  if (mdoFillSumHistograms) {
1321  idyt = b->idyt(mPair.DeltaYt(mass1,mass2));
1322  isyt = b->isyt(mPair.SigmaYt(mass1,mass2));
1323  iseta= b->iseta(mPair.SigmaEta(mass1,mass2));
1324 
1325  atytyt[icb][isyt].dyt[idyt] +=ytwgt;
1326  atnytyt[icb][isyt].dyt[idyt] +=1;
1327 
1328  jtsetadphi[icb][iseta].dphi[idphi]+=nwgt;
1329  jtnsetadphi[icb][iseta].dphi[idphi]+=1;
1330  if (!mdontFillMeanPt) {
1331  prjtsetadphi[icb][iseta].dphi[idphi] += pwgt;
1332  if (switchXX) {
1333  pajtsetadphi[icb][iseta].dphi[idphi] += sptb;
1334  pbjtsetadphi[icb][iseta].dphi[idphi] += spta;
1335  } else {
1336  pajtsetadphi[icb][iseta].dphi[idphi] += spta;
1337  pbjtsetadphi[icb][iseta].dphi[idphi] += sptb;
1338  }
1339  }
1340  }
1341 
1342 
1343  if (mFillQInv) {
1344  qinv[icb].q[b->iq(mPair.qInv())]+=nwgt;
1345  nqinv[icb].q[b->iq(mPair.qInv())]+=1;
1346  }
1347 
1348  // Note that pair density histograms are filled within the pair cut check block.
1349  // Could move end of that block earlier and recalculate everything we need for
1350  // histograms (in case pair failed cuts).
1351  int jden = cb->getPairDensityBin(jcb);
1352  int ipcb = -1;
1353  if (jden >= 0) {
1354  ipcb = jden*nZBin + iZBin;
1355  }
1356  if (mdoPairDensityHistograms && ipcb >= 0) {
1357  avgtsep[ipcb].sep[isavgt=b->isep(mPair.NominalTpcAvgXYSeparation())]+=nwgt;
1358  avgzsep[ipcb].sep[isavgz=b->isep(mPair.NominalTpcAvgZSeparation())]+=nwgt;
1359  float entXYSep = mPair.NominalTpcXYEntranceSeparation();
1360  enttsep[ipcb].sep[isentt=b->isep(entXYSep)]+=nwgt;
1361  entzsep[ipcb].sep[isentz=b->isep(mPair.NominalTpcZEntranceSeparation())]+=nwgt;
1362  float midXYSep = mPair.MidTpcXYSeparation();
1363  float midZSep = mPair.MidTpcZSeparation();
1364  midtsep[ipcb].sep[ismidt=b->isep(midXYSep)]+=nwgt;
1365  midzsep[ipcb].sep[ismidz=b->isep(midZSep)]+=nwgt;
1366  float exitXYSep = mPair.NominalTpcXYExitSeparation();
1367  exittsep[ipcb].sep[isexitt=b->isep(exitXYSep)]+=nwgt;
1368  exitzsep[ipcb].sep[isexitz=b->isep(mPair.NominalTpcZExitSeparation())]+=nwgt;
1369  float qual = mPair.quality();
1370  pairqual[ipcb].sep[iqual=b->iqual(qual)]+=nwgt;
1371 
1372  avgtz[ipcb][isavgt].sep[isavgz]+=nwgt;
1373  enttz[ipcb][isentt].sep[isentz]+=nwgt;
1374  midtz[ipcb][ismidt].sep[ismidz]+=nwgt;
1375  exittz[ipcb][isexitt].sep[isexitz]+=nwgt;
1376  entqz[ipcb][iqual].sep[isentz]+=nwgt;
1377  midqz[ipcb][iqual].sep[ismidz]+=nwgt;
1378  entqt[ipcb][iqual].sep[isentt]+=nwgt;
1379  midqt[ipcb][iqual].sep[ismidt]+=nwgt;
1380  entqzt[ipcb][isentz].sep[isentt]+=qual;
1381  midqzt[ipcb][ismidz].sep[ismidt]+=qual;
1382 
1383  // need to rearrange pair so that deltaPhi>0, avoiding iSwitch and symmetrize for simplicity
1384  float delpt; // my delta pt
1385  int idelpt; // bin index
1386  double q1;
1387  if (mPair.Track1()->Phi() - mPair.Track2()->Phi() >= 0) {
1388  delpt = mPair.Track1()->Pt() - mPair.Track2()->Pt();
1389  q1 = mPair.mTrack1->Charge();
1390  } else { // redefine delta as 2-1
1391  delpt = mPair.Track2()->Pt() - mPair.Track1()->Pt();
1392  q1 = mPair.mTrack2->Charge();
1393  }
1394  idelpt = b->idpt(delpt);
1395  enttd[ipcb][isentt].dpt[idelpt]+=nwgt;
1396  midtd[ipcb][ismidt].dpt[idelpt]+=nwgt;
1397  exittd[ipcb][isexitt].dpt[idelpt]+=nwgt;
1398  if (1 == j || 2 == j || 5 == j || 6 == j) { // US
1399  if (q1>=0) {
1400  midtp[ipcb].sep[ismidt]+=nwgt;
1401  midzp[ipcb].sep[ismidz]+=nwgt;
1402  } else {
1403  midtn[ipcb].sep[ismidt]+=nwgt;
1404  midzn[ipcb].sep[ismidz]+=nwgt;
1405  }
1406  } else { // LS
1407  if (delpt>=0) {
1408  midtp[ipcb].sep[ismidt]+=nwgt;
1409  midzp[ipcb].sep[ismidz]+=nwgt;
1410  } else {
1411  midtn[ipcb].sep[ismidt]+=nwgt;
1412  midzn[ipcb].sep[ismidz]+=nwgt;
1413  }
1414  }
1415  pairqual[ipcb].sep[ismidz]+=nwgt;
1416 
1417  // If tracks cross then the sign of the transverse distance at the endpoints
1418  // must be opposite. Plot mid separation XY versus Z for the crossing
1419  // candidates and noc-crossing separately.
1420  StThreeVectorF ent1 = mPair.Track1()->NominalTpcEntrancePoint();
1421  StThreeVectorF ent2 = mPair.Track2()->NominalTpcEntrancePoint();
1422  float sinEnt = (ent1.x()*ent2.y() - ent1.y()*ent2.x());
1423  StThreeVectorF exit1 = mPair.Track1()->NominalTpcExitPoint();
1424  StThreeVectorF exit2 = mPair.Track2()->NominalTpcExitPoint();
1425  float sinExit = (exit1.x()*exit2.y() - exit1.y()*exit2.x());
1426  if (sinEnt*sinExit < 0) {
1427  midtzc[ipcb][ismidt].sep[ismidz]+=nwgt;
1428  } else {
1429  midtznc[ipcb][ismidt].sep[ismidz]+=nwgt;
1430  }
1431  } // pair density
1432  };// check on if we include pair.
1433  };// pair cut
1434  };// iter2 loop
1435  };// iter 1 loop
1436 
1437  if(mtimer)mtimer->stop();
1438 }
1439 
1440 
1441 //--------------------------------------------------------------------------
1442 void StEStruct2ptCorrelations::debug_CheckHistograms() {
1443  //===========================================
1444  // Debugging. Seems I am often filling pp Mode3 twice (at least DEtaDPhi histograms.)
1445  if (!mHcb) {
1446  return;
1447  }
1448  StEStructBinning* b=StEStructBinning::Instance();
1449  int numCutBins=StEStructCutBin::Instance()->getNumBins();
1450  int nzb = 1;
1451  if (mZBufferCutBinning) {
1452  nzb = kNumBuffers;
1453  }
1454  for (int i=0;i<8;i++) {
1455  dphiBins** jtdetadphi = mJtDEtaDPhi[i];
1456  int totDEtaDPhi = 0;
1457  int totCB = 0;
1458  for(int y=0;y<numCutBins;y++){
1459  for (int z=0;z<nzb;z++) {
1460  int yz = y*nzb + z;
1461  for(int k=0;k<b->detaBins();k++){
1462  for(int j=0;j<b->dphiBins();j++){
1463  totDEtaDPhi += jtdetadphi[yz][k].dphi[j];
1464  }
1465  }
1466  }
1467  totCB += mHcb->GetBinContent(y+1,i+1);
1468  }
1469  cout << "Integral of jtdetadphi = " << totDEtaDPhi << ", for i = " << i << ". Compare to hcb = " << totCB << endl;
1470  }
1471  //===========================================
1472 }
1473 //
1474 //------------ Below are init, delete, write functions -------///
1475 //
1476 
1477 //--------------------------------------------------------------------------
1478 void StEStruct2ptCorrelations::fillHistograms() {
1479 
1480  float xv,yv;
1481  StEStructBinning* b=StEStructBinning::Instance();
1482  int numCutBins=StEStructCutBin::Instance()->getNumBins();
1483  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
1484  int nzb = 1;
1485  if (mZBufferCutBinning) {
1486  nzb = kNumBuffers;
1487  }
1488 
1489  TH1::AddDirectory(kFALSE);
1490 
1491  // Only care about limits of mHEtaPhi, but if we put 1 into contents I guess we can keep track
1492  // of number of files added together.
1493  mHEtaPhi->SetBinContent(1,1,1);
1494 
1495  for(int i=0; i<8; i++){
1496 
1497  phiBins** nphiphi = mNPhiPhi[i];
1498  qBins* qinv;
1499  qBins* nqinv;
1500  if (!mFillQInv) {
1501  qinv = mQinv[i];
1502  nqinv = mNQinv[i];
1503  }
1504 
1505  ytBins** ytyt;
1506  ytBins** nytyt;
1507  if (!mdontFillYtYt) {
1508  ytyt = mYtYt[i];
1509  nytyt = mNYtYt[i];
1510  }
1511 
1512  dphiBins** jtdetadphi = mJtDEtaDPhi[i];
1513  dphiBins** prjtdetadphi;
1514  dphiBins** pajtdetadphi;
1515  dphiBins** pbjtdetadphi;
1516  if (!mdontFillMeanPt) {
1517  prjtdetadphi = mPrJtDEtaDPhi[i];
1518  pajtdetadphi = mPaJtDEtaDPhi[i];
1519  pbjtdetadphi = mPbJtDEtaDPhi[i];
1520  }
1521 
1522  etaBins** etaeta;
1523  phiBins** phiphi;
1524  etaBins** pretaeta;
1525  phiBins** prphiphi;
1526  etaBins** paetaeta;
1527  phiBins** paphiphi;
1528  etaBins** pbetaeta;
1529  phiBins** pbphiphi;
1530 
1531  etaBins** etaetaSS;
1532  etaBins** pretaetaSS;
1533  etaBins** paetaetaSS;
1534  etaBins** pbetaetaSS;
1535  etaBins** etaetaAS;
1536  etaBins** pretaetaAS;
1537  etaBins** paetaetaAS;
1538  etaBins** pbetaetaAS;
1539 
1540  if (mdoFillEtaEta) {
1541  etaeta = mEtaEta[i];
1542  phiphi = mPhiPhi[i];
1543  if (mFillASSS) {
1544  etaetaSS = mEtaEtaSS[i];
1545  etaetaAS = mEtaEtaAS[i];
1546  }
1547  if (!mdontFillMeanPt) {
1548  pretaeta = mPrEtaEta[i];
1549  prphiphi = mPrPhiPhi[i];
1550  paetaeta = mPaEtaEta[i];
1551  paphiphi = mPaPhiPhi[i];
1552  pbetaeta = mPbEtaEta[i];
1553  pbphiphi = mPbPhiPhi[i];
1554 
1555  if (mFillASSS) {
1556  pretaetaSS = mPrEtaEtaSS[i];
1557  paetaetaSS = mPaEtaEtaSS[i];
1558  pbetaetaSS = mPbEtaEtaSS[i];
1559  pretaetaAS = mPrEtaEtaAS[i];
1560  paetaetaAS = mPaEtaEtaAS[i];
1561  pbetaetaAS = mPbEtaEtaAS[i];
1562  }
1563  }
1564  }
1565 
1566  dytBins** atytyt;
1567  dytBins** atnytyt;
1568  dphiBins** jtsetadphi;
1569  dphiBins** jtnsetadphi;
1570  dphiBins** prjtsetadphi;
1571  dphiBins** pajtsetadphi;
1572  dphiBins** pbjtsetadphi;
1573  if (mdoFillSumHistograms) {
1574  atytyt = mAtSYtDYt[i];
1575  atnytyt = mAtNSYtDYt[i];
1576  jtsetadphi = mJtSEtaDPhi[i];
1577  jtnsetadphi = mJtNSEtaDPhi[i];
1578  if (!mdontFillMeanPt) {
1579  prjtsetadphi = mPrJtSEtaDPhi[i];
1580  pajtsetadphi = mPaJtSEtaDPhi[i];
1581  pbjtsetadphi = mPbJtSEtaDPhi[i];
1582  }
1583  }
1584 
1585 
1586  for(int y=0;y<numCutBins;y++){
1587  for (int z=0;z<nzb;z++) {
1588  int yz = y*nzb + z;
1589 
1590  if (!mdontFillYtYt) {
1591  createHist2D(mHYtYt,"YtYt",i,y,z,yz,b->ytBins(),b->ytMin(),b->ytMax(),b->ytBins(),b->ytMin(),b->ytMax());
1592  createHist2D(mHNYtYt,"NYtYt",i,y,z,yz,b->ytBins(),b->ytMin(),b->ytMax(),b->ytBins(),b->ytMin(),b->ytMax());
1593  for(int k=0;k<b->ytBins();k++){
1594  for(int j=0;j<b->ytBins();j++){
1595  mHYtYt[i][yz]->Fill(b->ytVal(k),b->ytVal(j),ytyt[yz][k].yt[j]);
1596  mHNYtYt[i][yz]->Fill(b->ytVal(k),b->ytVal(j),nytyt[yz][k].yt[j]);
1597  }
1598  }
1599  delete [] ytyt[yz];
1600  delete [] nytyt[yz];
1601  }
1602 
1603  createHist2D(mHJtDEtaDPhi, "DEtaDPhiArr",i,y,z,yz,b->detaBins(),0.5,b->detaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1604  for(int k=0;k<b->detaBins();k++){
1605  for(int j=0;j<b->dphiBins();j++){
1606  // Symmetrize dEta,dPhi in StEStructHAdd of Support code.
1607  // here is just a copy of the array.
1608  mHJtDEtaDPhi[i][yz]->SetBinContent(k+1,j+1,jtdetadphi[yz][k].dphi[j]);
1609  }
1610  }
1611  delete [] jtdetadphi[yz];
1612  if (!mdontFillMeanPt) {
1613  createHist2D(mHPrJtDEtaDPhi,"PrDEtaDPhiArr",i,y,z,yz,b->detaBins(),0.5,b->detaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1614  createHist2D(mHPaJtDEtaDPhi,"PaDEtaDPhiArr",i,y,z,yz,b->detaBins(),0.5,b->detaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1615  createHist2D(mHPbJtDEtaDPhi,"PbDEtaDPhiArr",i,y,z,yz,b->detaBins(),0.5,b->detaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1616  for(int k=0;k<b->detaBins();k++){
1617  for(int j=0;j<b->dphiBins();j++){
1618  // Symmetrize dEta,dPhi in StEStructHAdd of Support code.
1619  // here is just a copy of the array.
1620  mHPrJtDEtaDPhi[i][yz]->SetBinContent(k+1,j+1,prjtdetadphi[yz][k].dphi[j]);
1621  mHPaJtDEtaDPhi[i][yz]->SetBinContent(k+1,j+1,pajtdetadphi[yz][k].dphi[j]);
1622  mHPbJtDEtaDPhi[i][yz]->SetBinContent(k+1,j+1,pbjtdetadphi[yz][k].dphi[j]);
1623  }
1624  }
1625  delete [] prjtdetadphi[yz];
1626  delete [] pajtdetadphi[yz];
1627  delete [] pbjtdetadphi[yz];
1628  }
1629 
1630  if (mdoFillEtaEta) {
1631  createHist2D(mHPhiPhi,"PhiPhi",i,y,z,yz,b->phiBins(),b->phiMin(),b->phiMax(),b->phiBins(),b->phiMin(),b->phiMax());
1632  createHist2D(mHNPhiPhi,"NPhiPhi",i,y,z,yz,b->phiBins(),b->phiMin(),b->phiMax(),b->phiBins(),b->phiMin(),b->phiMax());
1633  for(int k=0;k<b->phiBins();k++){
1634  for(int j=0;j<b->phiBins();j++){
1635  mHPhiPhi[i][yz]->Fill(xv=b->phiVal(k),yv=b->phiVal(j),phiphi[yz][k].phi[j]);
1636  mHNPhiPhi[i][yz]->Fill(xv,yv,nphiphi[yz][k].phi[j]);
1637  }
1638  }
1639  delete [] phiphi[yz];
1640  delete [] nphiphi[yz];
1641  if (!mdontFillMeanPt) {
1642  createHist2D(mHPrPhiPhi,"PrPhiPhi",i,y,z,yz,b->phiBins(),b->phiMin(),b->phiMax(),b->phiBins(),b->phiMin(),b->phiMax());
1643  createHist2D(mHPaPhiPhi,"PaPhiPhi",i,y,z,yz,b->phiBins(),b->phiMin(),b->phiMax(),b->phiBins(),b->phiMin(),b->phiMax());
1644  createHist2D(mHPbPhiPhi,"PbPhiPhi",i,y,z,yz,b->phiBins(),b->phiMin(),b->phiMax(),b->phiBins(),b->phiMin(),b->phiMax());
1645  for(int k=0;k<b->phiBins();k++){
1646  for(int j=0;j<b->phiBins();j++){
1647  mHPrPhiPhi[i][yz]->Fill(xv=b->phiVal(k),yv=b->phiVal(j),prphiphi[yz][k].phi[j]);
1648  mHPaPhiPhi[i][yz]->Fill(xv,yv,paphiphi[yz][k].phi[j]);
1649  mHPbPhiPhi[i][yz]->Fill(xv,yv,pbphiphi[yz][k].phi[j]);
1650  }
1651  }
1652  delete [] prphiphi[yz];
1653  delete [] paphiphi[yz];
1654  delete [] pbphiphi[yz];
1655  }
1656 
1657  createHist2D(mHEtaEta,"EtaEta",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1658  for(int k=0;k<b->etaBins();k++){
1659  for(int j=0;j<b->etaBins();j++){
1660  mHEtaEta[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),etaeta[yz][k].eta[j]);
1661  }
1662  }
1663  delete [] etaeta[yz];
1664  if (!mdontFillMeanPt) {
1665  createHist2D(mHPrEtaEta,"PrEtaEta",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1666  createHist2D(mHPaEtaEta,"PaEtaEta",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1667  createHist2D(mHPbEtaEta,"PbEtaEta",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1668  for(int k=0;k<b->etaBins();k++){
1669  for(int j=0;j<b->etaBins();j++){
1670  mHPrEtaEta[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),pretaeta[yz][k].eta[j]);
1671  mHPaEtaEta[i][yz]->Fill(xv,yv,paetaeta[yz][k].eta[j]);
1672  mHPbEtaEta[i][yz]->Fill(xv,yv,pbetaeta[yz][k].eta[j]);
1673  }
1674  }
1675  delete [] pretaeta[yz];
1676  delete [] paetaeta[yz];
1677  delete [] pbetaeta[yz];
1678  }
1679 
1680  if (mFillASSS) {
1681  createHist2D(mHEtaEtaSS,"EtaEtaSS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1682  for(int k=0;k<b->etaBins();k++){
1683  for(int j=0;j<b->etaBins();j++){
1684  mHEtaEtaSS[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),etaetaSS[yz][k].eta[j]);
1685  }
1686  }
1687  delete [] etaetaSS[yz];
1688  if (!mdontFillMeanPt) {
1689  createHist2D(mHPrEtaEtaSS,"PrEtaEtaSS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1690  createHist2D(mHPaEtaEtaSS,"PaEtaEtaSS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1691  createHist2D(mHPbEtaEtaSS,"PbEtaEtaSS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1692  for(int k=0;k<b->etaBins();k++){
1693  for(int j=0;j<b->etaBins();j++){
1694  mHPrEtaEtaSS[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),pretaetaSS[yz][k].eta[j]);
1695  mHPaEtaEtaSS[i][yz]->Fill(xv,yv,paetaetaSS[yz][k].eta[j]);
1696  mHPbEtaEtaSS[i][yz]->Fill(xv,yv,pbetaetaSS[yz][k].eta[j]);
1697  }
1698  }
1699  delete [] pretaetaSS[yz];
1700  delete [] paetaetaSS[yz];
1701  delete [] pbetaetaSS[yz];
1702  }
1703  createHist2D(mHEtaEtaAS,"EtaEtaAS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1704  for(int k=0;k<b->etaBins();k++){
1705  for(int j=0;j<b->etaBins();j++){
1706  mHEtaEtaAS[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),etaetaAS[yz][k].eta[j]);
1707  }
1708  }
1709  delete [] etaetaAS[yz];
1710  if (!mdontFillMeanPt) {
1711  createHist2D(mHPrEtaEtaAS,"PrEtaEtaAS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1712  createHist2D(mHPaEtaEtaAS,"PaEtaEtaAS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1713  createHist2D(mHPbEtaEtaAS,"PbEtaEtaAS",i,y,z,yz,b->etaBins(),b->etaMin(),b->etaMax(),b->etaBins(),b->etaMin(),b->etaMax());
1714  for(int k=0;k<b->etaBins();k++){
1715  for(int j=0;j<b->etaBins();j++){
1716  mHPrEtaEtaAS[i][yz]->Fill(xv=b->etaVal(k),yv=b->etaVal(j),pretaetaAS[yz][k].eta[j]);
1717  mHPaEtaEtaAS[i][yz]->Fill(xv,yv,paetaetaAS[yz][k].eta[j]);
1718  mHPbEtaEtaAS[i][yz]->Fill(xv,yv,pbetaetaAS[yz][k].eta[j]);
1719  }
1720  }
1721  delete [] pretaetaAS[yz];
1722  delete [] paetaetaAS[yz];
1723  delete [] pbetaetaAS[yz];
1724  }
1725  }
1726  }
1727 
1728  if (mdoFillSumHistograms) {
1729  createHist2D(mHAtSYtDYt, "SYtDYt",i,y,z,yz,b->sytBins(),b->sytMin(),b->sytMax(),b->dytBins(),b->dytMin(),b->dytMax());
1730  createHist2D(mHAtNSYtDYt,"NSYtDYt",i,y,z,yz,b->sytBins(),b->sytMin(),b->sytMax(),b->dytBins(),b->dytMin(),b->dytMax());
1731  for(int k=0;k<b->sytBins();k++){
1732  for(int j=0;j<b->dytBins();j++){
1733  mHAtSYtDYt[i][yz]->Fill(b->sytVal(k),b->dytVal(j),atytyt[yz][k].dyt[j]);
1734  mHAtNSYtDYt[i][yz]->Fill(b->sytVal(k),b->dytVal(j),atnytyt[yz][k].dyt[j]);
1735  }
1736  }
1737  delete [] atytyt[yz];
1738  delete [] atnytyt[yz];
1739 
1740  createHist2D(mHJtSEtaDPhi, "SEtaDPhiArr",i,y,z,yz,b->setaBins(),0.5,b->setaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1741  createHist2D(mHJtNSEtaDPhi, "NSEtaDPhiArr",i,y,z,yz,b->setaBins(),0.5,b->setaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1742  for(int k=0;k<b->setaBins();k++) {
1743  for(int j=0;j<b->dphiBins();j++) {
1744  mHJtSEtaDPhi[i][yz]->SetBinContent(k+1,j+1,jtsetadphi[yz][k].dphi[j]);
1745  mHJtNSEtaDPhi[i][yz]->SetBinContent(k+1,j+1,jtnsetadphi[yz][k].dphi[j]);
1746  }
1747  }
1748  delete [] jtsetadphi[yz];
1749  delete [] jtnsetadphi[yz];
1750  if (!mdontFillMeanPt) {
1751  createHist2D(mHPrJtSEtaDPhi,"PrSEtaDPhiArr",i,y,z,yz,b->setaBins(),0.5,b->setaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1752  createHist2D(mHPaJtSEtaDPhi,"PaSEtaDPhiArr",i,y,z,yz,b->setaBins(),0.5,b->setaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1753  createHist2D(mHPbJtSEtaDPhi,"PbSEtaDPhiArr",i,y,z,yz,b->setaBins(),0.5,b->setaBins()+0.5,b->dphiBins(),0.5,b->dphiBins()+0.5);
1754  for(int k=0;k<b->setaBins();k++) {
1755  for(int j=0;j<b->dphiBins();j++) {
1756  mHPrJtSEtaDPhi[i][yz]->SetBinContent(k+1,j+1,prjtsetadphi[yz][k].dphi[j]);
1757  mHPaJtSEtaDPhi[i][yz]->SetBinContent(k+1,j+1,pajtsetadphi[yz][k].dphi[j]);
1758  mHPbJtSEtaDPhi[i][yz]->SetBinContent(k+1,j+1,pbjtsetadphi[yz][k].dphi[j]);
1759  }
1760  }
1761  delete [] prjtsetadphi[yz];
1762  delete [] pbjtsetadphi[yz];
1763  delete [] pajtsetadphi[yz];
1764  }
1765  }
1766 
1767  if (mFillQInv) {
1768  createHist1D(mHQinv,"Qinv",i,y,z,yz,b->qBins(),b->qMin(),b->qMax());
1769  createHist1D(mHNQinv,"NQinv",i,y,z,yz,b->qBins(),b->qMin(),b->qMax());
1770  for(int k=0;k<b->qBins();k++){
1771  mHQinv[i][yz]->Fill(b->qVal(k),qinv[yz].q[k]);
1772  mHNQinv[i][yz]->Fill(b->qVal(k),nqinv[yz].q[k]);
1773  }
1774  // delete [] qinv[yz];
1775  }
1776 
1777  } // for z
1778  } // for y
1779  } // for i
1780 
1781  if(mdoPairDensityHistograms) {
1782  for(int i=0; i<8; i++){
1783  TPCSepBins* avgtsep = mTPCAvgTSep[i];
1784  TPCSepBins* avgzsep = mTPCAvgZSep[i];
1785  TPCSepBins* enttsep = mTPCEntTSep[i];
1786  TPCSepBins* entzsep = mTPCEntZSep[i];
1787  TPCSepBins* midtsep = mTPCMidTSep[i];
1788  TPCSepBins* midzsep = mTPCMidZSep[i];
1789  TPCSepBins* exittsep= mTPCExitTSep[i];
1790  TPCSepBins* exitzsep= mTPCExitZSep[i];
1791  TPCSepBins* midtp = mTPCMidTdptP[i];
1792  TPCSepBins* midtn = mTPCMidTdptN[i];
1793  TPCSepBins* midzp = mTPCMidZdptP[i];
1794  TPCSepBins* midzn = mTPCMidZdptN[i];
1795  TPCSepBins* pairqual = mTPCQuality[i];
1796  TPCSepBins** avgtz = mTPCAvgTZ[i];
1797  TPCSepBins** enttz = mTPCEntTZ[i];
1798  TPCSepBins** midtz = mTPCMidTZ[i];
1799  TPCSepBins** midtzc = mTPCMidTZC[i];
1800  TPCSepBins** midtznc = mTPCMidTZNC[i];
1801  TPCSepBins** exittz = mTPCExitTZ[i];
1802  TPCSepBins** entqz = mTPCEntQZ[i];
1803  TPCSepBins** midqz = mTPCMidQZ[i];
1804  TPCSepBins** entqt = mTPCEntQT[i];
1805  TPCSepBins** midqt = mTPCMidQT[i];
1806  TPCSepBins** entqzt = mTPCEntQZT[i];
1807  TPCSepBins** midqzt = mTPCMidQZT[i];
1808  dptBins** enttd = mTPCEntTdpt[i];
1809  dptBins** midtd = mTPCMidTdpt[i];
1810  dptBins** exittd = mTPCExitTdpt[i];
1811 
1812  for(int y=0;y<nden;y++){
1813  for (int z=0;z<nzb;z++) {
1814  int yz = y*nzb + z;
1815 
1816  createHist1D(mHTPCAvgTSep,"TPCAvgTSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1817  createHist1D(mHTPCAvgZSep,"TPCAvgZSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1818  createHist1D(mHTPCEntTSep,"TPCEntTSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1819  createHist1D(mHTPCEntZSep,"TPCEntZSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1820  createHist1D(mHTPCMidTSep,"TPCMidTSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1821  createHist1D(mHTPCMidZSep,"TPCMidZSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1822  createHist1D(mHTPCExitTSep,"TPCExitTSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1823  createHist1D(mHTPCExitZSep,"TPCExitZSep",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1824  createHist1D(mHTPCMidTdptP,"TPCMidTdptP",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1825  createHist1D(mHTPCMidTdptN,"TPCMidTdptN",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1826  createHist1D(mHTPCMidZdptP,"TPCMidZdptP",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1827  createHist1D(mHTPCMidZdptN,"TPCMidZdptN",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1828  createHist1D(mHTPCQuality,"TPCQuality",i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax());
1829  for(int k=0;k<b->TPCSepBins();k++) {
1830  mHTPCAvgTSep[i][yz]->Fill(xv=b->sepVal(k),avgtsep[yz].sep[k]);
1831  mHTPCAvgZSep[i][yz]->Fill(xv,avgzsep[yz].sep[k]);
1832  mHTPCEntTSep[i][yz]->Fill(xv,enttsep[yz].sep[k]);
1833  mHTPCEntZSep[i][yz]->Fill(xv,entzsep[yz].sep[k]);
1834  mHTPCMidTSep[i][yz]->Fill(xv,midtsep[yz].sep[k]);
1835  mHTPCMidZSep[i][yz]->Fill(xv,midzsep[yz].sep[k]);
1836  mHTPCExitTSep[i][yz]->Fill(xv,exittsep[yz].sep[k]);
1837  mHTPCExitZSep[i][yz]->Fill(xv,exitzsep[yz].sep[k]);
1838  mHTPCMidTdptP[i][yz]->Fill(xv,midtp[yz].sep[k]);
1839  mHTPCMidTdptN[i][yz]->Fill(xv,midtn[yz].sep[k]);
1840  mHTPCMidZdptP[i][yz]->Fill(xv,midzp[yz].sep[k]);
1841  mHTPCMidZdptN[i][yz]->Fill(xv,midzn[yz].sep[k]);
1842  }
1843  createHist1D(mHTPCQuality,"TPCQuality",i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax());
1844  for(int k=0;k<b->TPCQualityBins();k++) {
1845  mHTPCQuality[i][yz]->Fill(b->qualityVal(k),pairqual[yz].sep[k]);
1846  }
1847 /*
1848  * I think this is the right way to delete these 1D arrays, but until I can test it
1849  * I will live with possibly using a little extra memory.
1850  delete avgtsep[yz];
1851  delete avgzsep[yz];
1852  delete enttsep[yz];
1853  delete entzsep[yz];
1854  delete midtsep[yz];
1855  delete midzsep[yz];
1856  delete exittsep[yz];
1857  delete exitzsep[yz];
1858  delete midtpsep[yz];
1859  delete midtnsep[yz];
1860  delete midzpsep[yz];
1861  delete midznsep[yz];
1862  */
1863  createHist2D(mHTPCAvgTZ, "TPCAvgTZ", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1864  createHist2D(mHTPCEntTZ, "TPCEntTZ", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1865  createHist2D(mHTPCMidTZ, "TPCMidTZ", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1866  createHist2D(mHTPCMidTZC, "TPCMidTZC", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1867  createHist2D(mHTPCMidTZNC, "TPCMidTZNC", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1868  createHist2D(mHTPCExitTZ,"TPCExitTZ",i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1869  for(int k=0;k<b->TPCSepBins();k++) {
1870  for(int j=0;j<b->TPCSepBins();j++) {
1871  mHTPCAvgTZ[i][yz]->Fill(xv=b->sepVal(k),yv=b->sepVal(j),avgtz[yz][k].sep[j]);
1872  mHTPCEntTZ[i][yz]->Fill(xv,yv,enttz[yz][k].sep[j]);
1873  mHTPCMidTZ[i][yz]->Fill(xv,yv,midtz[yz][k].sep[j]);
1874  mHTPCMidTZC[i][yz]->Fill(xv,yv,midtzc[yz][k].sep[j]);
1875  mHTPCMidTZNC[i][yz]->Fill(xv,yv,midtznc[yz][k].sep[j]);
1876  mHTPCExitTZ[i][yz]->Fill(xv,yv,exittz[yz][k].sep[j]);
1877  }
1878  }
1879  delete [] avgtz[yz];
1880  delete [] enttz[yz];
1881  delete [] midtz[yz];
1882  delete [] midtzc[yz];
1883  delete [] midtznc[yz];
1884  delete [] exittz[yz];
1885 
1886  createHist2D(mHTPCEntQZ,"TPCEntQZ", i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1887  createHist2D(mHTPCMidQZ,"TPCMidQZ", i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1888  createHist2D(mHTPCEntQT,"TPCEntQT", i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1889  createHist2D(mHTPCMidQT,"TPCMidQT", i,y,z,yz,b->TPCQualityBins(),b->TPCQualityMin(),b->TPCQualityMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1890  createHist2D(mHTPCEntQZT,"TPCEntQZT", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1891  createHist2D(mHTPCMidQZT,"TPCMidQZT", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax());
1892  for(int k=0;k<b->TPCQualityBins();k++) {
1893  for(int j=0;j<b->TPCSepBins();j++) {
1894  mHTPCEntQZ[i][yz]->Fill(xv=b->qualityVal(k),yv=b->sepVal(j),entqz[yz][k].sep[j]);
1895  mHTPCMidQZ[i][yz]->Fill(xv,yv,midqz[yz][k].sep[j]);
1896  mHTPCEntQT[i][yz]->Fill(xv,yv,entqt[yz][k].sep[j]);
1897  mHTPCMidQT[i][yz]->Fill(xv,yv,midqt[yz][k].sep[j]);
1898  mHTPCEntQZT[i][yz]->Fill(xv=b->sepVal(k),yv,entqzt[yz][k].sep[j]);
1899  mHTPCMidQZT[i][yz]->Fill(xv,yv,midqzt[yz][k].sep[j]);
1900  }
1901  }
1902  delete [] entqz[yz];
1903  delete [] midqz[yz];
1904  delete [] entqt[yz];
1905  delete [] midqt[yz];
1906  delete [] entqzt[yz];
1907  delete [] midqzt[yz];
1908 
1909  createHist2D(mHTPCEntTdpt, "TPCEntTdpt", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->dptBins(),b->dptMin(),b->dptMax());
1910  createHist2D(mHTPCMidTdpt, "TPCMidTdpt", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->dptBins(),b->dptMin(),b->dptMax());
1911  createHist2D(mHTPCExitTdpt, "TPCExitTdpt", i,y,z,yz,b->TPCSepBins(),b->TPCSepMin(),b->TPCSepMax(),b->dptBins(),b->dptMin(),b->dptMax());
1912  for(int k=0;k<b->TPCSepBins();k++) {
1913  for(int j=0;j<b->dptBins();j++) {
1914  mHTPCEntTdpt[i][yz]->Fill(xv=b->sepVal(k),yv=b->dptVal(j),enttd[yz][k].dpt[j]);
1915  mHTPCMidTdpt[i][yz]->Fill(xv,yv,midtd[yz][k].dpt[j]);
1916  mHTPCExitTdpt[i][yz]->Fill(xv,yv,exittd[yz][k].dpt[j]);
1917  }
1918  }
1919  delete [] enttd[yz];
1920  delete [] midtd[yz];
1921  delete [] exittd[yz];
1922  }
1923  }
1924  }
1925  } // if pair density
1926 
1927 }
1928 
1929 //--------------------------------------------------------------------------
1930 void StEStruct2ptCorrelations::writeHistograms() {
1931 
1932  int numCutBins=StEStructCutBin::Instance()->getNumBins();
1933  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
1934  int numParentBins=StEStructCutBin::Instance()->getNumParentBins();
1935  int nZBins = 1;
1936  if (mZBufferCutBinning) {
1937  numCutBins *= kNumBuffers;
1938  nden *= kNumBuffers;
1939  numParentBins *= kNumBuffers;
1940  nZBins = kNumBuffers;
1941  }
1942  mHEtaPhi->Write();
1943  for (int j=0;j<nZBins;j++) {
1944  mHNEventsSib[j]->Write();
1945  mHNEventsMix[j]->Write();
1946  mHNEventsPosSib[j]->Write();
1947  mHNEventsPosMix[j]->Write();
1948  mHNEventsNegSib[j]->Write();
1949  mHNEventsNegMix[j]->Write();
1950  }
1951  mHMixZdN->Write();
1952  mHMixZN->Write();
1953  mHMixZdC->Write();
1954  mHMixZC->Write();
1955  mHMixZdZ->Write();
1956  mHMixdZdN->Write();
1957  mHMixdZN->Write();
1958  mHMixdZdC->Write();
1959  mHMixdZC->Write();
1960  mHMixNdC->Write();
1961  mHMixNC->Write();
1962  mHMixNdN->Write();
1963  mHMixdNdC->Write();
1964  mHMixdNC->Write();
1965  mHMixCdC->Write();
1966  mHcb->Write();
1967 
1968  mHptAll->Write();
1969  for(int j=0;j<numParentBins;j++){
1970  mHMeanPtTot[j]->Write();
1971  mHMeanPtP[j]->Write();
1972  mHMeanPtM[j]->Write();
1973  mHMeanYtTot[j]->Write();
1974  mHMeanYtP[j]->Write();
1975  mHMeanYtM[j]->Write();
1976  mHEtaTot[j]->Write();
1977  mHEtaP[j]->Write();
1978  mHEtaM[j]->Write();
1979  }
1980  mHPtTot[0]->Write();
1981  mHPtP[0]->Write();
1982  mHPtM[0]->Write();
1983  mHPtTot[1]->Write();
1984  mHPtP[1]->Write();
1985  mHPtM[1]->Write();
1986  mHPtTot[2]->Write();
1987  mHPtP[2]->Write();
1988  mHPtM[2]->Write();
1989  mHPtTot[3]->Write();
1990  mHPtP[3]->Write();
1991  mHPtM[3]->Write();
1992  mHYtTot[0]->Write();
1993  mHYtP[0]->Write();
1994  mHYtM[0]->Write();
1995  mHYtTot[1]->Write();
1996  mHYtP[1]->Write();
1997  mHYtM[1]->Write();
1998  mHYtTot[2]->Write();
1999  mHYtP[2]->Write();
2000  mHYtM[2]->Write();
2001  mHYtTot[3]->Write();
2002  mHYtP[3]->Write();
2003  mHYtM[3]->Write();
2004  mHPhiAssocTot->Write();
2005  mHPhiAssocP->Write();
2006  mHPhiAssocM->Write();
2007  mHPhiAssocPtTot->Write();
2008  mHPhiAssocPtP->Write();
2009  mHPhiAssocPtM->Write();
2010  mHPtTrigTot->Write();
2011  mHPtTrigP->Write();
2012  mHPtTrigM->Write();
2013  mHYtTrigTot->Write();
2014  mHYtTrigP->Write();
2015  mHYtTrigM->Write();
2016  for(int i=0;i<8;i++){
2017  for(int j=0;j<numCutBins;j++){
2018  if (!mdontFillYtYt) {
2019  mHYtYt[i][j]->Write();
2020  mHNYtYt[i][j]->Write();
2021  }
2022 
2023  mHJtDEtaDPhi[i][j]->Write();
2024  if (!mdontFillMeanPt) {
2025  mHPrJtDEtaDPhi[i][j]->Write();
2026  mHPaJtDEtaDPhi[i][j]->Write();
2027  mHPbJtDEtaDPhi[i][j]->Write();
2028  }
2029 
2030  if (mdoFillEtaEta) {
2031  mHPhiPhi[i][j]->Write();
2032  mHNPhiPhi[i][j]->Write();
2033  mHEtaEta[i][j]->Write();
2034  if (mFillASSS) {
2035  mHEtaEtaSS[i][j]->Write();
2036  mHEtaEtaAS[i][j]->Write();
2037  }
2038  if (!mdontFillMeanPt) {
2039  mHPrPhiPhi[i][j]->Write();
2040  mHPrEtaEta[i][j]->Write();
2041  mHPaPhiPhi[i][j]->Write();
2042  mHPaEtaEta[i][j]->Write();
2043  mHPbPhiPhi[i][j]->Write();
2044  mHPbEtaEta[i][j]->Write();
2045  if (mFillASSS) {
2046  mHPrEtaEtaSS[i][j]->Write();
2047  mHPaEtaEtaSS[i][j]->Write();
2048  mHPbEtaEtaSS[i][j]->Write();
2049  mHPrEtaEtaAS[i][j]->Write();
2050  mHPaEtaEtaAS[i][j]->Write();
2051  mHPbEtaEtaAS[i][j]->Write();
2052  }
2053  }
2054  }
2055 
2056  if (mdoFillSumHistograms) {
2057  mHAtSYtDYt[i][j]->Write();
2058  mHAtNSYtDYt[i][j]->Write();
2059  mHJtSEtaDPhi[i][j]->Write();
2060  mHJtNSEtaDPhi[i][j]->Write();
2061  if (!mdontFillMeanPt) {
2062  mHPrJtSEtaDPhi[i][j]->Write();
2063  mHPaJtSEtaDPhi[i][j]->Write();
2064  mHPbJtSEtaDPhi[i][j]->Write();
2065  }
2066  }
2067 
2068  if (mFillQInv) {
2069  mHQinv[i][j]->Write();
2070  mHNQinv[i][j]->Write();
2071  }
2072  }
2073 
2074  if(mdoPairDensityHistograms) {
2075  for (int j=0;j<nden;j++) {
2076  mHTPCAvgTSep[i][j]->Write();
2077  mHTPCAvgZSep[i][j]->Write();
2078  mHTPCEntTSep[i][j]->Write();
2079  mHTPCEntZSep[i][j]->Write();
2080  mHTPCMidTSep[i][j]->Write();
2081  mHTPCMidZSep[i][j]->Write();
2082  mHTPCExitTSep[i][j]->Write();
2083  mHTPCExitZSep[i][j]->Write();
2084 
2085  mHTPCMidTdptP[i][j]->Write();
2086  mHTPCMidTdptN[i][j]->Write();
2087  mHTPCMidZdptP[i][j]->Write();
2088  mHTPCMidZdptN[i][j]->Write();
2089 
2090  mHTPCQuality[i][j]->Write();
2091 
2092  mHTPCAvgTZ[i][j]->Write();
2093  mHTPCEntTZ[i][j]->Write();
2094  mHTPCMidTZ[i][j]->Write();
2095  mHTPCMidTZC[i][j]->Write();
2096  mHTPCMidTZNC[i][j]->Write();
2097  mHTPCExitTZ[i][j]->Write();
2098  mHTPCEntQZ[i][j]->Write();
2099  mHTPCMidQZ[i][j]->Write();
2100  mHTPCEntQT[i][j]->Write();
2101  mHTPCMidQT[i][j]->Write();
2102  mHTPCEntQZT[i][j]->Write();
2103  mHTPCMidQZT[i][j]->Write();
2104  mHTPCEntTdpt[i][j]->Write();
2105  mHTPCMidTdpt[i][j]->Write();
2106  mHTPCExitTdpt[i][j]->Write();
2107  }
2108  }
2109  }
2110 }
2111 
2112 //--------------------------------------------------------------------------
2113 void StEStruct2ptCorrelations::writeQAHists(TFile *qtf) {
2114 
2115  if(!mlocalQAHists) return;
2116 
2117  if(!qtf){
2118  cout<<" NO QA OUTPUT TFile TO WRITE TO ... giving up ..."<<endl;
2119  return;
2120  }
2121 
2122  mQAHists->writeTrackHistograms(qtf);
2123 
2124 
2125 }
2126 
2127 //--------------------------------------------------------------------------
2128 void StEStruct2ptCorrelations::initArrays(){
2129 
2130  int numCutBins=StEStructCutBin::Instance()->getNumBins();
2131  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
2132  if (mZBufferCutBinning) {
2133  numCutBins *= kNumBuffers;
2134  nden *= kNumBuffers;
2135  }
2136 
2137  for(int i=0;i<8;i++){
2138 
2139  if (!mdontFillYtYt) {
2140  mYtYt[i]=new ytBins*[numCutBins];
2141  mNYtYt[i]=new ytBins*[numCutBins];
2142  }
2143 
2144  mJtDEtaDPhi[i] =new dphiBins*[numCutBins];
2145  if (!mdontFillMeanPt) {
2146  mPaJtDEtaDPhi[i]=new dphiBins*[numCutBins];
2147  mPbJtDEtaDPhi[i]=new dphiBins*[numCutBins];
2148  mPrJtDEtaDPhi[i]=new dphiBins*[numCutBins];
2149  }
2150 
2151  if (mdoFillEtaEta) {
2152  mEtaEta[i]=new etaBins*[numCutBins];
2153  mPhiPhi[i]=new phiBins*[numCutBins];
2154  mNPhiPhi[i]=new phiBins*[numCutBins];
2155  if (mFillASSS) {
2156  mEtaEtaSS[i]=new etaBins*[numCutBins];
2157  mEtaEtaAS[i]=new etaBins*[numCutBins];
2158  }
2159  if (!mdontFillMeanPt) {
2160  mPrEtaEta[i]=new etaBins*[numCutBins];
2161  mPaEtaEta[i]=new etaBins*[numCutBins];
2162  mPbEtaEta[i]=new etaBins*[numCutBins];
2163  mPrPhiPhi[i]=new phiBins*[numCutBins];
2164  mPaPhiPhi[i]=new phiBins*[numCutBins];
2165  mPbPhiPhi[i]=new phiBins*[numCutBins];
2166  if (mFillASSS) {
2167  mPrEtaEtaSS[i]=new etaBins*[numCutBins];
2168  mPaEtaEtaSS[i]=new etaBins*[numCutBins];
2169  mPbEtaEtaSS[i]=new etaBins*[numCutBins];
2170  mPrEtaEtaAS[i]=new etaBins*[numCutBins];
2171  mPaEtaEtaAS[i]=new etaBins*[numCutBins];
2172  mPbEtaEtaAS[i]=new etaBins*[numCutBins];
2173  }
2174  }
2175  }
2176 
2177  if (mdoFillSumHistograms) {
2178  mAtSYtDYt[i]=new dytBins*[numCutBins];
2179  mAtNSYtDYt[i]=new dytBins*[numCutBins];
2180  mJtSEtaDPhi[i]=new dphiBins*[numCutBins];
2181  mJtNSEtaDPhi[i]=new dphiBins*[numCutBins];
2182  if (!mdontFillMeanPt) {
2183  mPaJtSEtaDPhi[i]=new dphiBins*[numCutBins];
2184  mPbJtSEtaDPhi[i]=new dphiBins*[numCutBins];
2185  mPrJtSEtaDPhi[i]=new dphiBins*[numCutBins];
2186  }
2187  }
2188 
2189  /* --- I cut out the ql,qo,qs
2190  if(mPair.doHbt3D()){
2191  mQlQs[i]=new qBins*[numCutBins];
2192  mQoQop[i]=new qBins*[numCutBins];
2193  }
2194  */
2195  if (mFillQInv) {
2196  mQinv[i]=new qBins[numCutBins];
2197  mNQinv[i]=new qBins[numCutBins];
2198  memset(mQinv[i],0,numCutBins*sizeof(qBins)); // do the memset here
2199  memset(mNQinv[i],0,numCutBins*sizeof(qBins)); // do the memset here
2200  }
2201 
2202  if(mdoPairDensityHistograms) {
2203  mTPCAvgTSep[i]=new TPCSepBins[nden]; //1D
2204  mTPCAvgZSep[i]=new TPCSepBins[nden];
2205  mTPCEntTSep[i]=new TPCSepBins[nden];
2206  mTPCEntZSep[i]=new TPCSepBins[nden];
2207  mTPCMidTSep[i]=new TPCSepBins[nden];
2208  mTPCMidZSep[i]=new TPCSepBins[nden];
2209  mTPCExitTSep[i]=new TPCSepBins[nden];
2210  mTPCExitZSep[i]=new TPCSepBins[nden];
2211  mTPCMidTdptP[i]=new TPCSepBins[nden];
2212  mTPCMidTdptN[i]=new TPCSepBins[nden];
2213  mTPCMidZdptP[i]=new TPCSepBins[nden];
2214  mTPCMidZdptN[i]=new TPCSepBins[nden];
2215  memset(mTPCAvgTSep[i], 0,nden*sizeof(TPCSepBins));
2216  memset(mTPCAvgZSep[i], 0,nden*sizeof(TPCSepBins));
2217  memset(mTPCEntTSep[i], 0,nden*sizeof(TPCSepBins));
2218  memset(mTPCEntZSep[i], 0,nden*sizeof(TPCSepBins));
2219  memset(mTPCMidTSep[i], 0,nden*sizeof(TPCSepBins));
2220  memset(mTPCMidZSep[i], 0,nden*sizeof(TPCSepBins));
2221  memset(mTPCExitTSep[i], 0,nden*sizeof(TPCSepBins));
2222  memset(mTPCExitZSep[i], 0,nden*sizeof(TPCSepBins));
2223  memset(mTPCMidTdptP[i],0,nden*sizeof(TPCSepBins));
2224  memset(mTPCMidTdptN[i],0,nden*sizeof(TPCSepBins));
2225  memset(mTPCMidZdptP[i],0,nden*sizeof(TPCSepBins));
2226  memset(mTPCMidZdptN[i],0,nden*sizeof(TPCSepBins));
2227  mTPCQuality[i]=new TPCSepBins[nden];
2228  memset(mTPCQuality[i], 0,nden*sizeof(TPCSepBins));
2229 
2230  mTPCAvgTZ[i]=new TPCSepBins*[nden]; //2D
2231  mTPCEntTZ[i]=new TPCSepBins*[nden];
2232  mTPCMidTZ[i]=new TPCSepBins*[nden];
2233  mTPCMidTZC[i]=new TPCSepBins*[nden];
2234  mTPCMidTZNC[i]=new TPCSepBins*[nden];
2235  mTPCExitTZ[i]=new TPCSepBins*[nden];
2236  mTPCEntQZ[i]=new TPCSepBins*[nden];
2237  mTPCMidQZ[i]=new TPCSepBins*[nden];
2238  mTPCEntQT[i]=new TPCSepBins*[nden];
2239  mTPCMidQT[i]=new TPCSepBins*[nden];
2240  mTPCEntQZT[i]=new TPCSepBins*[nden];
2241  mTPCMidQZT[i]=new TPCSepBins*[nden];
2242  mTPCEntTdpt[i]=new dptBins*[nden];
2243  mTPCMidTdpt[i]=new dptBins*[nden];
2244  mTPCExitTdpt[i]=new dptBins*[nden]; // Initialization to 0 done later.
2245  }
2246 
2247 
2248  for(int j=0;j<numCutBins;j++){
2249  if (!mdontFillYtYt) {
2250  mYtYt[i][j]=new ytBins[ESTRUCT_YT_BINS];
2251  memset(mYtYt[i][j],0,ESTRUCT_YT_BINS*sizeof(ytBins));
2252  mNYtYt[i][j]=new ytBins[ESTRUCT_YT_BINS];
2253  memset(mNYtYt[i][j],0,ESTRUCT_YT_BINS*sizeof(ytBins));
2254  }
2255 
2256  mJtDEtaDPhi[i][j]=new dphiBins[ESTRUCT_DETA_BINS];
2257  memset(mJtDEtaDPhi[i][j],0,ESTRUCT_DETA_BINS*sizeof(dphiBins));
2258  if (!mdontFillMeanPt) {
2259  mPaJtDEtaDPhi[i][j]=new dphiBins[ESTRUCT_DETA_BINS];
2260  memset(mPaJtDEtaDPhi[i][j],0,ESTRUCT_DETA_BINS*sizeof(dphiBins));
2261  mPbJtDEtaDPhi[i][j]=new dphiBins[ESTRUCT_DETA_BINS];
2262  memset(mPbJtDEtaDPhi[i][j],0,ESTRUCT_DETA_BINS*sizeof(dphiBins));
2263  mPrJtDEtaDPhi[i][j]=new dphiBins[ESTRUCT_DETA_BINS];
2264  memset(mPrJtDEtaDPhi[i][j],0,ESTRUCT_DETA_BINS*sizeof(dphiBins));
2265  }
2266 
2267  if (mdoFillEtaEta) {
2268  mEtaEta[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2269  memset(mEtaEta[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2270  mPhiPhi[i][j]=new phiBins[ESTRUCT_PHI_BINS];
2271  memset(mPhiPhi[i][j],0,ESTRUCT_PHI_BINS*sizeof(phiBins));
2272  mNPhiPhi[i][j]=new phiBins[ESTRUCT_PHI_BINS];
2273  memset(mNPhiPhi[i][j],0,ESTRUCT_PHI_BINS*sizeof(phiBins));
2274  if (mFillASSS) {
2275  mEtaEtaSS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2276  memset(mEtaEtaSS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2277  mEtaEtaAS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2278  memset(mEtaEtaAS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2279  }
2280  if (!mdontFillMeanPt) {
2281  mPrEtaEta[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2282  memset(mPrEtaEta[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2283  mPaEtaEta[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2284  memset(mPaEtaEta[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2285  mPbEtaEta[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2286  memset(mPbEtaEta[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2287  mPrPhiPhi[i][j]=new phiBins[ESTRUCT_PHI_BINS];
2288  memset(mPrPhiPhi[i][j],0,ESTRUCT_PHI_BINS*sizeof(phiBins));
2289  mPaPhiPhi[i][j]=new phiBins[ESTRUCT_PHI_BINS];
2290  memset(mPaPhiPhi[i][j],0,ESTRUCT_PHI_BINS*sizeof(phiBins));
2291  mPbPhiPhi[i][j]=new phiBins[ESTRUCT_PHI_BINS];
2292  memset(mPbPhiPhi[i][j],0,ESTRUCT_PHI_BINS*sizeof(phiBins));
2293  if (mFillASSS) {
2294  mPrEtaEtaSS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2295  memset(mPrEtaEtaSS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2296  mPaEtaEtaSS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2297  memset(mPaEtaEtaSS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2298  mPbEtaEtaSS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2299  memset(mPbEtaEtaSS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2300  mPrEtaEtaAS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2301  memset(mPrEtaEtaAS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2302  mPaEtaEtaAS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2303  memset(mPaEtaEtaAS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2304  mPbEtaEtaAS[i][j]=new etaBins[ESTRUCT_ETA_BINS];
2305  memset(mPbEtaEtaAS[i][j],0,ESTRUCT_ETA_BINS*sizeof(etaBins));
2306  }
2307  }
2308  }
2309 
2310  if (mdoFillSumHistograms) {
2311  mAtSYtDYt[i][j]=new dytBins[ESTRUCT_SYT_BINS];
2312  memset(mAtSYtDYt[i][j],0,ESTRUCT_SYT_BINS*sizeof(dytBins));
2313  mAtNSYtDYt[i][j]=new dytBins[ESTRUCT_SYT_BINS];
2314  memset(mAtNSYtDYt[i][j],0,ESTRUCT_SYT_BINS*sizeof(dytBins));
2315  mJtSEtaDPhi[i][j]=new dphiBins[ESTRUCT_SETA_BINS];
2316  memset(mJtSEtaDPhi[i][j],0,ESTRUCT_SETA_BINS*sizeof(dphiBins));
2317  mJtNSEtaDPhi[i][j]=new dphiBins[ESTRUCT_SETA_BINS];
2318  memset(mJtNSEtaDPhi[i][j],0,ESTRUCT_SETA_BINS*sizeof(dphiBins));
2319  if (!mdontFillMeanPt) {
2320  mPrJtSEtaDPhi[i][j]=new dphiBins[ESTRUCT_SETA_BINS];
2321  memset(mPrJtSEtaDPhi[i][j],0,ESTRUCT_SETA_BINS*sizeof(dphiBins));
2322  mPaJtSEtaDPhi[i][j]=new dphiBins[ESTRUCT_SETA_BINS];
2323  memset(mPaJtSEtaDPhi[i][j],0,ESTRUCT_SETA_BINS*sizeof(dphiBins));
2324  mPbJtSEtaDPhi[i][j]=new dphiBins[ESTRUCT_SETA_BINS];
2325  memset(mPbJtSEtaDPhi[i][j],0,ESTRUCT_SETA_BINS*sizeof(dphiBins));
2326  }
2327  }
2328  }
2329 
2330  if(mdoPairDensityHistograms) {
2331  for(int j=0;j<nden;j++){
2332  mTPCAvgTZ[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2333  memset(mTPCAvgTZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2334  mTPCEntTZ[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2335  memset(mTPCEntTZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2336  mTPCMidTZ[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2337  memset(mTPCMidTZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2338  mTPCMidTZC[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2339  memset(mTPCMidTZC[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2340  mTPCMidTZNC[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2341  memset(mTPCMidTZNC[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2342  mTPCExitTZ[i][j]=new TPCSepBins[ESTRUCT_TPCSEP_BINS];
2343  memset(mTPCExitTZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2344  mTPCEntQZ[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2345  memset(mTPCEntQZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2346  mTPCMidQZ[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2347  memset(mTPCMidQZ[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2348  mTPCEntQT[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2349  memset(mTPCEntQT[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2350  mTPCMidQT[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2351  memset(mTPCMidQT[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2352  mTPCEntQZT[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2353  memset(mTPCEntQZT[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2354  mTPCMidQZT[i][j]=new TPCSepBins[ESTRUCT_TPCQUALITY_BINS];
2355  memset(mTPCMidQZT[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(TPCSepBins));
2356 
2357  mTPCEntTdpt[i][j]=new dptBins[ESTRUCT_TPCSEP_BINS];
2358  memset(mTPCEntTdpt[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(dptBins));
2359  mTPCMidTdpt[i][j]=new dptBins[ESTRUCT_TPCSEP_BINS];
2360  memset(mTPCMidTdpt[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(dptBins));
2361  mTPCExitTdpt[i][j]=new dptBins[ESTRUCT_TPCSEP_BINS];
2362  memset(mTPCExitTdpt[i][j],0,ESTRUCT_TPCSEP_BINS*sizeof(dptBins));
2363  }
2364 
2365  }
2366  }
2367 }
2368 
2369 //--------------------------------------------------------------------------
2370 void StEStruct2ptCorrelations::deleteArrays(){
2371 
2372  int numCutBins=StEStructCutBin::Instance()->getNumBins();
2373  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
2374  if (mZBufferCutBinning) {
2375  numCutBins *= kNumBuffers;
2376  nden *= kNumBuffers;
2377  }
2378 
2379  for(int i=0;i<8;i++){
2380  for(int j=0;j<numCutBins;j++){ // Why were these all commented out?
2381 
2382  if (!mdontFillYtYt) {
2383  delete [] mYtYt[i][j];
2384  delete [] mNYtYt[i][j];
2385  }
2386 
2387  delete [] mJtDEtaDPhi[i][j];
2388  if (!mdontFillMeanPt) {
2389  delete [] mPrJtDEtaDPhi[i][j];
2390  delete [] mPaJtDEtaDPhi[i][j];
2391  delete [] mPbJtDEtaDPhi[i][j];
2392  }
2393 
2394  if (mdoFillEtaEta) {
2395  delete [] mEtaEta[i][j];
2396  delete [] mPhiPhi[i][j];
2397  delete [] mNPhiPhi[i][j];
2398  if (mFillASSS) {
2399  delete [] mEtaEtaSS[i][j];
2400  delete [] mEtaEtaAS[i][j];
2401  }
2402  if (!mdontFillMeanPt) {
2403  delete [] mPrEtaEta[i][j];
2404  delete [] mPrPhiPhi[i][j];
2405  delete [] mPaEtaEta[i][j];
2406  delete [] mPaPhiPhi[i][j];
2407  delete [] mPbEtaEta[i][j];
2408  delete [] mPbPhiPhi[i][j];
2409  if (mFillASSS) {
2410  delete [] mPrEtaEtaSS[i][j];
2411  delete [] mPaEtaEtaSS[i][j];
2412  delete [] mPbEtaEtaSS[i][j];
2413  delete [] mPrEtaEtaAS[i][j];
2414  delete [] mPaEtaEtaAS[i][j];
2415  delete [] mPbEtaEtaAS[i][j];
2416  }
2417  }
2418  }
2419 
2420  if (mdoFillSumHistograms) {
2421  delete [] mAtSYtDYt[i][j];
2422  delete [] mJtSEtaDPhi[i][j];
2423  delete [] mJtNSEtaDPhi[i][j];
2424  if (!mdontFillMeanPt) {
2425  delete [] mPrJtSEtaDPhi[i][j];
2426  delete [] mPaJtSEtaDPhi[i][j];
2427  delete [] mPbJtSEtaDPhi[i][j];
2428  }
2429  }
2430 
2431  }
2432 
2433  if(mdoPairDensityHistograms) { // These are the 2D arrays. Delete along one axis.
2434  for(int j=0;j<nden;j++){
2435  delete [] mTPCAvgTZ[i][j];
2436  delete [] mTPCEntTZ[i][j];
2437  delete [] mTPCMidTZ[i][j];
2438  delete [] mTPCMidTZC[i][j];
2439  delete [] mTPCMidTZNC[i][j];
2440  delete [] mTPCExitTZ[i][j];
2441  delete [] mTPCEntQZ[i][j];
2442  delete [] mTPCMidQZ[i][j];
2443  delete [] mTPCEntQT[i][j];
2444  delete [] mTPCMidQT[i][j];
2445  delete [] mTPCEntQZT[i][j];
2446  delete [] mTPCMidQZT[i][j];
2447  delete [] mTPCEntTdpt[i][j];
2448  delete [] mTPCMidTdpt[i][j];
2449  delete [] mTPCExitTdpt[i][j];
2450  }
2451  }
2452 
2453  if (!mdontFillYtYt) {
2454  delete [] mYtYt[i];
2455  delete [] mNYtYt[i];
2456  }
2457 
2458  delete [] mJtDEtaDPhi[i];
2459  if (!mdontFillMeanPt) {
2460  delete [] mPrJtDEtaDPhi[i];
2461  delete [] mPaJtDEtaDPhi[i];
2462  delete [] mPbJtDEtaDPhi[i];
2463  }
2464 
2465  if (mdoFillEtaEta) {
2466  delete [] mEtaEta[i];
2467  delete [] mPhiPhi[i];
2468  delete [] mNPhiPhi[i];
2469  if (mFillASSS) {
2470  delete [] mEtaEtaSS[i];
2471  delete [] mEtaEtaAS[i];
2472  }
2473  if (!mdontFillMeanPt) {
2474  delete [] mPrEtaEta[i];
2475  delete [] mPrPhiPhi[i];
2476  delete [] mPaEtaEta[i];
2477  delete [] mPaPhiPhi[i];
2478  delete [] mPbEtaEta[i];
2479  delete [] mPbPhiPhi[i];
2480  if (mFillASSS) {
2481  delete [] mPrEtaEtaSS[i];
2482  delete [] mPaEtaEtaSS[i];
2483  delete [] mPbEtaEtaSS[i];
2484  delete [] mPrEtaEtaAS[i];
2485  delete [] mPaEtaEtaAS[i];
2486  delete [] mPbEtaEtaAS[i];
2487  }
2488  }
2489  }
2490 
2491  if (mdoFillSumHistograms) {
2492  delete [] mAtSYtDYt[i];
2493  delete [] mJtSEtaDPhi[i];
2494  delete [] mJtNSEtaDPhi[i];
2495  if (!mdontFillMeanPt) {
2496  delete [] mPrJtSEtaDPhi[i];
2497  delete [] mPaJtSEtaDPhi[i];
2498  delete [] mPbJtSEtaDPhi[i];
2499  }
2500  }
2501 
2502  if (mFillQInv) {
2503  delete [] mQinv[i];
2504  delete [] mNQinv[i];
2505  }
2506 
2507  if(mdoPairDensityHistograms) {
2508  delete [] mTPCAvgTSep[i];
2509  delete [] mTPCAvgZSep[i];
2510  delete [] mTPCEntTSep[i];
2511  delete [] mTPCEntZSep[i];
2512  delete [] mTPCMidTSep[i];
2513  delete [] mTPCMidZSep[i];
2514  delete [] mTPCExitTSep[i];
2515  delete [] mTPCExitZSep[i];
2516  delete [] mTPCMidTdptP[i];
2517  delete [] mTPCMidTdptN[i];
2518  delete [] mTPCMidZdptP[i];
2519  delete [] mTPCMidZdptN[i];
2520  delete [] mTPCQuality[i];
2521  delete [] mTPCAvgTZ[i];
2522  delete [] mTPCEntTZ[i];
2523  delete [] mTPCMidTZ[i];
2524  delete [] mTPCMidTZC[i];
2525  delete [] mTPCMidTZNC[i];
2526  delete [] mTPCExitTZ[i];
2527  delete [] mTPCEntQZ[i];
2528  delete [] mTPCMidQZ[i];
2529  delete [] mTPCEntQT[i];
2530  delete [] mTPCMidQT[i];
2531  delete [] mTPCEntQZT[i];
2532  delete [] mTPCMidQZT[i];
2533  delete [] mTPCEntTdpt[i];
2534  delete [] mTPCMidTdpt[i];
2535  delete [] mTPCExitTdpt[i];
2536  }
2537  }
2538 
2539 }
2540 
2541 //--------------------------------------------------------------------------
2542 
2543 void StEStruct2ptCorrelations::initHistograms(){
2544 
2545  int numCutBins=StEStructCutBin::Instance()->getNumBins();
2546  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
2547  int numParentBins=StEStructCutBin::Instance()->getNumParentBins();
2548  if (mZBufferCutBinning) {
2549  numCutBins *= kNumBuffers;
2550  nden *= kNumBuffers;
2551  numParentBins *= kNumBuffers;
2552  }
2553 
2554  // Create histogram to store eta,phi limits (these need to be passed to StEStructHAdd).
2555  StEStructBinning* b = StEStructBinning::Instance();
2556  mHEtaPhi = new TH2D("EtaPhiRange","EtaPhiRange",1,b->etaMin(),b->etaMax(),1,b->phiMin(),b->phiMax());
2557 
2558  for(int i=0; i<8; i++){
2559 
2560  if (!mdontFillYtYt) {
2561  mHYtYt[i]=new TH2D*[numCutBins];
2562  mHNYtYt[i]=new TH2D*[numCutBins];
2563  }
2564 
2565  mHJtDEtaDPhi[i]=new TH2D*[numCutBins];
2566  if (!mdontFillMeanPt) {
2567  mHPrJtDEtaDPhi[i]=new TH2D*[numCutBins];
2568  mHPaJtDEtaDPhi[i]=new TH2D*[numCutBins];
2569  mHPbJtDEtaDPhi[i]=new TH2D*[numCutBins];
2570  }
2571 
2572  if (mdoFillEtaEta) {
2573  mHEtaEta[i]=new TH2D*[numCutBins];
2574  mHPhiPhi[i]=new TH2D*[numCutBins];
2575  mHNPhiPhi[i]=new TH2D*[numCutBins];
2576  if (mFillASSS) {
2577  mHEtaEtaSS[i]=new TH2D*[numCutBins];
2578  mHEtaEtaAS[i]=new TH2D*[numCutBins];
2579  }
2580  if (!mdontFillMeanPt) {
2581  mHPrEtaEta[i]=new TH2D*[numCutBins];
2582  mHPrPhiPhi[i]=new TH2D*[numCutBins];
2583  mHPaEtaEta[i]=new TH2D*[numCutBins];
2584  mHPaPhiPhi[i]=new TH2D*[numCutBins];
2585  mHPbEtaEta[i]=new TH2D*[numCutBins];
2586  mHPbPhiPhi[i]=new TH2D*[numCutBins];
2587  if (mFillASSS) {
2588  mHPrEtaEtaSS[i]=new TH2D*[numCutBins];
2589  mHPaEtaEtaSS[i]=new TH2D*[numCutBins];
2590  mHPbEtaEtaSS[i]=new TH2D*[numCutBins];
2591  mHPrEtaEtaAS[i]=new TH2D*[numCutBins];
2592  mHPaEtaEtaAS[i]=new TH2D*[numCutBins];
2593  mHPbEtaEtaAS[i]=new TH2D*[numCutBins];
2594  }
2595  }
2596  }
2597 
2598  if (mdoFillSumHistograms) {
2599  mHAtSYtDYt[i]=new TH2D*[numCutBins];
2600  mHAtNSYtDYt[i]=new TH2D*[numCutBins];
2601  mHJtSEtaDPhi[i]=new TH2D*[numCutBins];
2602  mHJtNSEtaDPhi[i]=new TH2D*[numCutBins];
2603  if (!mdontFillMeanPt) {
2604  mHPrJtSEtaDPhi[i]=new TH2D*[numCutBins];
2605  mHPaJtSEtaDPhi[i]=new TH2D*[numCutBins];
2606  mHPbJtSEtaDPhi[i]=new TH2D*[numCutBins];
2607  }
2608  }
2609 
2610  if (mFillQInv) {
2611  mHQinv[i]=new TH1D*[numCutBins];
2612  mHNQinv[i]=new TH1D*[numCutBins];
2613  }
2614 
2615  if(mdoPairDensityHistograms) {
2616  mHTPCAvgTSep[i]=new TH1D*[nden];
2617  mHTPCAvgZSep[i]=new TH1D*[nden];
2618  mHTPCEntTSep[i]=new TH1D*[nden];
2619  mHTPCEntZSep[i]=new TH1D*[nden];
2620  mHTPCMidTSep[i]=new TH1D*[nden];
2621  mHTPCMidZSep[i]=new TH1D*[nden];
2622  mHTPCExitTSep[i]=new TH1D*[nden];
2623  mHTPCExitZSep[i]=new TH1D*[nden];
2624  mHTPCMidTdptP[i]=new TH1D*[nden];
2625  mHTPCMidTdptN[i]=new TH1D*[nden];
2626  mHTPCMidZdptP[i]=new TH1D*[nden];
2627  mHTPCMidZdptN[i]=new TH1D*[nden];
2628  mHTPCQuality[i]=new TH1D*[nden];
2629  mHTPCAvgTZ[i]=new TH2D*[nden];
2630  mHTPCEntTZ[i]=new TH2D*[nden];
2631  mHTPCMidTZ[i]=new TH2D*[nden];
2632  mHTPCMidTZC[i]=new TH2D*[nden];
2633  mHTPCMidTZNC[i]=new TH2D*[nden];
2634  mHTPCExitTZ[i]=new TH2D*[nden];
2635  mHTPCEntQZ[i]=new TH2D*[nden];
2636  mHTPCMidQZ[i]=new TH2D*[nden];
2637  mHTPCEntQT[i]=new TH2D*[nden];
2638  mHTPCMidQT[i]=new TH2D*[nden];
2639  mHTPCEntQZT[i]=new TH2D*[nden];
2640  mHTPCMidQZT[i]=new TH2D*[nden];
2641  mHTPCEntTdpt[i]=new TH2D*[nden];
2642  mHTPCMidTdpt[i]=new TH2D*[nden];
2643  mHTPCExitTdpt[i]=new TH2D*[nden];
2644  }
2645 
2646  }
2647 }
2648 
2649 //--------------------------------------------------------------------------
2650 void StEStruct2ptCorrelations::deleteHistograms(){
2651 
2652  int numCutBins=StEStructCutBin::Instance()->getNumBins();
2653  int nden=StEStructCutBin::Instance()->getNumPairDensityBins();
2654  int numParentBins=StEStructCutBin::Instance()->getNumParentBins();
2655  if (mZBufferCutBinning) {
2656  numCutBins *= kNumBuffers;
2657  nden *= kNumBuffers;
2658  numParentBins *= kNumBuffers;
2659  }
2660 
2661  delete mHEtaPhi;
2662  for(int j=0;j<numParentBins;j++){
2663  delete mHMeanPtTot[j];
2664  delete mHMeanPtP[j];
2665  delete mHMeanPtM[j];
2666  delete mHMeanYtTot[j];
2667  delete mHMeanYtP[j];
2668  delete mHMeanYtM[j];
2669  delete mHEtaTot[j];
2670  delete mHEtaP[j];
2671  delete mHEtaM[j];
2672  }
2673  delete mHPtTot[0];
2674  delete mHPtP[0];
2675  delete mHPtM[0];
2676  delete mHPtTot[1];
2677  delete mHPtP[1];
2678  delete mHPtM[1];
2679  delete mHPtTot[2];
2680  delete mHPtP[2];
2681  delete mHPtM[2];
2682  delete mHPtTot[3];
2683  delete mHPtP[3];
2684  delete mHPtM[3];
2685  delete mHYtTot[0];
2686  delete mHYtP[0];
2687  delete mHYtM[0];
2688  delete mHYtTot[1];
2689  delete mHYtP[1];
2690  delete mHYtM[1];
2691  delete mHYtTot[2];
2692  delete mHYtP[2];
2693  delete mHYtM[2];
2694  delete mHYtTot[3];
2695  delete mHYtP[3];
2696  delete mHYtM[3];
2697  delete mHPhiAssocTot;
2698  delete mHPhiAssocP;
2699  delete mHPhiAssocM;
2700  delete mHPhiAssocPtTot;
2701  delete mHPhiAssocPtP;
2702  delete mHPhiAssocPtM;
2703  delete mHPtTrigTot;
2704  delete mHPtTrigP;
2705  delete mHPtTrigM;
2706  delete mHYtTrigTot;
2707  delete mHYtTrigP;
2708  delete mHYtTrigM;
2709  for(int i=0;i<8;i++){
2710  for(int j=0;j<numCutBins;j++){
2711  if (!mdontFillYtYt) {
2712  delete mHYtYt[i][j];
2713  delete mHNYtYt[i][j];
2714  }
2715 
2716  delete mHJtDEtaDPhi[i][j];
2717  if (!mdontFillMeanPt) {
2718  delete mHPrJtDEtaDPhi[i][j];
2719  delete mHPaJtDEtaDPhi[i][j];
2720  delete mHPbJtDEtaDPhi[i][j];
2721  }
2722 
2723  if (mdoFillEtaEta) {
2724  delete mHEtaEta[i][j];
2725  delete mHPhiPhi[i][j];
2726  delete mHNPhiPhi[i][j];
2727  if (mFillASSS) {
2728  delete mHEtaEtaSS[i][j];
2729  delete mHEtaEtaAS[i][j];
2730  }
2731  if (!mdontFillMeanPt) {
2732  delete mHPrEtaEta[i][j];
2733  delete mHPrPhiPhi[i][j];
2734  delete mHPaEtaEta[i][j];
2735  delete mHPaPhiPhi[i][j];
2736  delete mHPbEtaEta[i][j];
2737  delete mHPbPhiPhi[i][j];
2738  if (mFillASSS) {
2739  delete mHPrEtaEtaSS[i][j];
2740  delete mHPaEtaEtaSS[i][j];
2741  delete mHPbEtaEtaSS[i][j];
2742  delete mHPrEtaEtaAS[i][j];
2743  delete mHPaEtaEtaAS[i][j];
2744  delete mHPbEtaEtaAS[i][j];
2745  }
2746  }
2747  }
2748 
2749  if (mdoFillSumHistograms) {
2750  delete mHAtSYtDYt[i][j];
2751  delete mHAtNSYtDYt[i][j];
2752  delete mHJtSEtaDPhi[i][j];
2753  delete mHJtNSEtaDPhi[i][j];
2754  if (!mdontFillMeanPt) {
2755  delete mHPrJtSEtaDPhi[i][j];
2756  delete mHPaJtSEtaDPhi[i][j];
2757  delete mHPbJtSEtaDPhi[i][j];
2758  }
2759  }
2760 
2761  if (mFillQInv) {
2762  delete mHQinv[i][j];
2763  delete mHNQinv[i][j];
2764  }
2765  }
2766  if(mdoPairDensityHistograms) {
2767  for(int j=0;j<nden;j++){
2768  delete mHTPCAvgTSep[i][j];
2769  delete mHTPCAvgZSep[i][j];
2770  delete mHTPCEntTSep[i][j];
2771  delete mHTPCEntZSep[i][j];
2772  delete mHTPCMidTSep[i][j];
2773  delete mHTPCMidZSep[i][j];
2774  delete mHTPCExitTSep[i][j];
2775  delete mHTPCExitZSep[i][j];
2776  delete mHTPCMidTdptP[i][j];
2777  delete mHTPCMidTdptN[i][j];
2778  delete mHTPCMidZdptP[i][j];
2779  delete mHTPCMidZdptN[i][j];
2780  delete mHTPCQuality[i][j];
2781  delete mHTPCAvgTZ[i][j];
2782  delete mHTPCEntTZ[i][j];
2783  delete mHTPCMidTZ[i][j];
2784  delete mHTPCMidTZC[i][j];
2785  delete mHTPCMidTZNC[i][j];
2786  delete mHTPCExitTZ[i][j];
2787  delete mHTPCEntQZ[i][j];
2788  delete mHTPCMidQZ[i][j];
2789  delete mHTPCEntQT[i][j];
2790  delete mHTPCMidQT[i][j];
2791  delete mHTPCEntQZT[i][j];
2792  delete mHTPCMidQZT[i][j];
2793  delete mHTPCEntTdpt[i][j];
2794  delete mHTPCMidTdpt[i][j];
2795  delete mHTPCExitTdpt[i][j];
2796  }
2797 
2798  }
2799 
2800  if (!mdontFillYtYt) {
2801  delete [] mHYtYt[i];
2802  delete [] mHNYtYt[i];
2803  }
2804 
2805  delete [] mHJtDEtaDPhi[i];
2806  if (!mdontFillMeanPt) {
2807  delete [] mHPrJtDEtaDPhi[i];
2808  delete [] mHPaJtDEtaDPhi[i];
2809  delete [] mHPbJtDEtaDPhi[i];
2810  }
2811 
2812  if (mdoFillEtaEta) {
2813  delete [] mHEtaEta[i];
2814  delete [] mHPhiPhi[i];
2815  delete [] mHNPhiPhi[i];
2816  if (mFillASSS) {
2817  delete [] mHEtaEtaSS[i];
2818  delete [] mHEtaEtaAS[i];
2819  }
2820  if (!mdontFillMeanPt) {
2821  delete [] mHPrEtaEta[i];
2822  delete [] mHPrPhiPhi[i];
2823  delete [] mHPaEtaEta[i];
2824  delete [] mHPaPhiPhi[i];
2825  delete [] mHPbEtaEta[i];
2826  delete [] mHPbPhiPhi[i];
2827  if (mFillASSS) {
2828  delete [] mHPrEtaEtaSS[i];
2829  delete [] mHPaEtaEtaSS[i];
2830  delete [] mHPbEtaEtaSS[i];
2831  delete [] mHPrEtaEtaAS[i];
2832  delete [] mHPaEtaEtaAS[i];
2833  delete [] mHPbEtaEtaAS[i];
2834  }
2835  }
2836  }
2837 
2838  if (mdoFillSumHistograms) {
2839  delete [] mHAtSYtDYt[i];
2840  delete [] mHAtNSYtDYt[i];
2841  delete [] mHJtSEtaDPhi[i];
2842  delete [] mHJtNSEtaDPhi[i];
2843  if (!mdontFillMeanPt) {
2844  delete [] mHPrJtSEtaDPhi[i];
2845  delete [] mHPaJtSEtaDPhi[i];
2846  delete [] mHPbJtSEtaDPhi[i];
2847  }
2848  }
2849 
2850  if (mFillQInv) {
2851  delete [] mHQinv[i];
2852  delete [] mHNQinv[i];
2853  }
2854 
2855  if(mdoPairDensityHistograms) {
2856  delete [] mHTPCAvgTSep[i];
2857  delete [] mHTPCAvgZSep[i];
2858  delete [] mHTPCEntTSep[i];
2859  delete [] mHTPCEntZSep[i];
2860  delete [] mHTPCMidTSep[i];
2861  delete [] mHTPCMidZSep[i];
2862  delete [] mHTPCExitTSep[i];
2863  delete [] mHTPCExitZSep[i];
2864  delete [] mHTPCMidTdptP[i];
2865  delete [] mHTPCMidTdptN[i];
2866  delete [] mHTPCMidZdptP[i];
2867  delete [] mHTPCMidZdptN[i];
2868  delete [] mHTPCQuality[i];
2869  delete [] mHTPCAvgTZ[i];
2870  delete [] mHTPCEntTZ[i];
2871  delete [] mHTPCMidTZ[i];
2872  delete [] mHTPCMidTZC[i];
2873  delete [] mHTPCMidTZNC[i];
2874  delete [] mHTPCExitTZ[i];
2875  delete [] mHTPCEntQZ[i];
2876  delete [] mHTPCMidQZ[i];
2877  delete [] mHTPCEntQT[i];
2878  delete [] mHTPCMidQT[i];
2879  delete [] mHTPCEntQZT[i];
2880  delete [] mHTPCMidQZT[i];
2881  delete [] mHTPCEntTdpt[i];
2882  delete [] mHTPCMidTdpt[i];
2883  delete [] mHTPCExitTdpt[i];
2884  }
2885 
2886  }
2887 
2888 }
2889 
2890 //-----------------------------------------------------------------
2891 void StEStruct2ptCorrelations::createHist2D(TH2D*** h, const char* name, int iknd, int icut, int zcut, int ncut, int nx, float xmin, float xmax, int ny, float ymin, float ymax ){
2892 
2893  TString hname(bName[iknd]);
2894  hname+=name;
2895  hname += "_cutBin_";
2896  hname += icut;
2897  hname += "_zBuf_";
2898  hname += zcut;
2899  TString htitle(bTitle[iknd]);
2900  htitle+=name;
2901  h[iknd][ncut]=new TH2D(hname.Data(),htitle.Data(),nx,xmin,xmax,ny,ymin,ymax);
2902 
2903 }
2904 //-----------------------------------------------------------------
2905 void StEStruct2ptCorrelations::createHist1D(TH1D*** h, const char* name, int iknd, int icut, int zcut, int ncut, int nx, float xmin, float xmax){
2906 
2907  TString hname(bName[iknd]);
2908  hname+=name;
2909  hname += "_cutBin_";
2910  hname += icut;
2911  hname += "_zBuf_";
2912  hname += zcut;
2913  TString htitle(bTitle[iknd]);
2914  htitle+=name;
2915  h[iknd][ncut]=new TH1D(hname.Data(),htitle.Data(),nx,xmin,xmax);
2916 
2917 }
2918 //-----------------------------------------------------------------
2919 void StEStruct2ptCorrelations::createHist1D(TH1F*** h, const char* name, int iknd, int icut, int zcut, int ncut, int nx, float xmin, float xmax){
2920 
2921  TString hname(bName[iknd]);
2922  hname+=name;
2923  hname += "_cutBin_";
2924  hname += icut;
2925  hname += "_zBuf_";
2926  hname += zcut;
2927  TString htitle(bTitle[iknd]);
2928  htitle+=name;
2929  h[iknd][ncut]=new TH1F(hname.Data(),htitle.Data(),nx,xmin,xmax);
2930 
2931 }
2932 
2933 
2934 /***********************************************************************
2935  *
2936  * $Log: StEStruct2ptCorrelations.cxx,v $
2937  * Revision 1.31 2013/02/08 19:32:43 prindle
2938  * Added "Triggered" histograms in StEStruct2ptCorrelations.
2939  * Protected against using tracks cuts in StEStruct2ptCorrelations when reading EStruct format events.
2940  * Added comment in EventMaker/StEStructTrack.cxx pointing out need to set BField correctly
2941  * when reading EStruct format events. (This should be read from file somehow, but...)
2942  *
2943  * Revision 1.30 2012/11/16 21:22:27 prindle
2944  * 2ptCorrelations: SS, AS histograms. Get eta limits from cuts. Fit PtAll histogram. Add histograms to keep track of eta, phi limits. A few more histograms
2945  * Binning: Add quality cut.
2946  * CutBin: modify mode9
2947  * PairCuts: modify goodDeltaZ for case of one track leaving via endcap.
2948  *
2949  * Revision 1.29 2011/08/02 20:34:02 prindle
2950  * More detailed histograms for event mixing.
2951  * Buffer: increased mixed events to 4 (from 2)
2952  * CutBin: added mode 9 for exploration of p_t space, fixed place in mode 5 where
2953  * histogram was written before checking it existed.
2954  * OneBuffer: added ZDC coincidence rate to event sorting space.
2955  *
2956  * Revision 1.28 2010/09/02 21:24:07 prindle
2957  * 2ptCorrelations: Fill histograms for event mixing information
2958  * Option for common mixing buffer
2959  * Switch to selectively fill QInv histograms (which take a long time)
2960  * CutBin: Moved PID code to Track class from Pair class. Needed to update this code.
2961  * PairCuts: Moved PID code from here to Track class.
2962  * Removed unnecessary creation of StThreeVector which seem to take a long time
2963  * Add ToF momentum cuts, modify dEdx momentum cuts. (Now allow dEdx to be
2964  * be identified up to 15GeV/c, ToF up to 10GeV/c.)
2965  *
2966  * Revision 1.27 2010/03/02 21:45:27 prindle
2967  * Had a problem with pair cuts when one track exited via endplate
2968  * Calculate maxDEta properly
2969  * Warning if you try turning histograms for pair cuts on
2970  *
2971  * Revision 1.26 2009/11/09 21:32:41 prindle
2972  * Fix warnings about casting char * to a const char * by redeclaring as const char *.
2973  *
2974  * Revision 1.25 2009/05/08 00:09:54 prindle
2975  * In 2ptCorrelations we added switches to select blocks of histograms to fill.
2976  * (See constructor in StEStruct2ptCorrelations.cxx)
2977  * Use a brute force method for checking crossing cuts. I had too many corner
2978  * cases with my clever check.
2979  * In Binning, change Yt limit and add methods for accessing number of histogram bins
2980  * to use (used in Support)
2981  *
2982  * Revision 1.24 2008/12/02 23:45:04 prindle
2983  * Changed switchYt to switchXX (etc.) to better reflect function.
2984  * Change minYt to 1.0 in Binning so YtYt histogram doesn't have empty lower bin (pt = 0.164 for yt = 1.0)
2985  * In CutBin: remove initPtBin
2986  * add mode 8
2987  * add notSymmetrized (used in Support)
2988  * Added LUT (Look Up Table) for pair cuts. Experimental for now.
2989  * Modified cutMerging2 (to look at track separation at a few radii)
2990  * and cutCrossing2 so it doesn't accidentally reject almost back to back tracks.
2991  *
2992  * Revision 1.23 2008/03/19 22:06:00 prindle
2993  * Added doInvariantMass flag.
2994  * Added some plots in pairDensityHistograms.
2995  * SetZOffset used to only be done when doPairDensity was true.
2996  * Moved creating/copying pairDensity histograms to same place as other histograms.
2997  * Added cutBinHistMode
2998  * mode3 neck was defined as yt1<2.2 && yt2<2.2 (and not soft)
2999  * now is 1.8<yt1<2.2 && 1.8<yt2<2.2
3000  * Added gooddzdxy, Merging2 and Crossing2 to pair cuts.
3001  *
3002  * Revision 1.22 2007/11/26 19:55:22 prindle
3003  * In 2ptCorrelations: Support for keeping all z-bins of selected centralities
3004  * Change way \hat{p_t} is calculated for parent distributions in pid case.
3005  * Binning Added parent binning (for \hat{p_t}
3006  * CutBin: Mode 5 extensively modified.
3007  * Added invariant mass cuts (probably a bad idea in general.)
3008  *
3009  * Revision 1.21 2007/05/27 22:45:00 msd
3010  * Added new cut bin modes 2 (soft/hard SS/AS), 6 (z-vertex binning), and 7 (modes 2*6).
3011  * Fixed bug in merging cut.
3012  * Added a few histograms to 2pt corr.
3013  *
3014  * Revision 1.20 2007/01/26 17:17:04 msd
3015  * Implemented new binning scheme: dEta stored in array with bin centered at zero, dPhi array has bins centered at zero and pi. Final DEtaDPhi has 25x25 bins with dPhi bin width of pi/12 so all major angles are centered in bins.
3016  *
3017  * Revision 1.19 2006/10/02 22:20:51 prindle
3018  * Store only quadrant of eta_Delta - phi_Delta array/histogram.
3019  * Store half of eta_Sigma - phi_Delta array/histogram.
3020  * This required modifications in Binning.
3021  * I had a bug in the pair loop (which left +- not fully symmetrized)
3022  * and had to make changes in cut bins for mode 5 (and 3 I think)
3023  * when I fixed this.
3024  * Also change crossing cut to use only two parameters, the sign of
3025  * the magnetic field being taken from the MuDst.
3026  *
3027  * Revision 1.18 2006/05/03 18:14:43 msd
3028  * Fixed pair density again, removed references to iSwitch and symmetrize for simplicity
3029  *
3030  * Revision 1.17 2006/05/03 17:52:11 msd
3031  * Fixed pair density plots broken by recent symmetry changes
3032  *
3033  * Revision 1.16 2006/04/27 22:40:35 porter
3034  * 3 changes: 1) added npair hists for errors needed with eta_delta weighting
3035  * 2) commented out a few histograms to trim memory usage
3036  * 3) changed all hists to double precision (reflected in createHists member functions)
3037  *
3038  * Revision 1.15 2006/04/13 23:02:35 prindle
3039  * Added comment about not deleting output root file (which gets reported
3040  * as a memory leak but this is just before the program quits anayway.)
3041  *
3042  * Revision 1.14 2006/04/12 19:09:07 porter
3043  * added logic should z-vertex cut be ommitted (i.e. analysis of event generators)
3044  *
3045  * Revision 1.13 2006/04/10 23:42:32 porter
3046  * Added sameSide() & awaySide() methods to PairCut (so only defined in 1 place)
3047  * and added the eta_delta weighting as a binned correctin defined by the eta-limits in
3048  * the StEStructBinning object
3049  *
3050  * Revision 1.12 2006/04/06 01:01:11 prindle
3051  *
3052  * New mode in CutBin, 5, to do pid correlations. There is still an issue
3053  * of how to set the pt ranges allowed for the different particle types.
3054  * For data we probably want to restrict p to below 1GeV for pi and K, but
3055  * for Hijing and Pythia we can have perfect pid. Currently cuts are type
3056  * into the code (so you have to re-compile to change them.)
3057  *
3058  * In the Correlations code I split -+ from +- and am keeping track of
3059  * pt for each cut bin. These required changes in the Support code.
3060  *
3061  * Revision 1.11 2006/04/04 22:10:07 porter
3062  * a handful of changes (specific to correlations)
3063  * - added StEStructQAHists so that if NOT input frm Maker, each analysis has its own
3064  * - used ability to get any max,min val from the cut class - or z-vertex binning
3065  * - put z-vertex binning into 1 place
3066  * - switched back 1st line of pair cut method to keep pair if good, not to reject if bad.
3067  * - Pair cut object is now pointer in correlations
3068  * - some diagnostic printouts available from macro
3069  * - Duncan's delta-phi binning change
3070  *
3071  * Revision 1.10 2006/02/22 22:05:10 prindle
3072  * Removed all references to multRef (?)
3073  * Added cut mode 5 for particle identified correlations.
3074  * Other cut modes should be same as before
3075  *
3076  * Revision 1.8 2005/10/10 16:22:51 msd
3077  * fixing bug in buffer initialization
3078  *
3079  * Revision 1.7 2005/09/14 17:14:17 msd
3080  * Large update, added new pair-cut system, added pair density plots for new analysis mode (4), added event mixing cuts (rewrote buffer for this)
3081  *
3082  * Revision 1.6 2005/09/07 20:21:13 prindle
3083  *
3084  * 2ptCorrelations: Rearranged array/histogram initialization/destruction.
3085  * Now histograms are only allocated at end of job,
3086  * just before they are filled then written.
3087  *
3088  * Revision 1.5 2005/03/03 01:30:43 porter
3089  * updated StEStruct2ptCorrelations to include pt-correlations and removed
3090  * old version of pt-correlations from chunhuih (StEStruct2ptPtNbar)
3091  *
3092  * Revision 1.4 2004/09/16 23:44:05 chunhuih
3093  *
3094  * call init() from doEvent() instead of the constructor. This is needed
3095  * because we want to use the polymorphic feature of this virtual function.
3096  * The constructor would always call the local version of a virtual function.
3097  *
3098  * Revision 1.3 2004/07/01 00:34:52 porter
3099  * correct accounting for possible pairs in stats files
3100  *
3101  * Revision 1.2 2004/06/25 03:11:48 porter
3102  * New cut-binning implementation and modified pair-cuts for chunhui to review
3103  *
3104  * Revision 1.1 2003/10/15 18:20:46 porter
3105  * initial check in of Estruct Analysis maker codes.
3106  *
3107  *
3108  *********************************************************************/
3109 
3110 
Generated by   doxygen 1.8.5