StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StGlauberTree.cxx
1 //====================================================================================================
2 // $Id: StGlauberTree.cxx,v 1.2 2012/04/25 04:45:26 hmasui Exp $
3 // $Log: StGlauberTree.cxx,v $
4 // Revision 1.2 2012/04/25 04:45:26 hmasui
5 // Expand branches for eccentricity in the main tree, and deformation parameters in header
6 //
7 //====================================================================================================
8 
9 #include <assert.h>
10 
11 #include "TError.h"
12 #include "TFile.h"
13 #include "TMath.h"
14 #include "TTree.h"
15 #include "StMessMgr.h"
16 #include "StGlauberTree.h"
17 
18 ClassImp(StGlauberTree)
19 
20 //____________________________________________________________________________________________________
21 // Default constructor
22 StGlauberTree::StGlauberTree(const UInt_t mode)
23  : mMode(mode)
24 {
25  switch ( mMode ){
26  case 0: LOG_INFO << "StGlauberTree read mode" << endm; break ;
27  case 1: LOG_INFO << "StGlauberTree write mode" << endm; break ;
28  default:
29  Error("StGlauberTree", "Mode should 0(read) or 1(write). abort");
30  assert(0);
31  }
32 
33  mFile = 0 ;
34  mTree = 0 ;
35 
36  Clear() ;
37 
38  // Initialize header stuffs only once
39  sprintf(mNameNucleusA, "%s", "") ;
40  sprintf(mNameNucleusB, "%s", "") ;
41  mMassNumberA = 0 ;
42  mMassNumberB = 0 ;
43  mRadiusA = -9999. ;
44  mRadiusA = -9999. ;
45  mSkinDepthA = -9999. ;
46  mSkinDepthB = -9999. ;
47  mBeta2A = -9999. ;
48  mBeta4A = -9999. ;
49  mBeta2B = -9999. ;
50  mBeta4B = -9999. ;
51  mSigmaNN = -9999. ;
52  mSqrtSNN = -9999. ;
53  mRepulsionD = -9999. ;
54  mTotalXsec = -9999. ;
55  mTotalXsecError = 0.0 ;
56  mSmearHardCore = 0 ;
57  mSmearGaussian = 0 ;
58  mCollisionHardCore = 0 ;
59  mCollisionGaussian = 0 ;
60  mBMax = -9999. ;
61  mNeventsAccept = 0 ;
62  mNeventsThrow = 0 ;
63  mNpp = -9999. ;
64  mK = -9999. ;
65  mX = -9999. ;
66  mEfficiency = -9999. ;
67  mIsConstEfficiency = 0 ;
68  mVersion = 0 ;
69 }
70 
71 //____________________________________________________________________________________________________
72 // Default destructor
73 StGlauberTree::~StGlauberTree()
74 {
75 }
76 
77 //____________________________________________________________________________________________________
78 Int_t StGlauberTree::Clear()
79 {
81  mB = -9999. ;
82  mNpart = 0 ;
83  mNcoll = 0 ;
84  mMultiplicity = 0 ;
85 
86  for(UInt_t i=0;i<2;i++){
87  mTheta[i] = 0.0 ;
88  mPhi[i] = 0.0 ;
89  }
90 
91  for(UInt_t i=0;i<4;i++){
92  mSumX[i] = 0.0 ;
93  mSumY[i] = 0.0 ;
94  mSumX2[i] = 0.0 ;
95  mSumY2[i] = 0.0 ;
96  mSumXY[i] = 0.0 ;
97  mEccRP2[i] = -9999. ;
98  mEccPP2[i] = -9999. ;
99  mEccPP3[i] = -9999. ;
100  mEccPP4[i] = -9999. ;
101  mPP2[i] = -9999. ;
102  mPP3[i] = -9999. ;
103  mPP4[i] = -9999. ;
104  }
105 
106  return 1 ;
107 }
108 
109 //____________________________________________________________________________________________________
110 Int_t StGlauberTree::InitBranch()
111 {
112  // event-wise tree
113  mTree->SetMakeClass(1);
114  mTree->SetBranchAddress("b", &mB, &b_b);
115  mTree->SetBranchAddress("npart", &mNpart, &b_npart);
116  mTree->SetBranchAddress("ncoll", &mNcoll, &b_ncoll);
117  mTree->SetBranchAddress("mult", &mMultiplicity, &b_mult);
118  mTree->SetBranchAddress("theta", mTheta, &b_theta);
119  mTree->SetBranchAddress("phi", mPhi, &b_phi);
120  mTree->SetBranchAddress("sumx", mSumX, &b_sumx);
121  mTree->SetBranchAddress("sumy", mSumY, &b_sumy);
122  mTree->SetBranchAddress("sumx2", mSumX2, &b_sumx2);
123  mTree->SetBranchAddress("sumy2", mSumY2, &b_sumy2);
124  mTree->SetBranchAddress("sumxy", mSumXY, &b_sumxy);
125  mTree->SetBranchAddress("eccrp2", mEccRP2, &b_eccrp2);
126  mTree->SetBranchAddress("eccpp2", mEccPP2, &b_eccpp2);
127  mTree->SetBranchAddress("eccpp3", mEccPP3, &b_eccpp3);
128  mTree->SetBranchAddress("eccpp4", mEccPP4, &b_eccpp4);
129  mTree->SetBranchAddress("pp2", mPP2, &b_pp2);
130  mTree->SetBranchAddress("pp3", mPP3, &b_pp3);
131  mTree->SetBranchAddress("pp4", mPP4, &b_pp4);
132 
133  // header
134  mHeader->SetMakeClass(1);
135  mHeader->SetBranchAddress("nameA", mNameNucleusA, &b_nameA);
136  mHeader->SetBranchAddress("nameB", mNameNucleusB, &b_nameB);
137  mHeader->SetBranchAddress("massNumberA", &mMassNumberA, &b_massNumberA);
138  mHeader->SetBranchAddress("massNumberB", &mMassNumberB, &b_massNumberB);
139  mHeader->SetBranchAddress("radiusA", &mRadiusA, &b_radiusA);
140  mHeader->SetBranchAddress("radiusB", &mRadiusB, &b_radiusB);
141  mHeader->SetBranchAddress("skinDepthA", &mSkinDepthA, &b_skinDepthA);
142  mHeader->SetBranchAddress("skinDepthB", &mSkinDepthB, &b_skinDepthB);
143  mHeader->SetBranchAddress("beta2A", &mBeta2A, &b_beta2A);
144  mHeader->SetBranchAddress("beta4A", &mBeta4A, &b_beta4A);
145  mHeader->SetBranchAddress("beta2B", &mBeta2B, &b_beta2B);
146  mHeader->SetBranchAddress("beta4B", &mBeta4B, &b_beta4B);
147  mHeader->SetBranchAddress("sigmaNN", &mSigmaNN, &b_sigmaNN);
148  mHeader->SetBranchAddress("sqrtSNN", &mSqrtSNN, &b_sqrtSNN);
149  mHeader->SetBranchAddress("repulsionD", &mRepulsionD, &b_repulsionD);
150  mHeader->SetBranchAddress("totalXsec", &mTotalXsec, &b_totalXsec);
151  mHeader->SetBranchAddress("totalXsecError", &mTotalXsecError, &b_totalXsecError);
152  mHeader->SetBranchAddress("smearHardCore", &mSmearHardCore, &b_smearHardCore);
153  mHeader->SetBranchAddress("smearGaussian", &mSmearGaussian, &b_smearGaussian);
154  mHeader->SetBranchAddress("collisionHardCore", &mCollisionHardCore, &b_collisionHardCore);
155  mHeader->SetBranchAddress("collisionGaussian", &mCollisionGaussian, &b_collisionGaussian);
156  mHeader->SetBranchAddress("maxB", &mBMax, &b_maxB);
157  mHeader->SetBranchAddress("neventsAccept", &mNeventsAccept, &b_neventsAccept);
158  mHeader->SetBranchAddress("neventsThrow", &mNeventsThrow, &b_neventsThrow);
159  mHeader->SetBranchAddress("npp", &mNpp, &b_npp);
160  mHeader->SetBranchAddress("k", &mK, &b_k);
161  mHeader->SetBranchAddress("x", &mX, &b_x);
162  mHeader->SetBranchAddress("efficiency", &mEfficiency, &b_efficiency);
163  mHeader->SetBranchAddress("isConstEfficiency", &mIsConstEfficiency, &b_isConstEfficiency);
164  mHeader->SetBranchAddress("version", &mVersion, &b_version);
165 
166  return 1;
167 }
168 
169 
170 //____________________________________________________________________________________________________
171 Int_t StGlauberTree::Open(const TString filename)
172 {
173  if(mMode==0){
174  // Read mode
175  mFile = TFile::Open(filename);
176  if(!mFile || !mFile->IsOpen()){
177  Error("StGlauberTree::Open", "can't open %s", filename.Data());
178  assert(0);
179  }
180  LOG_INFO << "StGlauberTree::Open Open input ROOT file: " << mFile->GetName() << endm;
181 
182  // Get tree
183  mTree = (TTree*) mFile->Get("tree");
184  if(!mTree){
185  Error("StGlauberTree::Open", "No tree found in %s", filename.Data());
186  assert(0);
187  }
188 
189  mHeader = (TTree*) mFile->Get("header");
190  if(!mHeader){
191  Error("StGlauberTree::Open", "No header found in %s", filename.Data());
192  assert(0);
193  }
194 
195  InitBranch() ;
196  }
197  else{
198  // Make sure mFile/mTree/mHeader are NULL
199  if(mFile) delete mFile ;
200  if(mTree) delete mTree ;
201  if(mHeader) delete mHeader ;
202 
203  // Write mode
204  mFile = TFile::Open(filename, "recreate", "", 9); // max compression
205  if(!mFile || !mFile->IsOpen()){
206  Error("StGlauberTree::Open", "can't open %s", filename.Data());
207  assert(0);
208  }
209  LOG_INFO << "StGlauberTree::Open Open output ROOT file: " << mFile->GetName() << endm;
210 
211  // Initialization of tree
212  LOG_INFO << "StGlauberTree::Open Init tree" << endm;
213  mTree = new TTree("tree", "Event-wise information for Fast MC glauber tree");
214  mTree->Branch("b", &mB, "b/D") ; // Impact parameter
215  mTree->Branch("npart", &mNpart, "npart/i") ; // Number of participants
216  mTree->Branch("ncoll", &mNcoll, "ncoll/i") ; // Number of collicipants
217  mTree->Branch("mult", &mMultiplicity, "mult/i") ; // Multiplicity
218  mTree->Branch("theta", mTheta, "theta[2]/D") ; // Polar angle rotation for deformed nuclei
219  mTree->Branch("phi", mPhi, "phi[2]/D") ; // Azimuthal angle rotation for deformed nuclei
220  mTree->Branch("sumx", mSumX, "sumx[4]/D") ; // <x>
221  mTree->Branch("sumy", mSumY, "sumy[4]/D") ; // <y>
222  mTree->Branch("sumx2", mSumX2, "sumx2[4]/D") ; // <x^2>
223  mTree->Branch("sumy2", mSumY2, "sumy2[4]/D") ; // <y^2>
224  mTree->Branch("sumxy", mSumXY, "sumxy[4]/D") ; // <xy>
225  mTree->Branch("eccrp2", mEccRP2, "eccrp2[4]/D") ; // 2nd order Reaction plane eccentricity
226  mTree->Branch("eccpp2", mEccPP2, "eccpp2[4]/D") ; // 2nd order Participant plane eccentricity
227  mTree->Branch("eccpp3", mEccPP3, "eccpp3[4]/D") ; // 3rd order Participant plane eccentricity
228  mTree->Branch("eccpp4", mEccPP4, "eccpp4[4]/D") ; // 4th order Participant plane eccentricity
229  mTree->Branch("pp2", mPP2, "pp2[4]/D") ; // 2nd order participant plane
230  mTree->Branch("pp3", mPP3, "pp3[4]/D") ; // 3rd order participant plane
231  mTree->Branch("pp4", mPP4, "pp4[4]/D") ; // 4th order participant plane
232 
233  mHeader = new TTree("header", "Constants used in the Fast MC glauber");
234  mHeader->Branch("nameA", mNameNucleusA, "nameA[2]/C");
235  mHeader->Branch("nameB", mNameNucleusB, "nameB[2]/C");
236  mHeader->Branch("massNumberA", &mMassNumberA, "massNumberA/i");
237  mHeader->Branch("massNumberB", &mMassNumberB, "massNumberB/i");
238  mHeader->Branch("radiusA", &mRadiusA, "radiusA/F");
239  mHeader->Branch("radiusB", &mRadiusB, "radiusB/F");
240  mHeader->Branch("skinDepthA", &mSkinDepthA, "skinDepthA/F");
241  mHeader->Branch("skinDepthB", &mSkinDepthB, "skinDepthB/F");
242  mHeader->Branch("beta2A", &mBeta2A, "beta2A/F");
243  mHeader->Branch("beta4A", &mBeta4A, "beta4A/F");
244  mHeader->Branch("beta2B", &mBeta2B, "beta2B/F");
245  mHeader->Branch("beta4B", &mBeta4B, "beta4B/F");
246  mHeader->Branch("sigmaNN", &mSigmaNN, "sigmaNN/F");
247  mHeader->Branch("sqrtSNN", &mSqrtSNN, "sqrtSNN/F");
248  mHeader->Branch("repulsionD", &mRepulsionD, "repulsionD/F");
249  mHeader->Branch("totalXsec", &mTotalXsec, "totalXsec/F");
250  mHeader->Branch("totalXsecError", &mTotalXsecError, "totalXsecError/F");
251  mHeader->Branch("smearHardCore", &mSmearHardCore, "smearHardCore/i");
252  mHeader->Branch("smearGaussian", &mSmearGaussian, "smearGaussian/i");
253  mHeader->Branch("collisionHardCore", &mCollisionHardCore, "collisionHardCore/i");
254  mHeader->Branch("collisionGaussian", &mCollisionGaussian, "collisionGaussian/i");
255  mHeader->Branch("maxB", &mBMax, "maxB/F");
256  mHeader->Branch("neventsAccept", &mNeventsAccept, "neventsAccept/i");
257  mHeader->Branch("neventsThrow", &mNeventsThrow, "neventsThrow/i");
258  mHeader->Branch("npp", &mNpp, "npp/F");
259  mHeader->Branch("k", &mK, "k/F");
260  mHeader->Branch("x", &mX, "x/F");
261  mHeader->Branch("efficiency", &mEfficiency, "efficiency/F");
262  mHeader->Branch("isConstEfficiency", &mIsConstEfficiency, "isConstEfficiency/i");
263  mHeader->Branch("version", &mVersion, "version/i");
264  }
265 
266  return 1;
267 }
268 
269 //____________________________________________________________________________________________________
270 void StGlauberTree::Sort()
271 {
273  if( mMode == 1 ) mFile->GetList()->Sort() ;
274 }
275 
276 //____________________________________________________________________________________________________
277 Int_t StGlauberTree::Fill()
278 {
279  return mTree->Fill();
280 }
281 
282 //____________________________________________________________________________________________________
283 Int_t StGlauberTree::FillHeader()
284 {
285  return mHeader->Fill();
286 }
287 
288 //____________________________________________________________________________________________________
289 Int_t StGlauberTree::Close()
290 {
291  LOG_INFO << "StGlauberTree::Close close ROOT file : " << mFile->GetName() << endm;
292  if(mMode==1) mFile->Write();
293  mFile->Close();
294 
295  return 1;
296 }
297 
298 //____________________________________________________________________________________________________
299 Int_t StGlauberTree::GetEntries() const
300 {
301  return mTree->GetEntries() ;
302 }
303 
304 //____________________________________________________________________________________________________
305 Int_t StGlauberTree::GetEntry(const Int_t ievent)
306 {
307  return mTree->GetEntry(ievent);
308 }
309 
310 //____________________________________________________________________________________________________
311 Double_t StGlauberTree::GetSigmaA2(const Double_t a2, const Double_t a) const
312 {
313  const Double_t val = a2 - a*a ;
314  if( val < 0.0 ){
315  Error("GetSigmaA2", "{a^2}-{a}^2 < 0 (= %1.3f)", val);
316  return -9999. ;
317  }
318 
319  return a2 - a*a ;
320 }
321 
322 //____________________________________________________________________________________________________
323 Double_t StGlauberTree::GetSigmaXY(const Double_t xy, const Double_t x, const Double_t y) const
324 {
325  return xy - x*y ;
326 }
327 
328 //____________________________________________________________________________________________________
329 Double_t StGlauberTree::GetSRP(const UInt_t id) const
330 {
332 
333  if( id >= 4 ){
334  Error("StGlauberTree::GetSRP", "Unknown id, id=%3d. abort", id);
335  assert(0);
336  }
337 
338  const Double_t sigmax2 = GetSigmaA2(mSumX2[id], mSumX[id]);
339  const Double_t sigmay2 = GetSigmaA2(mSumY2[id], mSumY[id]);
340  const Double_t area = TMath::Pi() * TMath::Sqrt(sigmax2*sigmay2) ;
341 
342  return area ;
343 }
344 
345 //____________________________________________________________________________________________________
346 Double_t StGlauberTree::GetSPP(const UInt_t id) const
347 {
349 
350  if( id >= 4 ){
351  Error("StGlauberTree::GetSPP", "Unknown id, id=%3d. abort", id);
352  assert(0);
353  }
354 
355  const Double_t sigmax2 = GetSigmaA2(mSumX2[id], mSumX[id]);
356  const Double_t sigmay2 = GetSigmaA2(mSumY2[id], mSumY[id]);
357  const Double_t sigmaxy = GetSigmaXY(mSumXY[id], mSumX[id], mSumY[id]);
358  const Double_t term = sigmax2*sigmay2 - sigmaxy*sigmaxy ;
359  if( term < 0 ) return -9999. ;
360 
361  return TMath::Pi() * TMath::Sqrt(term) ;
362 }
363 
364 //____________________________________________________________________________________________________
365 void StGlauberTree::SetNameNucleusA(const Char_t* val)
366 {
367  sprintf(mNameNucleusA, "%s", val);
368 }
369 
370 //____________________________________________________________________________________________________
371 void StGlauberTree::SetNameNucleusB(const Char_t* val)
372 {
373  sprintf(mNameNucleusB, "%s", val);
374 }
375 
StGlauberTree::Clear
Int_t Clear()
Default destructor.
Definition: StGlauberTree.cxx:78
StGlauberTree::GetSPP
Double_t GetSPP(const UInt_t id) const
Definition: StGlauberTree.cxx:346
Error
Definition: FJcore.h:367
StGlauberTree::GetEntries
Int_t GetEntries() const
Close ROOT file.
Definition: StGlauberTree.cxx:299
StGlauberTree::Sort
void Sort()
Open ROOT file, initialize tree.
Definition: StGlauberTree.cxx:270
StGlauberTree::FillHeader
Int_t FillHeader()
Fill event-wise tree.
Definition: StGlauberTree.cxx:283
StGlauberTree::GetEntry
Int_t GetEntry(const Int_t ievent)
Get entries.
Definition: StGlauberTree.cxx:305
StGlauberTree::Close
Int_t Close()
Fill header tree.
Definition: StGlauberTree.cxx:289
StGlauberTree::~StGlauberTree
virtual ~StGlauberTree()
Default constructor.
Definition: StGlauberTree.cxx:73
StGlauberTree::GetSRP
Double_t GetSRP(const UInt_t id) const
Definition: StGlauberTree.cxx:329
StGlauberTree::Open
Int_t Open(const TString filename)
Clear data members.
Definition: StGlauberTree.cxx:171
StGlauberTree::Fill
Int_t Fill()
Sort outputs in ROOT file.
Definition: StGlauberTree.cxx:277
Generated by   doxygen 1.8.5