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"
22 Double_t getRPhiRatio(StSPtrVecFgtHitConstIterator hitIterBegin, StSPtrVecFgtHitConstIterator hitIterEnd)
24 Short_t quad, disc, strip;
26 Bool_t stripDead=
false;
30 StSPtrVecFgtHitConstIterator hitIter=hitIterBegin;
31 for(;hitIter!=hitIterEnd;hitIter++)
33 StFgtGeom::decodeGeoId((*hitIter)->getCentralStripGeoId(),disc, quad, layer, strip);
41 return (numR-numPhi)/((Double_t)(numR+numPhi));
53 eventPtr = (
StEvent*)GetInputDS(
"StEvent");
67 Double_t ratio=getRPhiRatio(tmpClusterCol->getHitVec().begin(),tmpClusterCol->getHitVec().end());
68 rPhiRatioPlots[i]->Fill(ratio);
72 const StSPtrVecFgtHit &hitVecD1=clusterColD1->getHitVec();
73 const StSPtrVecFgtHit &hitVecD6=clusterColD6->getHitVec();
76 set<Int_t> usedPoints;
77 Double_t D1Pos=StFgtGeom::getDiscZ(0);
78 Double_t D6Pos=StFgtGeom::getDiscZ(5);
79 Double_t zArm=D6Pos-D1Pos;
80 StSPtrVecFgtHitConstIterator hitIterD1,hitIterD6, hitIterD1R, hitIterD6R, hitIter, hitIter2;
81 cout <<
" there are " << hitVecD1.size() <<
" hits in d1 "<<endl;
82 for(hitIterD1=hitVecD1.begin();hitIterD1 != hitVecD1.end();hitIterD1++)
84 Short_t quad, disc, strip;
86 Bool_t stripDead=
false;
87 StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
89 cout <<
"found phi d1 " <<endl;
91 cout <<
"found r d1 " <<endl;
93 cout <<
" there are " << hitVecD6.size() <<
" hits in d6 "<<endl;
94 for(hitIterD1=hitVecD6.begin();hitIterD1 != hitVecD6.end();hitIterD1++)
96 Short_t quad, disc, strip;
98 Bool_t stripDead=
false;
99 StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
101 cout <<
"found phi d6 " <<endl;
103 cout <<
"found r d6 " <<endl;
107 for(
int iSeed1=0;iSeed1<5;iSeed1++)
109 for(
int iSeed2=iSeed1+1;iSeed2<6;iSeed2++)
113 const StSPtrVecFgtHit &hitVecSeed1=clusterColSeed1->getHitVec();
114 const StSPtrVecFgtHit &hitVecSeed2=clusterColSeed2->getHitVec();
115 D1Pos=StFgtGeom::getDiscZ(iSeed1);
116 D6Pos=StFgtGeom::getDiscZ(iSeed2);
121 for(hitIterD1=hitVecSeed1.begin();hitIterD1 != hitVecSeed1.end();hitIterD1++)
123 Short_t quad, disc, strip;
125 Bool_t stripDead=
false;
126 StFgtGeom::decodeGeoId((*hitIterD1)->getCentralStripGeoId(),disc, quad, layer, strip);
132 Int_t geoIdD1=(*hitIterD1)->getCentralStripGeoId();
133 StFgtGeom::decodeGeoId(geoIdD1,disc, quad, layer, strip);
137 Float_t phiD1=(*hitIterD1)->getPositionPhi();
138 for(hitIterD1R=hitVecSeed1.begin();hitIterD1R != hitVecSeed1.end();hitIterD1R++)
140 Int_t geoIdR=(*hitIterD1R)->getCentralStripGeoId();
141 StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
145 Float_t rD1=(*hitIterD1R)->getPositionR();
146 Float_t xD1=rD1*cos(phiD1);
147 Float_t yD1=rD1*sin(phiD1);
148 for(hitIterD6=hitVecSeed2.begin();hitIterD6 != hitVecSeed2.end();hitIterD6++)
150 Int_t geoIdD6=(*hitIterD6)->getCentralStripGeoId();
151 StFgtGeom::decodeGeoId(geoIdD6,disc, quad, layer, strip);
154 Float_t phiD6=(*hitIterD6)->getPositionPhi();
155 for(hitIterD6R=hitVecSeed2.begin();hitIterD6R != hitVecSeed2.end();hitIterD6R++)
157 Int_t geoIdR=(*hitIterD6R)->getCentralStripGeoId();
158 StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
162 Float_t rD6=(*hitIterD6R)->getPositionR();
163 Double_t xD6=rD6*cos(phiD6);
164 Double_t yD6=rD6*sin(phiD6);
167 vector< pair<Int_t,Double_t> > v_x;
168 vector< pair<Int_t,Double_t> > v_y;
169 vector< pair<Int_t,Double_t> > v_r;
171 vector< pair<Int_t,Double_t> > v_xFail;
172 vector< pair<Int_t,Double_t> > v_yFail;
173 vector< pair<Int_t,Double_t> > v_rFail;
175 vector< pair< Int_t, Double_t> > v_rCharge;
176 vector< pair< Int_t, Double_t> > v_phiCharge;
178 vector<Int_t> v_clusterSizeR;
179 vector<Int_t> v_clusterSizePhi;
181 vector<Int_t> v_geoIDsR;
182 vector<Int_t> v_geoIDsPhi;
188 Double_t xIpExp=xD1+(xD6-xD1)*(-D1Pos)/zArm;
189 Double_t yIpExp=yD1+(yD6-yD1)*(-D1Pos)/zArm;
191 Double_t zIpExpX0=(D6Pos-(xD6/xD1)*D1Pos)*1/(1-xD6/xD1);
193 Double_t zIpExpY0=(D6Pos-yD6/yD1*D1Pos)*1/(1-yD6/yD1);
196 for(
int iD=0;iD<6;iD++)
198 if(iD==iSeed1 || iD==iSeed2)
204 Double_t diskZ=StFgtGeom::getDiscZ(iD);
207 Double_t xPosExp=xD1+(xD6-xD1)*(diskZ-D1Pos)/zArm;
208 Double_t yPosExp=yD1+(yD6-yD1)*(diskZ-D1Pos)/zArm;
209 Double_t rPosExp=rD1+(rD6-rD1)*(diskZ-D1Pos)/zArm;
217 const StSPtrVecFgtHit &hitVec=clusterCol->getHitVec();
219 for(hitIter=hitVec.begin();hitIter!=hitVec.end();hitIter++)
222 Int_t geoIdPhi=(*hitIter)->getCentralStripGeoId();
223 StFgtGeom::decodeGeoId(geoIdPhi,disc, quad, layer, strip);
224 if(usedPoints.find(geoIdPhi)!=usedPoints.end())
228 Float_t phi=(*hitIter)->getPositionPhi();
229 Double_t phiCharge=(*hitIter)->charge();
230 Int_t clusterSizePhi=(*hitIter)->getStripWeightMap().size();
231 for(hitIter2=hitVec.begin();hitIter2!=hitVec.end();hitIter2++)
233 Int_t geoIdR=(*hitIter2)->getCentralStripGeoId();
234 StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);
235 if(usedPoints.find(geoIdR)!=usedPoints.end())
239 Float_t r=(*hitIter2)->getPositionR();
240 Double_t rCharge=(*hitIter2)->charge();
241 Int_t clusterSizeR=(*hitIter2)->getStripWeightMap().size();
243 if(fabs(r-rPosExp)<0.5)
247 v_r.push_back(pair<Int_t,Double_t> (iD,r));
253 if(fabs(x-xPosExp) < 0.5 && fabs(y-yPosExp)<0.5)
256 cout <<
"found! " <<
" pushing back: iD: " << iD <<
"x: " << x <<
" y "<< y <<endl;
257 v_x.push_back(pair<Int_t,Double_t>(iD,x));
258 v_y.push_back(pair<Int_t,Double_t>(iD,y));
259 v_rCharge.push_back(pair<Int_t, Double_t>(iD,rCharge));
260 v_phiCharge.push_back(pair<Int_t, Double_t>(iD,phiCharge));
262 v_geoIDsPhi.push_back(geoIdPhi);
263 v_geoIDsR.push_back(geoIdR);
265 v_clusterSizeR.push_back(clusterSizeR);
266 v_clusterSizePhi.push_back(clusterSizePhi);
276 v_xFail.push_back(pair<Int_t,Double_t>(iD,xPosExp));
277 v_yFail.push_back(pair<Int_t,Double_t>(iD,yPosExp));
284 v_rFail.push_back(pair<Int_t, Double_t>(iD,rPosExp));
292 hIp->Fill(xIpExp,yIpExp);
293 hIpZAtX0->Fill(zIpExpX0);
294 hIpZAtY0->Fill(zIpExpY0);
300 if(iFound>=2 && fabs(zIpExpX0)<100 && fabs(zIpExpY0)<100)
302 if(v_x.size()>iFound)
304 cout<<
"more hits than disks hit!!!" <<endl;
305 cout <<
" vx_size: " << v_x.size() <<
" ifound: " << iFound <<endl;
311 for(
int i=0;i<v_x.size();i++)
313 usedPoints.insert(v_geoIDsPhi[i]);
314 usedPoints.insert(v_geoIDsR[i]);
315 Int_t disk=v_x[i].first;
316 Double_t x=v_x[i].second;
317 Double_t y=v_y[i].second;
318 Double_t rCharge=v_rCharge[i].second;
319 Double_t phiCharge=v_phiCharge[i].second;
321 radioPlotsEff[disk]->Fill(x,y);
322 chargeCorr[disk]->Fill(rCharge,phiCharge);
323 h_clusterSizeR[disk]->Fill(v_clusterSizeR[i]);
324 h_clusterSizePhi[disk]->Fill(v_clusterSizePhi[i]);
325 h_clusterChargeR[disk]->Fill(rCharge);
326 h_clusterChargePhi[disk]->Fill(phiCharge);
328 for(
int i=0;i<v_xFail.size();i++)
330 Int_t disk=v_xFail[i].first;
331 Double_t x=v_xFail[i].second;
332 Double_t y=v_yFail[i].second;
334 radioPlotsNonEff[disk]->Fill(x,y);
342 if(v_r.size()>iFound)
348 for(
int i=0;i<v_r.size();i++)
350 Int_t disk=v_r[i].first;
351 Double_t r=v_r[i].second;
355 for(
int i=0;i<v_rFail.size();i++)
357 Int_t disk=v_rFail[i].first;
358 Double_t r=v_rFail[i].second;
360 rNonEff[disk]->Fill(r);
389 StFgtAVEfficiencyMaker::StFgtAVEfficiencyMaker(
const Char_t* name):runningEvtNr(0),hitCounter(0),hitCounterR(0)
391 cout <<
"AVE constructor!!" <<endl;
396 StFgtAVEfficiencyMaker::~StFgtAVEfficiencyMaker()
404 gStyle->SetPalette(1);
405 cout <<
"cluster plotter finish funciton " <<endl;
408 TCanvas* cRadio=
new TCanvas(
"radioPlots",
"radioPlot",1000,1500);
409 TCanvas* cRadioHits=
new TCanvas(
"radioPlotsHits",
"radioPlotHits",1000,1500);
410 TCanvas* cRadioNonHits=
new TCanvas(
"radioPlotsNonHits",
"radioPlotNonHits",1000,1500);
412 cRadioHits->Divide(2,3);
413 cRadioNonHits->Divide(2,3);
414 TCanvas* cRPRatio=
new TCanvas(
"rPhiRatio",
"rPhiRatios",1000,1500);
415 cRPRatio->Divide(2,3);
417 TCanvas* cREff=
new TCanvas(
"crEff",
"crEff",1000,1500);
421 for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
424 cRadioHits->cd(iD+1);
425 radioPlotsEff[iD]->Draw(
"colz");
426 cRadioNonHits->cd(iD+1);
427 radioPlotsNonEff[iD]->Draw(
"colz");
430 cRadioHits->SaveAs(
"radioPlotsHits.png");
431 cRadioNonHits->SaveAs(
"radioPlotsNonHits.png");
434 TCanvas* cClusterSizeR=
new TCanvas(
"clusterSizeR",
"clusterSizeR",1000,1500);
435 cClusterSizeR->Divide(2,3);
436 TCanvas* cClusterSizePhi=
new TCanvas(
"clusterSizePhi",
"clusterSizePhi",1000,1500);
437 cClusterSizePhi->Divide(2,3);
438 TCanvas* cChargeCorr=
new TCanvas(
"chargeCorr",
"chargeCorr",1000,1500);
439 cChargeCorr->Divide(2,3);
441 TCanvas* cClusterChargePhi=
new TCanvas(
"clusterChargePhi",
"clusterChargePhi",1000,1500);
442 cClusterChargePhi->Divide(2,3);
443 TCanvas* cClusterChargeR=
new TCanvas(
"clusterChargeR",
"clusterChargeR",1000,1500);
444 cClusterChargeR->Divide(2,3);
448 cIPProj.SaveAs(
"ipProj.png");
451 cIPProj.SaveAs(
"ipProjAtX0.png");
454 cIPProj.SaveAs(
"ipProjAtY0.png");
457 for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
459 cClusterSizeR->cd(iD+1);
460 h_clusterSizeR[iD]->Draw();
461 cClusterSizePhi->cd(iD+1);
462 h_clusterSizePhi[iD]->Draw();
463 cChargeCorr->cd(iD+1);
464 chargeCorr[iD]->Draw(
"colz");
466 cClusterChargeR->cd(iD+1);
467 h_clusterChargeR[iD]->Draw();
468 cClusterChargePhi->cd(iD+1);
469 h_clusterChargePhi[iD]->Draw();
472 cClusterSizeR->SaveAs(
"clusterSizeR.png");
473 cClusterSizePhi->SaveAs(
"clusterSizePhi.png");
474 cChargeCorr->SaveAs(
"chargeCorrelation.png");
476 cClusterChargeR->SaveAs(
"clusterChargeR.png");
477 cClusterChargePhi->SaveAs(
"clusterChargePhi.png");
479 cout <<
"saving .." <<endl;
481 for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
485 TH2D* tmpAllCounts=(TH2D*)radioPlotsEff[iD]->
Clone(
"tmp");
486 radioPlotsEff[iD]->Add(radioPlotsNonEff[iD]);
488 for(
int nx=0;nx<radioPlotsEff[iD]->GetNbinsX();nx++)
490 for(
int ny=0;ny<radioPlotsEff[iD]->GetNbinsY();ny++)
492 Double_t denom=radioPlotsEff[iD]->GetBinContent(nx,ny);
494 radioPlotsEff[iD]->SetBinContent(nx,ny,tmpAllCounts->GetBinContent(nx,ny)/denom);
496 radioPlotsEff[iD]->SetBinContent(nx,ny,0.0);
500 radioPlotsEff[iD]->Draw(
"colz");
501 TArc *innerArc =
new TArc(0,0,103.5+11.5,0.0,90.0);
502 TArc *outerArc =
new TArc(0,0,394.0-11.5,0.0,90.0);
503 TLine *left =
new TLine( 11.5, 103.5+11.5, 11.5, 394.0-11.5 );
504 TLine *right =
new TLine( 103.5+11.5, 11.5, 394.0-11.5, 11.5 );
505 TLine *center =
new TLine(
506 (103.5+11.5)*cos( 0.785398163 ),
507 (103.5+11.5)*sin( 0.785398163 ),
508 (394.0-11.5)*cos( 0.785398163 ),
509 (394.0-11.5)*sin( 0.785398163 )
511 TLine *weird =
new TLine(
512 (394.0-11.5)*cos( 0.190240888 ),
513 (394.0-11.5)*sin( 0.190240888 ),
514 (394.0-11.5)*cos( 0.891863248 ),
515 (394.0-11.5)*sin( 0.891863248 )
517 innerArc->SetFillStyle(0);
518 outerArc->SetFillStyle(0);
519 innerArc->SetLineWidth(2);
520 outerArc->SetLineWidth(2);
521 left->SetLineWidth(2);
522 right->SetLineWidth(2);
523 center->SetLineWidth(2);
524 weird->SetLineWidth(2);
525 outerArc->Draw(
"only");
526 innerArc->Draw(
"only");
537 rPhiRatioPlots[iD]->Draw();
540 TH1D* tmpR=(TH1D*)rEff[iD]->
Clone(
"tmpR");
541 rEff[iD]->Add(rNonEff[iD]);
542 for(
int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
544 Double_t denom=rEff[iD]->GetBinContent(nx);
546 rEff[iD]->SetBinContent(nx,tmpR->GetBinContent(nx)/denom);
548 rEff[iD]->SetBinContent(nx,0.0);
552 cRadio->SaveAs(
"radioPlotsEff.png");
553 cRadio->SaveAs(
"radioPlotsEff.pdf");
556 cREff->SaveAs(
"rEff.png");
557 cREff->SaveAs(
"rEff.pdf");
559 cRPRatio->SaveAs(
"rpRatio.png");
560 cRPRatio->SaveAs(
"rpRatio.pdf");
562 cout <<
"returning after finish" <<endl;
572 cout <<
"AVE!!" <<endl;
573 myRootFile=
new TFile(
"clusterEff.root",
"RECREATE");
582 radioPlotsEff=
new TH2D*[kFgtNumDiscs];
583 radioPlotsNonEff=
new TH2D*[kFgtNumDiscs];
584 rPhiRatioPlots=
new TH1D*[kFgtNumDiscs];
585 rEff=
new TH1D*[kFgtNumDiscs];
586 rNonEff=
new TH1D*[kFgtNumDiscs];
589 chargeCorr=
new TH2D*[kFgtNumDiscs];
590 h_clusterSizeR=
new TH1D*[kFgtNumDiscs];
591 h_clusterSizePhi=
new TH1D*[kFgtNumDiscs];
593 h_clusterChargeR=
new TH1D*[kFgtNumDiscs];
594 h_clusterChargePhi=
new TH1D*[kFgtNumDiscs];
597 hIp=
new TH2D(
"Proj_to_IP",
"Proj_to_Ip",50,-100,100,50,-100,100);
598 hIpZAtX0=
new TH1D(
"Proj_to_IPAtX0",
"Proj_toIPAtX0",50,-100,100);
599 hIpZAtY0=
new TH1D(
"Proj_to_IPAtY0",
"Proj_toIPAtY0",50,-100,100);
601 for(
int iD=0;iD<kFgtNumDiscs;iD++)
605 sprintf(buffer,
"radioDiskEff_%d",iD);
606 radioPlotsEff[iD]=
new TH2D(buffer,buffer,20,-50,50,20,-50,50);
609 sprintf(buffer,
"rEff_%d",iD);
610 rEff[iD]=
new TH1D(buffer,buffer,100,0,50);
611 sprintf(buffer,
"rNonEff_%d",iD);
612 rNonEff[iD]=
new TH1D(buffer,buffer,100,0,50);
613 sprintf(buffer,
"clusterSizeR_Disk_%d",iD);
614 h_clusterSizeR[iD]=
new TH1D(buffer,buffer,20,0,20);
615 sprintf(buffer,
"clusterSizePhi_Disk_%d",iD);
616 h_clusterSizePhi[iD]=
new TH1D(buffer,buffer,20,0,20);
617 sprintf(buffer,
"clusterChargeR_Disk_%d",iD);
618 h_clusterChargeR[iD]=
new TH1D(buffer,buffer,100,0,50000);
619 sprintf(buffer,
"clusterChargePhi_Disk_%d",iD);
620 h_clusterChargePhi[iD]=
new TH1D(buffer,buffer,100,0,50000);
622 for(
int nx=0;nx<radioPlotsEff[iD]->GetNbinsX();nx++)
624 for(
int ny=0;ny<radioPlotsEff[iD]->GetNbinsY();ny++)
630 sprintf(buffer,
"r_phi_ChargeCorrelationInDisk_%d",iD);
631 chargeCorr[iD]=
new TH2D(buffer,buffer,50,0,50000,50,0,50000);
632 sprintf(buffer,
"radioDiskNonEff_%d",iD);
633 radioPlotsNonEff[iD]=
new TH2D(buffer,buffer,20,-50,50,20,-50,50);
634 for(
int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
639 sprintf(buffer,
"rPhiRatio_%d",iD);
640 rPhiRatioPlots[iD]=
new TH1D(buffer,buffer,100,-2,10);
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone