StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtAVEfficiencyMaker.cxx
1 #include "StFgtAVEfficiencyMaker.h"
3 #include "StRoot/StEvent/StFgtCollection.h"
4 #include "StRoot/StEvent/StFgtHitCollection.h"
5 #include "StRoot/StEvent/StFgtHit.h"
6 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
7 #include "StRoot/StEvent/StEvent.h"
8 #include "StRoot/StEvent/StEventInfo.h"
9 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
10 #include <TH2D.h>
11 #include <TROOT.h>
12 #include <TStyle.h>
13 #include <TCanvas.h>
14 #include <utility>
15 #include <TArc.h>
16 #include <TLine.h>
17 #include <set>
18 
19 
20 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
21 
22 #include "StRoot/StEvent/StEvent.h"
23 #include "StRoot/StEvent/StFgtCollection.h"
24 
25 //Double_t StFgtAVEfficiencyMaker::getRPhiRatio()
26 //{
27 
28 //};
29 Double_t getRPhiRatio(StSPtrVecFgtHitConstIterator hitIterBegin, StSPtrVecFgtHitConstIterator hitIterEnd)
30 {
31  Short_t quad, disc, strip;
32  Char_t layer;
33  Int_t numR=0;
34  Int_t numPhi=0;
35 
36  StSPtrVecFgtHitConstIterator hitIter=hitIterBegin;
37  for(;hitIter!=hitIterEnd;hitIter++)
38  {
39  StFgtGeom::decodeGeoId((*hitIter)->getCentralStripGeoId(),disc, quad, layer, strip);
40  if(layer=='R')
41  numR++;
42  else
43  numPhi++;
44  }
45 
46  if(numR+numPhi>0)
47  return (numR-numPhi)/((Double_t)(numR+numPhi));
48  else
49  return -1;
50 };
51 
55 Int_t StFgtAVEfficiencyMaker::Make()
56 {
57  Int_t ierr = kStOk;
58 
59  StEvent* eventPtr = 0;
60  eventPtr = (StEvent*)GetInputDS("StEvent");
61 
62  if( !eventPtr ) {
63  LOG_ERROR << "Error getting pointer to StEvent from '" << ClassName() << "'" << endm;
64  ierr = kStErr;
65  };
66 
67  mFgtCollectionPtr = 0;
68 
69  if( eventPtr ) {
70  mFgtCollectionPtr=eventPtr->fgtCollection();
71  };
72 
73  if( !mFgtCollectionPtr) {
74  LOG_ERROR << "Error getting pointer to StFgtCollection from '" << ClassName() << "'" << endm;
75  ierr = kStErr;
76  };
77 
78  // cout <<" in eff make " <<endl;
79  Float_t x;
80  Float_t y;
81  Int_t prvGeoId=-1;
82 
83  StFgtHitCollection* clusterColD1=mFgtCollectionPtr->getHitCollection(0);
84  StFgtHitCollection* clusterColD6=mFgtCollectionPtr->getHitCollection(5);
85 
86  //look at the r phi ratio for each disk
87  for(int i=0;i<6;i++)
88  {
89  StFgtHitCollection* tmpClusterCol=mFgtCollectionPtr->getHitCollection(i);
90  Double_t ratio=getRPhiRatio(tmpClusterCol->getHitVec().begin(),tmpClusterCol->getHitVec().end());
91  rPhiRatioPlots[i]->Fill(ratio);
92  // cout << "ratio for disk: " << i << " is " << ratio <<" disk has: " << tmpClusterCol->getHitVec().size() << "hits" <<endl;
93  }
94 
95  const StSPtrVecFgtHit &hitVecD1=clusterColD1->getHitVec();
96  const StSPtrVecFgtHit &hitVecD6=clusterColD6->getHitVec();
97 
98 
99  set<Int_t> usedPoints;//saves the points that have been used for tracks (and shouldn't be reused)
100  Double_t D1Pos=StFgtGeom::getDiscZ(0);
101  Double_t D6Pos=StFgtGeom::getDiscZ(5);
102  Double_t zArm=D6Pos-D1Pos;
103  StSPtrVecFgtHitConstIterator hitIterD1,hitIterD6, hitIterD1R, hitIterD6R, hitIter, hitIter2;
104  // cout <<" there are " << hitVecD1.size() << " hits in d1 "<<endl;
105  for(hitIterD1=hitVecD1.begin();hitIterD1 != hitVecD1.end();hitIterD1++)
106  {
107  Short_t quad, disc, strip;
108  Char_t layer;
109  StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
110  // if(layer=='P')
111  // cout <<"found phi d1 " <<endl;
112  // else
113  // cout <<"found r d1 " <<endl;
114  }
115  // cout <<" there are " << hitVecD6.size() << " hits in d6 "<<endl;
116  for(hitIterD1=hitVecD6.begin();hitIterD1 != hitVecD6.end();hitIterD1++)
117  {
118  Short_t quad, disc, strip;
119  Char_t layer;
120  StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
121  // if(layer=='P')
122  // cout <<"found phi d6 " <<endl;
123  // else
124  // cout <<"found r d6 " <<endl;
125  }
126 
128  for(int iSeed1=0;iSeed1<5;iSeed1++)
129  {
130  for(int iSeed2=iSeed1+1;iSeed2<6;iSeed2++)
131  {
132  StFgtHitCollection* clusterColSeed1=mFgtCollectionPtr->getHitCollection(iSeed1);
133  StFgtHitCollection* clusterColSeed2=mFgtCollectionPtr->getHitCollection(iSeed2);
134  const StSPtrVecFgtHit &hitVecSeed1=clusterColSeed1->getHitVec();
135  const StSPtrVecFgtHit &hitVecSeed2=clusterColSeed2->getHitVec();
136  D1Pos=StFgtGeom::getDiscZ(iSeed1);
137  D6Pos=StFgtGeom::getDiscZ(iSeed2);
138  zArm=D6Pos-D1Pos;
139 
140  // cout <<"looking at " << hitVecD1.size() << " hits in D1 and " << hitVecD6.size() << " in D6 " <<endl;
141  // for(hitIterD1=hitVecD1.begin();hitIterD1 != hitVecD1.end();hitIterD1++)
142  for(hitIterD1=hitVecSeed1.begin();hitIterD1 != hitVecSeed1.end();hitIterD1++)
143  {
144  Short_t quad, disc, strip;
145  Char_t layer;
146  StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
148  if(quad>2)
149  continue;
150 
151  //do 1D 'fit' with r strips and the (x,y) thing
152  Int_t geoIdD1=(*hitIterD1)->getCentralStripGeoId();
153  StFgtGeom::decodeGeoId(geoIdD1,disc, quad, layer, strip);
154  if(layer!='P')
155  continue;
156 
157  Float_t phiD1=(*hitIterD1)->getPositionPhi();
158  for(hitIterD1R=hitVecSeed1.begin();hitIterD1R != hitVecSeed1.end();hitIterD1R++)
159  {
160  Int_t geoIdR=(*hitIterD1R)->getCentralStripGeoId();
161  StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
162  if(layer!='R')
163  continue;
164 
165  Float_t rD1=(*hitIterD1R)->getPositionR();
166  Float_t xD1=rD1*cos(phiD1);
167  Float_t yD1=rD1*sin(phiD1);
168  for(hitIterD6=hitVecSeed2.begin();hitIterD6 != hitVecSeed2.end();hitIterD6++)
169  {
170  Int_t geoIdD6=(*hitIterD6)->getCentralStripGeoId();
171  StFgtGeom::decodeGeoId(geoIdD6,disc, quad, layer, strip);
172  if(layer!='P')
173  continue;
174  Float_t phiD6=(*hitIterD6)->getPositionPhi();
175  for(hitIterD6R=hitVecSeed2.begin();hitIterD6R != hitVecSeed2.end();hitIterD6R++)
176  {
177  Int_t geoIdR=(*hitIterD6R)->getCentralStripGeoId();
178  StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
179  if(layer!='R')
180  continue;
181 
182  Float_t rD6=(*hitIterD6R)->getPositionR();
183  Double_t xD6=rD6*cos(phiD6);
184  Double_t yD6=rD6*sin(phiD6);
185 
187  vector< pair<Int_t,Double_t> > v_x;
188  vector< pair<Int_t,Double_t> > v_y;
189  vector< pair<Int_t,Double_t> > v_r;
190 
191  vector< pair<Int_t,Double_t> > v_xFail;
192  vector< pair<Int_t,Double_t> > v_yFail;
193  vector< pair<Int_t,Double_t> > v_rFail;
194 
195  vector< pair< Int_t, Double_t> > v_rCharge;
196  vector< pair< Int_t, Double_t> > v_phiCharge;
197 
198  vector<Int_t> v_clusterSizeR;
199  vector<Int_t> v_clusterSizePhi;
200 
201  vector<Int_t> v_geoIDsR;
202  vector<Int_t> v_geoIDsPhi;
203 
204  int iFound=0;
205  int iFoundR=0;
206 
207  //zarm is d6 position - D1
208  Double_t xIpExp=xD1+(xD6-xD1)*(-D1Pos)/zArm;
209  Double_t yIpExp=yD1+(yD6-yD1)*(-D1Pos)/zArm;
211  Double_t zIpExpX0=(D6Pos-(xD6/xD1)*D1Pos)*1/(1-xD6/xD1);
213  Double_t zIpExpY0=(D6Pos-yD6/yD1*D1Pos)*1/(1-yD6/yD1);
214  // cout <<"
215 
216  for(int iD=0;iD<6;iD++)
217  {
218  if(iD==iSeed1 || iD==iSeed2)
219  continue;
220 
221  Bool_t found=false;
222  Bool_t foundR=false;
223  //check for hit
224  Double_t diskZ=StFgtGeom::getDiscZ(iD);
225  //expected
226 
227  Double_t xPosExp=xD1+(xD6-xD1)*(diskZ-D1Pos)/zArm;
228  Double_t yPosExp=yD1+(yD6-yD1)*(diskZ-D1Pos)/zArm;
229  Double_t rPosExp=rD1+(rD6-rD1)*(diskZ-D1Pos)/zArm;
230 
231  //at z=0;
232 
233 
234  // cout <<"x1: " << xD1 << " y1: " << yD1 <<" x6: " << xD6 <<" y6: " << yD6 << " zpos: " << diskZ <<" arm: " << zArm<<endl;
235  // cout <<"expect hit at : " << xPosExp <<" / " << yPosExp <<" r: " << rPosExp <<endl;
236  StFgtHitCollection* clusterCol=mFgtCollectionPtr->getHitCollection(iD);
237  const StSPtrVecFgtHit &hitVec=clusterCol->getHitVec();
238  // cout <<" there are " << hitVec.size() << " hits in d" << iD <<endl;
239  for(hitIter=hitVec.begin();hitIter!=hitVec.end();hitIter++)
240  {
241  //do 1D 'fit' with r strips and the (x,y) thing
242  Int_t geoIdPhi=(*hitIter)->getCentralStripGeoId();
243  StFgtGeom::decodeGeoId(geoIdPhi,disc, quad, layer, strip);
244  if(usedPoints.find(geoIdPhi)!=usedPoints.end())
245  continue;
246  if(layer!='P')
247  continue;
248  Float_t phi=(*hitIter)->getPositionPhi();
249  Double_t phiCharge=(*hitIter)->charge();
250  Int_t clusterSizePhi=(*hitIter)->getStripWeightMap().size();
251  for(hitIter2=hitVec.begin();hitIter2!=hitVec.end();hitIter2++)
252  {
253  Int_t geoIdR=(*hitIter2)->getCentralStripGeoId();
254  StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
255  if(usedPoints.find(geoIdR)!=usedPoints.end())
256  continue;
257  if(layer!='R')
258  continue;
259  Float_t r=(*hitIter2)->getPositionR();
260  Double_t rCharge=(*hitIter2)->charge();
261  Int_t clusterSizeR=(*hitIter2)->getStripWeightMap().size();
262  // cout <<"checking with r:" << r <<endl;
263  if(fabs(r-rPosExp)<0.5)
264  {
265  foundR=true;
266  // cout <<"found r: " << r <<endl;
267  v_r.push_back(pair<Int_t,Double_t> (iD,r));
268  }
269 
270  x=r*cos(phi);
271  y=r*sin(phi);
272  // cout <<"checking with x: " << x << " y: " << y <<endl;
273  if(fabs(x-xPosExp) < 0.5 && fabs(y-yPosExp)<0.5) //found hit
274  {
275  found=true;
276  // cout <<"found! " <<" pushing back: iD: " << iD << "x: " << x << " y "<< y <<endl;
277  v_x.push_back(pair<Int_t,Double_t>(iD,x));
278  v_y.push_back(pair<Int_t,Double_t>(iD,y));
279  v_rCharge.push_back(pair<Int_t, Double_t>(iD,rCharge));
280  v_phiCharge.push_back(pair<Int_t, Double_t>(iD,phiCharge));
281  //to avoid double counting
282  v_geoIDsPhi.push_back(geoIdPhi);
283  v_geoIDsR.push_back(geoIdR);
284  // cout <<" we have rcharge: " << rCharge<<" phi: " << phiCharge <<endl;
285  v_clusterSizeR.push_back(clusterSizeR);
286  v_clusterSizePhi.push_back(clusterSizePhi);
287  }
288  }
289  }
290  //only one per disk
291  if(found)
292  iFound++;
293  else
294  {
295  // cout <<"failed to find, pushing back " << xPosExp << " y: " << yPosExp <<endl;
296  v_xFail.push_back(pair<Int_t,Double_t>(iD,xPosExp));
297  v_yFail.push_back(pair<Int_t,Double_t>(iD,yPosExp));
298  }
299  if(foundR)
300  iFoundR++;
301  else
302  {
303  // cout <<"failed to find r " << rPosExp<<endl;
304  v_rFail.push_back(pair<Int_t, Double_t>(iD,rPosExp));
305  }
306 
307  }
308 
309  // cout << " Ifound: " << iFound <<endl;
310  if(iFound>=2)
311  {
312  hIp->Fill(xIpExp,yIpExp);
313  hIpZAtX0->Fill(zIpExpX0);
314  hIpZAtY0->Fill(zIpExpY0);
315  // cout <<" we have zIPX0: " << zIpExpX0 <<" and y: " << zIpExpY0 <<endl;
316  }
317 
318 
319  // if(iFound>=2 && fabs(xIpExp)<20 && fabs(yIpExp)<20 && fabs(zIpExpX0)<40 && fabs(zIpExpY0)<40) //at least one hit plus seed and pointing roughly to the interaction region
320  if(iFound>=2 && fabs(zIpExpX0)<100 && fabs(zIpExpY0)<100) //at least one hit plus seed and pointing roughly to the interaction region
321  {
322  if(v_x.size()>iFound)
323  {
324  cout<<"more hits than disks hit!!!" <<endl;
325  cout <<" vx_size: " << v_x.size() <<" ifound: " << iFound <<endl;
326  }
327  else
328  {
329 
330  // cout <<"filled hip with " << xIpExp << " / " << yIpExp <<endl;
331  for(int i=0;i<v_x.size();i++)
332  {
333  usedPoints.insert(v_geoIDsPhi[i]);
334  usedPoints.insert(v_geoIDsR[i]);
335  Int_t disk=v_x[i].first;
336  Double_t x=v_x[i].second;
337  Double_t y=v_y[i].second;
338  Double_t rCharge=v_rCharge[i].second;
339  Double_t phiCharge=v_phiCharge[i].second;
340  // cout <<"filling disk: " << disk <<" with: " << x <<" / " <<y <<endl;
341  radioPlotsEff[disk]->Fill(x,y);
342  chargeCorr[disk]->Fill(rCharge,phiCharge);
343  h_clusterSizeR[disk]->Fill(v_clusterSizeR[i]);
344  h_clusterSizePhi[disk]->Fill(v_clusterSizePhi[i]);
345  h_clusterChargeR[disk]->Fill(rCharge);
346  h_clusterChargePhi[disk]->Fill(phiCharge);
347  }
348  for(int i=0;i<v_xFail.size();i++)
349  {
350  Int_t disk=v_xFail[i].first;
351  Double_t x=v_xFail[i].second;
352  Double_t y=v_yFail[i].second;
353  // cout <<"filling disk fail: " << disk <<" with: " << x <<" / " <<y <<endl;
354  radioPlotsNonEff[disk]->Fill(x,y);
355  }
356  hitCounter++;
357  }
358  }
359 
360  if(iFoundR>=1) //at least one hit plus seed
361  {
362  if(v_r.size()>iFound)
363  {
364  // cout<<"more r hits than disks hit!!!" <<endl;
365  }
366  else
367  {
368  for(int i=0;i<v_r.size();i++)
369  {
370  Int_t disk=v_r[i].first;
371  Double_t r=v_r[i].second;
372  // cout <<"filling r disk: " << disk <<" with: " << r <<endl;
373  rEff[disk]->Fill(r);
374  }
375  for(int i=0;i<v_rFail.size();i++)
376  {
377  Int_t disk=v_rFail[i].first;
378  Double_t r=v_rFail[i].second;
379  // cout <<"filling r disk fail: " << disk <<" with: " << r <<endl;
380  rNonEff[disk]->Fill(r);
381  }
382  hitCounterR++;
383  }
384  }
385 
386 
387  //start over
388  iFound=0;
389  iFoundR=0;
390  v_x.clear();
391  v_y.clear();
392  v_r.clear();
393  v_xFail.clear();
394  v_yFail.clear();
395  v_rFail.clear();
396 
397  }
398  }
399  }
400  }
401 
402  }
403  }
404 
405  return ierr;
406 
407 };
408 
409 StFgtAVEfficiencyMaker::StFgtAVEfficiencyMaker( const Char_t* name): StMaker( name ),runningEvtNr(0),hitCounter(0),hitCounterR(0)
410 {
411  cout <<"AVE constructor!!" <<endl;
412 
413 };
414 
415 StFgtAVEfficiencyMaker::~StFgtAVEfficiencyMaker()
416 {
417 
418  //delete histogram arrays
419 };
420 
421 
422 Int_t StFgtAVEfficiencyMaker::Finish(){
423  gStyle->SetPalette(1);
424  cout <<"cluster plotter finish funciton " <<endl;
425  Int_t ierr = kStOk;
426 
427  TCanvas* cRadio=new TCanvas("radioPlots","radioPlot",1000,1500);
428  TCanvas* cRadioHits=new TCanvas("radioPlotsHits","radioPlotHits",1000,1500);
429  TCanvas* cRadioNonHits=new TCanvas("radioPlotsNonHits","radioPlotNonHits",1000,1500);
430  cRadio->Divide(2,3); //6 discs
431  cRadioHits->Divide(2,3); //6 discs
432  cRadioNonHits->Divide(2,3); //6 discs
433  TCanvas* cRPRatio=new TCanvas("rPhiRatio","rPhiRatios",1000,1500);
434  cRPRatio->Divide(2,3); //6 discs
435 
436  TCanvas* cREff=new TCanvas("crEff","crEff",1000,1500);
437 
438  cREff->Divide(2,3); //6 discs
439 
440  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
441  {
442  // cRadio->cd(iD+1)->SetLogz();
443  cRadioHits->cd(iD+1);
444  radioPlotsEff[iD]->Draw("colz");
445  cRadioNonHits->cd(iD+1);
446  radioPlotsNonEff[iD]->Draw("colz");
447 
448  }
449  cRadioHits->SaveAs("radioPlotsHits.png");
450  cRadioNonHits->SaveAs("radioPlotsNonHits.png");
451 
452 
453  TCanvas* cClusterSizeR=new TCanvas("clusterSizeR","clusterSizeR",1000,1500);
454  cClusterSizeR->Divide(2,3);
455  TCanvas* cClusterSizePhi=new TCanvas("clusterSizePhi","clusterSizePhi",1000,1500);
456  cClusterSizePhi->Divide(2,3);
457  TCanvas* cChargeCorr=new TCanvas("chargeCorr","chargeCorr",1000,1500);
458  cChargeCorr->Divide(2,3);
459 
460  TCanvas* cClusterChargePhi=new TCanvas("clusterChargePhi","clusterChargePhi",1000,1500);
461  cClusterChargePhi->Divide(2,3);
462  TCanvas* cClusterChargeR=new TCanvas("clusterChargeR","clusterChargeR",1000,1500);
463  cClusterChargeR->Divide(2,3);
464 
465  TCanvas cIPProj;
466  hIp->Draw("colz");
467  cIPProj.SaveAs("ipProj.png");
468 
469  hIpZAtX0->Draw();
470  cIPProj.SaveAs("ipProjAtX0.png");
471 
472  hIpZAtY0->Draw();
473  cIPProj.SaveAs("ipProjAtY0.png");
474 
475 
476  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
477  {
478  cClusterSizeR->cd(iD+1);
479  h_clusterSizeR[iD]->Draw();
480  cClusterSizePhi->cd(iD+1);
481  h_clusterSizePhi[iD]->Draw();
482  cChargeCorr->cd(iD+1);
483  chargeCorr[iD]->Draw("colz");
484 
485  cClusterChargeR->cd(iD+1);
486  h_clusterChargeR[iD]->Draw();
487  cClusterChargePhi->cd(iD+1);
488  h_clusterChargePhi[iD]->Draw();
489 
490  }
491  cClusterSizeR->SaveAs("clusterSizeR.png");
492  cClusterSizePhi->SaveAs("clusterSizePhi.png");
493  cChargeCorr->SaveAs("chargeCorrelation.png");
494 
495  cClusterChargeR->SaveAs("clusterChargeR.png");
496  cClusterChargePhi->SaveAs("clusterChargePhi.png");
497 
498  cout <<"saving .." <<endl;
499 
500  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
501  {
502  cRadio->cd(iD+1);
503 
504  TH2D* tmpAllCounts=(TH2D*)radioPlotsEff[iD]->Clone("tmp");
505  radioPlotsEff[iD]->Add(radioPlotsNonEff[iD]);//all counts
506  // radioPlotsEff[iD]->Add(radioPlotsNonEff[iD],-1); //subtract non eff
507  for(int nx=0;nx<radioPlotsEff[iD]->GetNbinsX();nx++)
508  {
509  for(int ny=0;ny<radioPlotsEff[iD]->GetNbinsY();ny++)
510  {
511  Double_t denom=radioPlotsEff[iD]->GetBinContent(nx,ny);
512  if(denom>0)
513  radioPlotsEff[iD]->SetBinContent(nx,ny,tmpAllCounts->GetBinContent(nx,ny)/denom);
514  else
515  radioPlotsEff[iD]->SetBinContent(nx,ny,0.0);
516  }
517  }
518  // radioPlotsEff[iD]->Divide(tmpAllCounts);
519  radioPlotsEff[iD]->Draw("colz");
520  TArc *innerArc = new TArc(0,0,103.5+11.5,0.0,90.0);
521  TArc *outerArc = new TArc(0,0,394.0-11.5,0.0,90.0);
522  TLine *left = new TLine( 11.5, 103.5+11.5, 11.5, 394.0-11.5 );
523  TLine *right = new TLine( 103.5+11.5, 11.5, 394.0-11.5, 11.5 );
524  TLine *center = new TLine(
525  (103.5+11.5)*cos( 0.785398163 ),
526  (103.5+11.5)*sin( 0.785398163 ),
527  (394.0-11.5)*cos( 0.785398163 ),
528  (394.0-11.5)*sin( 0.785398163 )
529  );
530  TLine *weird = new TLine(
531  (394.0-11.5)*cos( 0.190240888 ),
532  (394.0-11.5)*sin( 0.190240888 ),
533  (394.0-11.5)*cos( 0.891863248 ),
534  (394.0-11.5)*sin( 0.891863248 )
535  );
536  innerArc->SetFillStyle(0);
537  outerArc->SetFillStyle(0);
538  innerArc->SetLineWidth(2);
539  outerArc->SetLineWidth(2);
540  left->SetLineWidth(2);
541  right->SetLineWidth(2);
542  center->SetLineWidth(2);
543  weird->SetLineWidth(2);
544  outerArc->Draw("only");
545  innerArc->Draw("only");
546  left->Draw();
547  right->Draw();
548  center->Draw();
549  weird->Draw();
550 
551 
552 
553 
554 
555  cRPRatio->cd(iD+1);
556  rPhiRatioPlots[iD]->Draw();
557  cREff->cd(iD+1);
558 
559  TH1D* tmpR=(TH1D*)rEff[iD]->Clone("tmpR");
560  rEff[iD]->Add(rNonEff[iD]);
561  for(int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
562  {
563  Double_t denom=rEff[iD]->GetBinContent(nx);
564  if(denom>0)
565  rEff[iD]->SetBinContent(nx,tmpR->GetBinContent(nx)/denom);
566  else
567  rEff[iD]->SetBinContent(nx,0.0);
568  }
569  rEff[iD]->Draw();
570  }
571  cRadio->SaveAs("radioPlotsEff.png");
572  cRadio->SaveAs("radioPlotsEff.pdf");
573 
574 
575  cREff->SaveAs("rEff.png");
576  cREff->SaveAs("rEff.pdf");
577 
578  cRPRatio->SaveAs("rpRatio.png");
579  cRPRatio->SaveAs("rpRatio.pdf");
580 
581  cout <<"returning after finish" <<endl;
582  return ierr;
583 };
584 
585 
590 Int_t StFgtAVEfficiencyMaker::Init(){
591  cout <<"AVE!!" <<endl;
592  myRootFile=new TFile("clusterEff.root","RECREATE");
593  // outTxtFile=new ofstream;
594  // outTxtFile->open("clusters.txt");
595 
596 
597  Int_t ierr = kStOk;
598 
599  char buffer[100];
600 
601  radioPlotsEff=new TH2D*[kFgtNumDiscs];
602  radioPlotsNonEff=new TH2D*[kFgtNumDiscs];
603  rPhiRatioPlots=new TH1D*[kFgtNumDiscs];
604  rEff=new TH1D*[kFgtNumDiscs];
605  rNonEff=new TH1D*[kFgtNumDiscs];
606 
607 
608  chargeCorr=new TH2D*[kFgtNumDiscs];
609  h_clusterSizeR=new TH1D*[kFgtNumDiscs];
610  h_clusterSizePhi=new TH1D*[kFgtNumDiscs];
611 
612  h_clusterChargeR=new TH1D*[kFgtNumDiscs];
613  h_clusterChargePhi=new TH1D*[kFgtNumDiscs];
614 
615 
616  hIp=new TH2D("Proj_to_IP","Proj_to_Ip",50,-100,100,50,-100,100);
617  hIpZAtX0=new TH1D("Proj_to_IPAtX0","Proj_toIPAtX0",50,-100,100);
618  hIpZAtY0=new TH1D("Proj_to_IPAtY0","Proj_toIPAtY0",50,-100,100);
619 
620  for(int iD=0;iD<kFgtNumDiscs;iD++)
621  {
622  // cout <<"id: " << iD <<endl;
623 
624  sprintf(buffer,"radioDiskEff_%d",iD);
625  radioPlotsEff[iD]=new TH2D(buffer,buffer,20,-50,50,20,-50,50);
626 
627  // cout <<"1" <<endl;
628  sprintf(buffer,"rEff_%d",iD);
629  rEff[iD]=new TH1D(buffer,buffer,100,0,50);
630  sprintf(buffer,"rNonEff_%d",iD);
631  rNonEff[iD]=new TH1D(buffer,buffer,100,0,50);
632  sprintf(buffer,"clusterSizeR_Disk_%d",iD);
633  h_clusterSizeR[iD]=new TH1D(buffer,buffer,20,0,20);
634  sprintf(buffer,"clusterSizePhi_Disk_%d",iD);
635  h_clusterSizePhi[iD]=new TH1D(buffer,buffer,20,0,20);
636  sprintf(buffer,"clusterChargeR_Disk_%d",iD);
637  h_clusterChargeR[iD]=new TH1D(buffer,buffer,100,0,50000);
638  sprintf(buffer,"clusterChargePhi_Disk_%d",iD);
639  h_clusterChargePhi[iD]=new TH1D(buffer,buffer,100,0,50000);
640  // cout <<"2" <<endl;
641  for(int nx=0;nx<radioPlotsEff[iD]->GetNbinsX();nx++)
642  {
643  for(int ny=0;ny<radioPlotsEff[iD]->GetNbinsY();ny++)
644  {
645  // radioPlotsEff[iD]->SetBinContent(nx,ny,0.1);//so that there is no divide by zero
646  }
647  }
648  // cout <<"3" <<endl;
649  sprintf(buffer,"r_phi_ChargeCorrelationInDisk_%d",iD);
650  chargeCorr[iD]=new TH2D(buffer,buffer,50,0,50000,50,0,50000);
651  sprintf(buffer,"radioDiskNonEff_%d",iD);
652  radioPlotsNonEff[iD]=new TH2D(buffer,buffer,20,-50,50,20,-50,50);
653  for(int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
654  {
655  // rEff[iD]->SetBinContent(nx,0.1);
656  }
657  // cout <<"4" <<endl;
658  sprintf(buffer,"rPhiRatio_%d",iD);
659  rPhiRatioPlots[iD]=new TH1D(buffer,buffer,100,-2,10);
660  }
661  return ierr;
662 };
663 ClassImp(StFgtAVEfficiencyMaker);
TDataSet::Clone
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone
Definition: TDataSet.cxx:308
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5