StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WriteOutPIDTableMacro.C
1 #include <float.h>
2 
3 #include <fstream.h>
4 
5 #include "/afs/rhic.bnl.gov/star/replicas/DEV/StRoot/StEventUtilities/BetheBlochFunction.hh"
6 
7 #include "/star/data05/scratch/aihong/pidamp_dAu/StRoot/StPidAmpMaker/StPidProbabilityConst.hh"
8 
9 
10 void WriteOutPIDTableMacro( char* myOutputName){
11 
13 
14  char* mInputAmpFileName[mMultiplicityBins][mNDcaBins][mNChargeBins];
15 
16  char* mInputSigmaFileName[mMultiplicityBins][mNDcaBins][mNChargeBins];
17  char* mInputCalibFileName[mMultiplicityBins][mNDcaBins][mNChargeBins];
18 
19  char* mOutputFileName;
20 
21  double* mSigmaOfSingleTrail;
22  mSigmaOfSingleTrail = new double[mNEtaBins];
23 
24  TF1* EBandCenter;
25  TF1* PiBandCenter;
26  TF1* KBandCenter;
27  TF1* PBandCenter;
28 
29  TF1* BandCenterPtr;
30 
31 
32  double **BBPrePar;
33  double **BBTurnOver;
34  double **BBOffSetPar;
35  double **BBScalePar;
36  double **BBSaturate;
37 
38 
39  BBPrePar = new double*[mNEtaBins];
40  BBTurnOver = new double*[mNEtaBins];
41  BBOffSetPar = new double*[mNEtaBins];
42  BBScalePar = new double*[mNEtaBins];
43  BBSaturate = new double*[mNEtaBins];
44 
45 
46 
47  for (int ii=0; ii<mNEtaBins; ii++) {
48  BBPrePar[ii] = new double[mNNHitsBins];
49  BBTurnOver[ii] = new double[mNNHitsBins];
50  BBOffSetPar[ii] = new double[mNNHitsBins];
51  BBScalePar[ii] = new double[mNNHitsBins];
52  BBSaturate[ii] = new double[mNNHitsBins];
53  }
54 
55 
56  //allocate functions...
57  EBandCenter
58  =new TF1("EBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
59  PiBandCenter
60  =new TF1("PiBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
61  KBandCenter
62  =new TF1("KBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
63  PBandCenter
64  =new TF1("PBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
65 
66 
67 
68  Double_t electronPars[7]
69  ={ 1.072, 0.3199, 1.66032e-07, 1, 0.511e-3, 2.71172e-07, 0.0005 };
70  Double_t pionPars[7]
71  ={ 1.072, 0.3199, 1.66032e-07, 1, 0.13957, 2.71172e-07, 0.0005 };
72  Double_t kaonPars[7]
73  ={ 1.072, 0.3199, 1.66032e-07, 1, 0.49368, 2.71172e-07, 0.0005 };
74  Double_t antiprotonPars[7]
75  ={ 1.072, 0.3199, 1.66032e-07, 1, 0.93827, 2.71172e-07, 0.0005 };
76 
77 
78  EBandCenter->SetParameters(&electronPars[0]);
79  PiBandCenter->SetParameters(&pionPars[0]);
80  KBandCenter->SetParameters(&kaonPars[0]);
81  PBandCenter->SetParameters(&antiprotonPars[0]);
82 
83 
85 
86 
87 
88  TString OutPutName(myOutputName);
89 
90  char* PidTag;
91 
92  if (OutPutName.Contains("PiPIDTable"))
93  {PidTag="Pi"; BandCenterPtr=PiBandCenter;}
94  else if (OutPutName.Contains("EPIDTable"))
95  {PidTag="E"; BandCenterPtr=EBandCenter;}
96  else if (OutPutName.Contains("KPIDTable"))
97  {PidTag="K"; BandCenterPtr=KBandCenter;}
98  else if (OutPutName.Contains("PPIDTable"))
99  {PidTag="P"; BandCenterPtr=PBandCenter;}
100 
101  //set file names
102  for (int i=0; i<mMultiplicityBins; i++)
103  for (int j=0; j<mNDcaBins; j++)
104  for (int k=0; k<mNChargeBins; k++) {
105 
106 
107 
108  //-------------amp file
109  mInputAmpFileName[i][j][k]= new char[80];
110  sprintf(mInputAmpFileName[i][j][k],"./PidHistoAmp_%d%d%d.root",i,j,k);
111 
112  //------------sigma file
113  mInputSigmaFileName[i][j][k]= new char[80];
114  sprintf(mInputSigmaFileName[i][j][k],"PidSigmaOfSingleTrail_%d%d%d_basedOn_%d00.txt",i,j,k,i);
115 
116  //------------calib file
117  mInputCalibFileName[i][j][k]= new char[80];
118  sprintf(mInputCalibFileName[i][j][k], "./PhaseSpaceCalib%d%d%dButItisbasedOn_%d01.txt",i,j,k,i);
119 
120  //-----------out put file
121  mOutputFileName=myOutputName;
122  }
123 
124 
126  TFile outFile(mOutputFileName,"RECREATE");
127 
128  outFile.cd();
129 
130 
131  int thisMultBins=mMultiplicityBins;
132  int thisDcaBins=mNDcaBins;
133  int thisChargeBins=mNChargeBins;
134 
135  int thisDedxBins=200; //fake
136  int thisPBins=200;
137  int thisEtaBins=mNEtaBins;
138  int thisNHitsBins=mNNHitsBins;
139 
141 
142  double thisDedxStart=0.;
143  double thisDedxEnd=0.53e-5;
144  double thisPStart=1e-12;
145  double thisPEnd=2.0;
146  double thisEtaStart=mEtaStart;
147  double thisEtaEnd=mEtaEnd;
148  double thisNHitsStart=mNNHitsStart;
149  double thisNHitsEnd=mNNHitsEnd;
150 
152 //p, eta, nhits can be equavily dividable.
153 
154  TVectorD* EqualyDividableRangeStartSet = new TVectorD(20);
155 
156  (*EqualyDividableRangeStartSet)(0)=thisDedxStart;
157  (*EqualyDividableRangeStartSet)(1)=thisPStart;
158  (*EqualyDividableRangeStartSet)(2)=thisEtaStart;
159  (*EqualyDividableRangeStartSet)(3)=thisNHitsStart;
160 
161  EqualyDividableRangeStartSet->Write("EqualyDividableRangeStartSet",TObject::kOverwrite | TObject::kSingleKey);
162 //-------------
163  TVectorD* EqualyDividableRangeEndSet = new TVectorD(20);
164 
165  (*EqualyDividableRangeEndSet)(0)=thisDedxEnd;
166  (*EqualyDividableRangeEndSet)(1)=thisPEnd;
167  (*EqualyDividableRangeEndSet)(2)=thisEtaEnd;
168  (*EqualyDividableRangeEndSet)(3)=thisNHitsEnd;
169 
170  EqualyDividableRangeEndSet->Write("EqualyDividableRangeEndSet",TObject::kOverwrite | TObject::kSingleKey);
171 
172 //-------------
173 
174  TVectorD* EqualyDividableRangeNBinsSet = new TVectorD(20);
175 
176  (*EqualyDividableRangeNBinsSet)(0)=thisDedxBins;
177  (*EqualyDividableRangeNBinsSet)(1)=thisPBins;
178  (*EqualyDividableRangeNBinsSet)(2)=thisEtaBins;
179  (*EqualyDividableRangeNBinsSet)(3)=thisNHitsBins;
180 
181  EqualyDividableRangeNBinsSet->Write("EqualyDividableRangeNBinsSet",TObject::kOverwrite | TObject::kSingleKey);
182 
184 // multi, dca, charge , not easy to be equavily dividable
185 
186  TVectorD* NoEqualyDividableRangeNBinsSet = new TVectorD(20);
187 
188  (*NoEqualyDividableRangeNBinsSet)(0)=thisMultBins;
189  (*NoEqualyDividableRangeNBinsSet)(1)=thisDcaBins;
190  (*NoEqualyDividableRangeNBinsSet)(2)=thisChargeBins;
191 
192  NoEqualyDividableRangeNBinsSet->Write("NoEqualyDividableRangeNBinsSet",TObject::kOverwrite | TObject::kSingleKey);
193 
194 //----------
195  TVectorD* MultiBinEdgeSet = new TVectorD(20);
196  // for dAu
197  (*MultiBinEdgeSet)(0) = 1.; //
198  (*MultiBinEdgeSet)(1) =.4;
199  (*MultiBinEdgeSet)(2) =.2; //20% central
200 
201  MultiBinEdgeSet->Write("MultiBinEdgeSet",TObject::kOverwrite | TObject::kSingleKey);
202 //----------
203  TVectorD* DcaBinEdgeSet = new TVectorD(20);
204 
205  (*DcaBinEdgeSet)(0) = 0.0000; //
206  (*DcaBinEdgeSet)(1) = 3.0; // 3cm
207 
208  DcaBinEdgeSet->Write("DcaBinEdgeSet",TObject::kOverwrite | TObject::kSingleKey);
209 
210 //----------charge bin edge is obviouse, charge<0, bin0. charge>0 bin1. so forget it.
211 
212 
213  outFile.Write();
214 
215 
216  int objAryIdx=0;
217 
218  TVectorD* theAmpTable;
219  TVectorD* theCenterTable;
220  TVectorD* theSigmaTable;
221 
222 
223  TVectorD* theBBPreParTable;
224  TVectorD* theBBTurnOverTable;
225  TVectorD* theBBScaleTable;
226  TVectorD* theBBOffSetTable;
227  TVectorD* theBBSaturateTable;
228 
229 
230 
233 
234 
235  //-------fill and write eAmp
236  objAryIdx=0;
237  theAmpTable = new TVectorD(thisMultBins*thisDcaBins*thisChargeBins*thisPBins*thisEtaBins*thisNHitsBins);
238  theCenterTable = new TVectorD(thisMultBins*thisDcaBins*thisChargeBins*thisPBins*thisEtaBins*thisNHitsBins);
239  theSigmaTable = new TVectorD(thisMultBins*thisDcaBins*thisChargeBins*thisPBins*thisEtaBins*thisNHitsBins);
240 
241 
242 
243  theBBPreParTable = new TVectorD(thisEtaBins*thisNHitsBins);
244  theBBTurnOverTable = new TVectorD(thisEtaBins*thisNHitsBins);
245  theBBScaleTable = new TVectorD(thisEtaBins*thisNHitsBins);
246  theBBOffSetTable = new TVectorD(thisEtaBins*thisNHitsBins);
247  theBBSaturateTable = new TVectorD(thisEtaBins*thisNHitsBins);
248 
249 
250  for (int multIdx=0; multIdx<thisMultBins; multIdx++)
251  for (int dcaIdx=0; dcaIdx<thisDcaBins; dcaIdx++)
252  for (int chargeIdx=0; chargeIdx<thisChargeBins; chargeIdx++){
253 
254  TFile AmpFile(mInputAmpFileName[multIdx][dcaIdx][chargeIdx],"READ");
255 
256  ifstream sigmaFile(mInputSigmaFileName[multIdx][dcaIdx][chargeIdx]);
257  ifstream calibFile(mInputCalibFileName[multIdx][dcaIdx][chargeIdx]);
258 
259 
260  cout<<mInputSigmaFileName[multIdx][dcaIdx][chargeIdx]<<endl;
261 
262  int ie,in;
263 //read in calib parameters.
264  for (int h=0; h<thisEtaBins*thisNHitsBins;h++){
265  calibFile>>ie; calibFile>>in;
266 
267 
268  calibFile>>BBPrePar[ie][in];
269  calibFile>>BBTurnOver[ie][in];
270  calibFile>>BBOffSetPar[ie][in];
271  calibFile>>BBScalePar[ie][in];
272  calibFile>>BBSaturate[ie][in];
273 
274 
275  (*theBBPreParTable)(h)=BBPrePar[ie][in];
276  (*theBBTurnOverTable)(h)=BBTurnOver[ie][in];
277  (*theBBOffSetTable)(h)=BBOffSetPar[ie][in];
278  (*theBBScaleTable)(h)=BBScalePar[ie][in];
279  (*theBBSaturateTable)(h)=BBSaturate[ie][in];
280 
281  }
282 
283  calibFile.close();
284 
285  for (int h=0; h<thisEtaBins; h++){
286  sigmaFile>>ie; sigmaFile>>mSigmaOfSingleTrail[h];
287  }
288 
289  sigmaFile.close();
290 
291 
292  for (int pIdx=0; pIdx<thisPBins; pIdx++)
293  for (int etaIdx=0; etaIdx<thisEtaBins; etaIdx++)
294  for (int nhitsIdx=0; nhitsIdx<thisNHitsBins; nhitsIdx++){
295 
296 
297  double pPosition= (pIdx+0.5)*(thisPEnd-thisPStart)/float(thisPBins);
298  double etaPosition= (etaIdx+0.5)*(thisEtaEnd-thisEtaStart)/float(thisEtaBins);
299  double nhitsPosition=double((nhitsIdx+0.5)*(float(thisNHitsEnd-thisNHitsStart)/float(thisNHitsBins)));
300 
301  BandCenterPtr->SetParameter(0,BBPrePar[etaIdx][nhitsIdx]);
302  BandCenterPtr->SetParameter(1,BBTurnOver[etaIdx][nhitsIdx]);
303  BandCenterPtr->SetParameter(2,BBOffSetPar[etaIdx][nhitsIdx]);
304  BandCenterPtr->SetParameter(5,BBScalePar[etaIdx][nhitsIdx]);
305  BandCenterPtr->SetParameter(6,BBSaturate[etaIdx][nhitsIdx]);
306 
307  double theAmp;
308  double theCenter;
309  double theSigma;
310 
311  char *theAmpName = new char[80];
312  sprintf(theAmpName,"%sAmp%d%d",PidTag,etaIdx,nhitsIdx);
313 
314  TH1F* theAmpHist=(TH1F *)AmpFile.Get(theAmpName);
315 
316  if (theAmpName) delete theAmpName;
317 
318  if (theAmpHist) {
319 
320 
321 
322  if (OutPutName.Contains("PiPIDTable"))
323  {PidTag="Pi"; BandCenterPtr=PiBandCenter;}
324  if (OutPutName.Contains("EPIDTable"))
325  {PidTag="E"; BandCenterPtr=EBandCenter;}
326  if (OutPutName.Contains("KPIDTable"))
327  {PidTag="K"; BandCenterPtr=KBandCenter;}
328  if (OutPutName.Contains("PPIDTable"))
329  {PidTag="P"; BandCenterPtr=PBandCenter;}
330 
331 
332 
333  if (OutPutName.Contains("PiPIDTable")){//3 fcn for Pi
334 
335  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
336  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
337 
338  TF1* PiFcnCenter = theAmpHist->GetFunction("PiFcnCenter");
339  TF1* PiFcnLeft = theAmpHist->GetFunction("PiFcnLeft");
340  TF1* PiFcnRight = theAmpHist->GetFunction("PiFcnRight");
341 
342  if (PiFcnRight){
343  PiFcnRight->SetRange(PiFcnRight->GetXmin(),mPEnd);
344  if ( pPosition<=PiFcnRight->GetXmax()
345  && pPosition >PiFcnRight->GetXmin()&& PiFcnRight->GetParameter(1)<0 )
346  theAmp = PiFcnRight->Eval(pPosition,0.,0.);
347  }
348 
349  if (PiFcnCenter){
350  if ( pPosition<=PiFcnCenter->GetXmax()
351  && pPosition >PiFcnCenter->GetXmin() )
352  theAmp = PiFcnCenter->Eval(pPosition,0.,0.);
353  }
354 
355  if (PiFcnLeft){
356  if ( pPosition <= PiFcnLeft->GetXmax()
357  && pPosition > PiFcnLeft->GetXmin() )
358  theAmp = PiFcnLeft->Eval(pPosition,0.,0.);
359  }
360 
361 
362  } else if (OutPutName.Contains("EPIDTable")) { //1 fcn for E
363 
364  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
365  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
366 
367  TF1* EFcnLeft = theAmpHist->GetFunction("EFcnLeft");
368  if (EFcnLeft){
369  if ( pPosition<=EFcnLeft->GetXmax()
370  && pPosition >EFcnLeft->GetXmin() && EFcnLeft->GetParameter(1)<0 )
371  theAmp = EFcnLeft->Eval(pPosition,0.,0.);
372  }
373 
374  TF1* EFcnRight = theAmpHist->GetFunction("EFcnRight");
375  if (EFcnRight){
376  // EFcnRight->SetRange(EFcnRight->GetXmin(), mPEnd);
377  EFcnRight->SetRange(0.15, mPEnd); //extrapolate it to lower pt
378  if ( pPosition<=EFcnRight->GetXmax()
379  && pPosition >EFcnRight->GetXmin() && EFcnRight->GetParameter(1)<0 )
380  theAmp = EFcnRight->Eval(pPosition,0.,0.);
381  }
382 
383  if (nhitsPosition<=5) //ndedx<5, E fcn fitting not good. use hist.
384  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
385  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
386 
387  } else if (OutPutName.Contains("KPIDTable")) {//3 fcn for K
388 
389  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
390  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
391 
392 
393  TF1* KFcnCenter = theAmpHist->GetFunction("KFcnCenter");
394  TF1* KFcnLeft = theAmpHist->GetFunction("KFcnLeft");
395  TF1* KFcnRight = theAmpHist->GetFunction("KFcnRight");
396 
397  if (KFcnRight){
398  if ( pPosition<=KFcnRight->GetXmax()
399  && pPosition >KFcnRight->GetXmin() )
400  theAmp = KFcnRight->Eval(pPosition,0.,0.);
401  }
402 
403  if (KFcnCenter){
404  if ( pPosition<=KFcnCenter->GetXmax()
405  && pPosition >KFcnCenter->GetXmin() )
406  theAmp = KFcnCenter->Eval(pPosition,0.,0.);
407  }
408 
409  if (KFcnLeft){
410  if ( pPosition<=KFcnLeft->GetXmax()
411  && pPosition >KFcnLeft->GetXmin() )
412  theAmp = KFcnLeft->Eval(pPosition,0.,0.);
413  }
414 
415  if (nhitsPosition<=5) //ndedx<5, K fcn fitting not good. use hist.
416  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
417  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
418 
419  } else if (OutPutName.Contains("PPIDTable")) { // 1 fcn for P
420 
421  theAmp = (theAmpHist->GetBinContent(pIdx+1)>0) ?
422  (theAmpHist->GetBinContent(pIdx+1)) : 0.;
423 
424  TF1* PFcnRight = theAmpHist->GetFunction("PFcnRight");
425  TF1* PFcnLeft = theAmpHist->GetFunction("PFcnLeft");
426 
427  if (PFcnRight){
428  if ( pPosition<=PFcnRight->GetXmax()
429  && pPosition >PFcnRight->GetXmin() )
430  theAmp = PFcnRight->Eval(pPosition,0.,0.);
431  }
432 
433  if (PFcnLeft){
434  if ( pPosition<=PFcnLeft->GetXmax()
435  && pPosition >PFcnLeft->GetXmin() )
436  theAmp = PFcnLeft->Eval(pPosition,0.,0.);
437  }
438 
439  }
440 
441 
442  } else theAmp=0;
443 
444 
445  theCenter=BandCenterPtr->Eval(pPosition);
446  theSigma=
447  mSigmaOfSingleTrail[etaIdx]*theCenter/TMath::Sqrt(nhitsPosition);
448 
449 
450 
451  theAmp= (theAmp<FLT_MAX && theAmp>0.) ? theAmp : 0.;
452  theCenter = (theCenter<FLT_MAX && theCenter>0.) ? theCenter : 0.;
453  theSigma = (theSigma<FLT_MAX && theSigma>0. ) ? theSigma : 0.;
454 
455  (*theAmpTable)(objAryIdx)=(theAmp);
456  (*theCenterTable)(objAryIdx)=(theCenter);
457  (*theSigmaTable)(objAryIdx)=(theSigma);
458 
459  // cout<<" mSigmaOfSingleTrail["<<etaIdx<<"] "<<mSigmaOfSingleTrail[etaIdx]<<" theAmp "<<theAmp<<" theCenter "<<theCenter<<" theSigma "<<theSigma<<endl;
460 
461  objAryIdx++;
462 
463  // cout<<objAryIdx<<endl;
464 
465  }
466 
467  AmpFile.Close();
468 
469  }
470 
471  outFile.cd();
472 
473 
474  char *ampNm = new char[80]; // name of amp. table
475  sprintf(ampNm,"%sAmp",PidTag);
476  char *ctrNm = new char[80]; // name of center table
477  sprintf(ctrNm,"%sCenter",PidTag);
478  char *sigmaNm = new char[80];// name of sigma table
479  sprintf(sigmaNm,"%sSigma",PidTag);
480 
481  theAmpTable->Write(ampNm,TObject::kOverwrite | TObject::kSingleKey);
482  theCenterTable->Write(ctrNm,TObject::kOverwrite | TObject::kSingleKey);
483  theSigmaTable->Write(sigmaNm,TObject::kOverwrite | TObject::kSingleKey);
484 
485  theBBPreParTable->Write("BBPrePar",TObject::kOverwrite | TObject::kSingleKey);
486  theBBTurnOverTable->Write("BBTurnOver",TObject::kOverwrite | TObject::kSingleKey);
487  theBBOffSetTable->Write("BBOffSet",TObject::kOverwrite | TObject::kSingleKey);
488  theBBScaleTable->Write("BBScale",TObject::kOverwrite | TObject::kSingleKey);
489  theBBSaturateTable->Write("BBSaturate",TObject::kOverwrite | TObject::kSingleKey);
490 
491 
492  outFile.Write();
493 
494  if (theAmpTable) delete theAmpTable;
495  if (theCenterTable) delete theCenterTable;
496  if (theSigmaTable) delete theSigmaTable;
497 
498  if (theBBScaleTable) delete theBBScaleTable;
499  if (theBBOffSetTable) delete theBBOffSetTable;
500 
501  if (ampNm) delete ampNm;
502  if (ctrNm) delete ctrNm;
503  if (sigmaNm) delete sigmaNm;
504 
505  outFile.Close();
506 }
507 
508 
509 
510 
511 
512 
513 
514 
515 
516 
517 
518 
519 
520 
521 
522 
523 
524 
525 
526 
527 
528 
529 
530 
531 
532 
533 
534 
535 
536 
537 //submit.pl event /star/data12/reco/ProductionMinBias/ReversedFullField/P01gl/2001/269/ /star/data10/GC/aihong/P01glPIDFlowPicoDST/
Generated by   doxygen 1.8.5