StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtQAMaker.cxx
1 
11 #include "StFgtQAMaker.h"
12 #include "StEventTypes.h"
13 #include "StMessMgr.h"
14 #include "StDcaGeometry.h"
15 #include "StThreeVectorF.hh"
16 #include "StFgtQAMaker.h"
17 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
18 #include "StRoot/StFgtDbMaker/StFgtDbMaker.h"
19 #include "StRoot/StFgtDbMaker/StFgtDb.h"
20 #include "StRoot/StFgtPool/StFgtClusterTools/StFgtGeneralBase.h"
21 #include "StRoot/StFgtPool/StFgtClusterTools/StFgtStraightTrackMaker.h"
22 #include "StarClassLibrary/StThreeVectorF.hh"
23 
24 #include <TMath.h>
25 #include <TGraphAsymmErrors.h>
26 
27 static const int mDebug=0; // debug mesasge level
28 static const int NTB=15;
29 static const char* cquad[kFgtNumQuads]={"A","B","C","D"};
30 
31 inline int globalapv(int rdo, int arm, int apv){
32  return (rdo-1)*6*2 + arm*2 + apv/12;
33  // if(apv<10) { return (rdo-1)*6*2*10 + arm*2*10 + apv;}
34  //else { return (rdo-1)*6*2*10 + arm*2*10 + apv - 2;}
35 }
36 
37 inline void getPRC(int rdo, int arm, int apv, int& page, int& row, int& col){
38  page = ((rdo-1)*6 + arm)/2;
39  row = (((rdo-1)*6 + arm)%2)*4 + (apv/12)*2 + (apv%12)/5;
40  col = (apv%12)%5;
41  //printf("rdo=%1d arm=%1d apv=%2d page=%d row=%d col=%d\n",rdo,arm,apv,page,row,col);
42 }
43 
44 inline double localphi(int quad, double phi){
45  static const double mPi = TMath::Pi();
46  static const double phioff[kFgtNumQuads] = {-15.0*mPi/180.0, -105.0*mPi/180.0, 165.0*mPi/180.0, 165.0*mPi/180.0};
47  double local=phi - phioff[quad];
48  while(local<0.0) {local+=2.0*mPi;}
49  while(local>2.0*mPi) {local-=2.0*mPi;}
50  return local;
51 }
52 
53 ClassImp(StFgtQAMaker);
54 
55 StFgtQAMaker::StFgtQAMaker(const Char_t *name) : StMaker(name),mEventCounter(0),mRunNumber(0),mSigmaCut(10),mAdcCut(500){
56 }
57 
58 Int_t StFgtQAMaker::Init(){
59  bookHist();
60  for(int i=0; i<7; i++){
61  char c[40];
62  sprintf(c,"QA%d",i);
63  mCanvas[i] = new TCanvas(c,c,50,0,700,800);
64  TH2F* frame;
65  sprintf(c,"ScopeTrace%d",i);
66  if(i==6) { frame = new TH2F(c,c,1,0,NTB*4,1,-1000,5000*6); frame->GetXaxis()->SetNdivisions(NTB*100+4,kFALSE);}
67  else { frame = new TH2F(c,c,1,0,NTB*5,1,-1000,5000*8); frame->GetXaxis()->SetNdivisions(NTB*100+5,kFALSE);}
68  frame->SetStats(0); frame->Draw();
69  }
70  memset(mNTrace,0,sizeof(mNTrace));
71  memset(mNTrace2,0,sizeof(mNTrace2));
72  return kStOK;
73 }
74 
75 Int_t StFgtQAMaker::InitRun(Int_t runnum){
76  LOG_INFO << "StFgtQAMaker::InitRun for " << runnum << endm;
77  StFgtDbMaker *fgtDbMkr = static_cast<StFgtDbMaker * >( GetMaker("fgtDb"));
78  if( !fgtDbMkr ){
79  LOG_FATAL << "StFgtDb not provided and error finding StFgtDbMaker" << endm;
80  return kStFatal;
81  } else {
82  mDb = fgtDbMkr->getDbTables();
83  if( !mDb ){
84  LOG_FATAL << "StFgtDb not provided and error retrieving pointer from StFgtDbMaker '"
85  << fgtDbMkr->GetName() << endm;
86  return kStFatal;
87  }
88  }
89  return kStOK;
90 }
91 
92 Int_t StFgtQAMaker::Make() {
93  fillHist();
94  mEventCounter++; //if(mEventCounter%100==0) printf("Event=%d\n",mEventCounter);
95  return kStOK;
96 }
97 
98 Int_t StFgtQAMaker::Finish() {
99  gMessMgr->Info() << "StFgtQAMaker::Finish()" << endm;
100  saveHist();
101  saveTrace();
102  return kStOK;
103 }
104 
105 void StFgtQAMaker::saveTrace(){
106  char cq[4][3]={"DA","AB","BC","CD"};
107  char c[20];
108 
109  char cthr[100];
110  sprintf(cthr,"MaxADC>%2.0f*sigma and %4.0f",mSigmaCut,mAdcCut);
111  TText *tthr = new TText(0.6,0.95,cthr);
112  tthr->SetNDC(); tthr->SetTextSize(0.03);
113 
114  for(int rdo=1; rdo<=kFgtNumRdos; rdo++){
115  for(int arm=0; arm<kFgtNumArms; arm++){
116  for(int apv=0; apv<kFgtMaxApvId; apv++){
117  if(apv==10 || apv==11 || apv==22 || apv==23) continue;
118  int page,row,col;
119  getPRC(rdo,arm,apv,page,row,col);
120  mCanvas[page]->cd();
121 
122  sprintf(c,"Apv%1d",apv);
123  float x = 2 + col*NTB;
124  float y = 4000 + (7-row)*5000;
125  TText *t1 = new TText(x,y,c);
126  t1->SetTextSize(0.02); t1->Draw();
127 
128  if(apv==9|| apv==21){
129  Short_t disc, quad, strip;
130  Char_t layer;
131  int geoid=mDb->getGeoIdFromElecCoord(rdo,arm,apv-9,0);
132  StFgtGeom::decodeGeoId(geoid,disc,quad,layer,strip);
133  sprintf(c,"Rdo%1dArm%1d/%1d%2s",rdo,arm,disc+1,cq[quad]);
134  float x = (col+1.2)*NTB;
135  float y = 1000 + (7-row)*5000;
136  TText *t1 = new TText(x,y,c);
137  t1->SetTextSize(0.022); t1->SetTextAngle(90); t1->Draw();
138  }
139  }
140  }
141  }
142  for(int i=0; i<6; i++){
143  mCanvas[i]->cd();
144  tthr->Draw();
145  char ff[50],fname[50]="fgtScopeTrace.pdf";
146  if(mRunNumber>0) {
147  sprintf(fname,"%d/fgtScopeTrace_%d.pdf",mRunNumber/1000,mRunNumber);
148  }
149  if(i==0) {
150  cout << "Writing " << fname << endl;
151  sprintf(ff,"%s(",fname);
152  }else if(i==5){
153  sprintf(ff,"%s)",fname);
154  }else{
155  sprintf(ff,"%s",fname);
156  }
157  mCanvas[i]->Update();
158  mCanvas[i]->Print(ff,"pdf");
159  }
160  mCanvas[6]->cd();
161  tthr->Draw();
162  for(int disc=0; disc<kFgtNumDiscs; disc++){
163  sprintf(c,"Disc%1d",disc+1);
164  float y=0.82-disc*0.13;
165  TText *t1 = new TText(0.91,y,c);
166  t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw();
167  }
168  for(int quad=0; quad<kFgtNumQuads; quad++){
169  sprintf(c,"Quad%1c",'A'+quad);
170  float x=0.15 + quad*0.20;
171  TText *t1 = new TText(x,0.91,c);
172  t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw();
173  }
174  char fname[50]="fgtScopeTrace.png";
175  if(mRunNumber>0) {
176  sprintf(fname,"%d/fgtScopeTrace_%d.png",mRunNumber/1000,mRunNumber);
177  }
178  cout << "Writing " << fname << endl;
179  mCanvas[6]->Update();
180  mCanvas[6]->SaveAs(fname);
181 }
182 
183 void StFgtQAMaker::bookHist(){
184  char c[50];
185  int ns=kFgtNumStrips;
186  const int nbin= 512;
187  const int max=30720;
188  const int binwid=max/nbin;
189  const float min=0.5-binwid;
190  const float max2=max+0.5;
191  const char* cHist[NHist]={"MaxTimeBin","ADC","DataSize", "ZSdata", "10sigma", "Nevt", "LandauChi2","LandauSig","LandauMpv", "Mpv-3Sig"};
192  const int nHist[NHist]={ NTB, 200, nbin+1, kFgtNumElecIds, kFgtNumElecIds, 1, 40, 40, 60, 60};
193  const float lHist[NHist]={ 0, 0, min, 0, 0, 0, 0, 0, 0, -3};
194  const float hHist[NHist]={ float(NTB), 4200, max2, kFgtNumElecIds, kFgtNumElecIds, 1, 20, 4, float(NTB),float(NTB-3)};
195  const char* c1dHist[N1dHist]={"NHitStrip", "PhiHit", "RHit", "NCluster","ClusterSize","ClusterCharge","MaxADC","ChargeAsy","CluChargeT","MaxADCT","ChargeAsyTrk","LandauN"};
196  const int n1dHist[N1dHist]={ 50, ns, ns, 25, 15, 50, 50, 25, 50, 50, 25, 50};
197  const float l1dHist[N1dHist]={ 0, 0, 0, 0, 0, 0, 0, -0.3, 0, 0, -0.3, 0};
198  const float h1dHist[N1dHist]={ 100,float(ns),float(ns), 25.0, 15, 30000, 4000, 0.3, 30000, 4000, 0.3, 30000};
199  const char* c2dHist[N2dHist] ={ "XY", "ADCvsTB", "APVTB", "XYT"};
200  const int xn2dHist[N2dHist]={ 50, NTB, 60, 50};
201  const float xl2dHist[N2dHist]={ -40, 0, 0, -40};
202  const float xh2dHist[N2dHist]={ 40,float(NTB),float(NTB), 40};
203  const int yn2dHist[N2dHist]={ 50, 200, 24, 50};
204  const float yl2dHist[N2dHist]={ -40, 0, 0, -40};
205  const float yh2dHist[N2dHist]={ 40, 4000, 24, 40};
206  const char* cTHist[NTrkHist]={ "NTrk","NHitTrk","DCA","ZVTX","Chi2","HitDisc"};
207  const int nTHist[NTrkHist]={ 10, 7, 60, 50, 50, 6};
208  const float lTHist[NTrkHist]={ 0, 0, 0, -100, 0, 0.5};
209  const float hTHist[NTrkHist]={ 10, 7, 4, 100, 0.015, 6.5};
210  //histos for whole FGT
211  for(int i=0; i<NHist; i++){
212  hist0[i]=new TH1F(cHist[i],cHist[i],nHist[i],lHist[i],hHist[i]);
213  }
214  //1d histos per disc/quad
215  for(int disc=0; disc<kFgtNumDiscs; disc++){
216  for(int quad=0; quad<kFgtNumQuads; quad++){
217  for(int i=0; i<N1dHist; i++){
218  sprintf(c,"%1d%1s-%s",disc+1,cquad[quad],c1dHist[i]);
219  hist1[disc][quad][i]=new TH1F(c,c,n1dHist[i],l1dHist[i],h1dHist[i]);
220  }
221  }
222  }
223  //2d histos per disc
224  for(int disc=0; disc<kFgtNumDiscs; disc++){
225  for(int i=0; i<N2dHist; i++){
226  sprintf(c,"Disc%1d%s",disc+1,c2dHist[i]);
227  hist2[disc][i]=new TH2F(c,c,xn2dHist[i],xl2dHist[i],xh2dHist[i],yn2dHist[i],yl2dHist[i],yh2dHist[i]);
228  }
229  }
230  //1d histos per quad
231  for(int quad=0; quad<kFgtNumQuads; quad++){
232  for(int i=0; i<NTrkHist; i++){
233  sprintf(c,"Quad%1s-%s",cquad[quad],cTHist[i]);
234  histTrk[quad][i]=new TH1F(c,c,nTHist[i],lTHist[i],hTHist[i]);
235  }
236  }
237 }
238 
239 void StFgtQAMaker::saveHist(){
240  char fname[50]="fgtQA.root";
241  if(mRunNumber>0) {
242  sprintf(fname,"%d/fgtQA_%d.root",mRunNumber/1000,mRunNumber);
243  }
244  cout << "Writing " << fname << endl;
245  TFile *hfile = new TFile(fname,"RECREATE");
246  for(int i=0; i<NHist; i++){
247  hist0[i]->Write();
248  }
249  //1d histos per disc/quad
250  for(int disc=0; disc<kFgtNumDiscs; disc++){
251  for(int quad=0; quad<kFgtNumQuads; quad++){
252  for(int i=0; i<N1dHist; i++){
253  hist1[disc][quad][i]->Write();
254  }
255  }
256  }
257  //2d histos per disc
258  for(int disc=0; disc<kFgtNumDiscs; disc++){
259  for(int i=0; i<N2dHist; i++){
260  hist2[disc][i]->Write();
261  }
262  }
263  //1d histos per quad
264  for(int quad=0; quad<kFgtNumQuads; quad++){
265  for(int i=0; i<NTrkHist; i++){
266  histTrk[quad][i]->Write();
267  }
268  }
269  hfile->Close();
270 }
271 
272 void StFgtQAMaker::fillHist(){
273  StEvent* eventPtr = (StEvent*)GetInputDS("StEvent");
274  fgtCollectionPtr=0;
275  if( eventPtr ) {
276  //LOG_Info << "Found StFgtCollection from StEvent in '" << ClassName() << "'" << endm;
277  fgtCollectionPtr = eventPtr->fgtCollection();
278  }else{
279  TObjectSet *os = (TObjectSet*)GetDataSet("FGTCOLLECTION");
280  if (os) {
281  fgtCollectionPtr = (StFgtCollection*)os->GetObject();
282  //LOG_Info << "Found StFgtCollection from FGTCOLLECTION in '" << ClassName() << "'" << endm;
283  }
284  }
285  if( !fgtCollectionPtr) {
286  LOG_ERROR << "Error getting pointer to StFgtCollection in '" << ClassName() << "'" << endm;
287  return;
288  };
289 
290  mNTimeBin=fgtCollectionPtr->getNumTimeBins();
291  int datasize=0;
292  if(mNTimeBin==0) {hist0[2]->Fill(float(datasize)); return; }
293  for (int disc=0; disc<kFgtNumDiscs; disc++){
294  //Strips
295  int nhit[kFgtNumQuads]; memset(nhit,0,sizeof(nhit));
296  StFgtStripCollection *cstrip = fgtCollectionPtr->getStripCollection(disc);
297  StSPtrVecFgtStrip &strip=cstrip->getStripVec();
298  for(StSPtrVecFgtStripIterator it=strip.begin(); it!=strip.end(); it++) {
299  datasize++;
300  int geoid =(*it)->getGeoId();
301  int maxadc =(*it)->getMaxAdc();
302  float ped =(*it)->getPed();
303  float pederr=(*it)->getPedErr();
304  char layer;
305  short idisc,iquad,ipr,istr;
306  StFgtGeom::decodeGeoId(geoid,idisc,iquad,layer,istr);
307  if(layer=='P') {ipr=0;}
308  if(layer=='R') {ipr=1;}
309  Int_t rdo, arm, apv, chan;
310  (*it)->getElecCoords(rdo, arm, apv, chan);
311  int gapv = globalapv(rdo,arm,apv);
312  int eid = StFgtGeom::encodeElectronicId(rdo, arm, apv, chan);
313  hist0[3]->Fill(float(eid));
314  if(maxadc>mSigmaCut*pederr && maxadc>mAdcCut){
315  hist0[4]->Fill(float(eid));
316  if(mNTrace2[rdo-1][arm][apv]<MAXTRACE) {
317  //printf("ntrace %d %d %d %d\n",rdo,arm,apv,mNTrace2[rdo-1][arm][apv]);
318  mNTrace2[rdo-1][arm][apv]++;
319  TGraph* g = new TGraph();
320  int page,row,col;
321  getPRC(rdo,arm,apv,page,row,col);
322  for(int tb=0; tb<mNTimeBin; tb++){
323  g->SetPoint(tb,float(tb+col*NTB),float((*it)->getAdc(tb)+(7-row)*5000));
324  }
325  mCanvas[page]->cd();
326  g->SetLineColor(kBlue); g->SetLineWidth(0.3); g->Draw("C");
327  }
328  if(mNTrace[idisc][iquad]<MAXTRACE) {
329  mNTrace[idisc][iquad]++;
330  TGraph* g = new TGraph();
331  for(int tb=0; tb<mNTimeBin; tb++){
332  g->SetPoint(tb,float(tb+iquad*NTB),float((*it)->getAdc(tb)+25000-idisc*5000));
333  }
334  mCanvas[6]->cd();
335  g->SetLineColor(kBlue); g->SetLineWidth(0.5); g->Draw("C");
336  }
337  }
338  hist0[1]->Fill(float(maxadc)+ped);
339  if(maxadc>4*pederr && maxadc>160) {
340  nhit[iquad]++;
341  if(ipr==1) {hist1[disc][iquad][1+ipr]->Fill(float(istr));}
342  else {hist1[disc][iquad][1+ipr]->Fill(float(istr/2 + (istr%2)*360));}
343  }
344  if(maxadc>mSigmaCut*pederr && maxadc>mAdcCut){
345  float max=0;
346  int maxtb=0;
347  for(int tb=0; tb<mNTimeBin; tb++){
348  float a=float((*it)->getAdc(tb));
349  if(a>200) hist2[idisc][1]->Fill(float(tb),a);
350  if(a>max) {max=a; maxtb=tb;}
351  hist2[0][2]->Fill(float(tb), float(gapv), a);
352  }
353  hist2[1][2]->Fill(float(maxtb), float(gapv));
354  }
355  }
356  for (int quad=0; quad<kFgtNumQuads; quad++) hist1[disc][quad][0]->Fill(float(nhit[quad]));
357  }
358  hist0[2]->Fill(float(datasize));
359  if(datasize>0) hist0[5]->Fill(0);
360  if(mEventCounter==100 || mEventCounter==500 || mEventCounter==1000){
361  printf("Event=%d Updating canvas...\n",mEventCounter); mCanvas[6]->Update(); printf("...Done\n");
362  }
363 
364  //Clusters
365  for (int disc=0; disc<kFgtNumDiscs; disc++){
366  int ncluster[kFgtNumQuads]; memset(ncluster,0,sizeof(ncluster));
367  StFgtHitCollection *chit = fgtCollectionPtr->getHitCollection(disc);
368  StSPtrVecFgtHit &hit=chit->getHitVec();
369  for(StSPtrVecFgtHitIterator it=hit.begin(); it!=hit.end(); it++) {
370  int geoid=(*it)->getCentralStripGeoId();
371  char layer;
372  short idisc,iquad,istr;
373  StFgtGeom::decodeGeoId(geoid,idisc,iquad,layer,istr);
374  ncluster[iquad]++;
375  hist1[disc][iquad][4]->Fill(float((*it)->getNstrip()));
376  hist1[disc][iquad][5]->Fill((*it)->charge());
377  hist1[disc][iquad][6]->Fill((*it)->getMaxAdc());
378  if((*it)->getMaxAdc() > 500){
379  float tb = float((*it)->getMaxTimeBin());
380  //float tb = (*it)->getLandauPeak();
381  hist0[0]->Fill(tb);
382  int rdo,arm,apv,ch;
383  mDb->getElecCoordFromGeoId(geoid,rdo,arm,apv,ch);
384  int gapv=globalapv(rdo,arm,apv);
385  hist2[2][2]->Fill(tb, float(gapv));
386  hist2[3][2]->Fill((*it)->getLandauMpv(), float(gapv));
387  hist2[4][2]->Fill((*it)->getLandauMpv()-(*it)->getLandauSigma()*3, float(gapv));
388  }
389  }
390  for (int quad=0; quad<kFgtNumQuads; quad++) hist1[disc][quad][3]->Fill(float(ncluster[quad]));
391  }
392  //Points
393  StFgtPointCollection *cpoint = fgtCollectionPtr->getPointCollection();
394  StSPtrVecFgtPoint &point = cpoint->getPointVec();
395  for(StSPtrVecFgtPointIterator it=point.begin(); it!=point.end(); it++) {
396  int idisc=(*it)->getDisc();
397  int iquad=(*it)->getQuad();
398  float phi = (*it)->getPositionPhi();
399  float r = (*it)->getPositionR();
400  TVector3 xyz;
401  mDb->getStarXYZ(idisc,iquad,r,phi,xyz);
402  hist2[idisc][0]->Fill(xyz.X(),xyz.Y());
403  hist1[idisc][iquad][7]->Fill((*it)->getChargeAsymmetry());
404  }
405 
406  //Tracks
407  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
408  if (fgtSTracker){
409  vector<AVTrack>& tracks=fgtSTracker->getTracks();
410  int ntrk[kFgtNumQuads]; memset(ntrk,0,sizeof(ntrk));
411  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
412  if(mDebug>0) cout<<Form("Trk chi2=%8.3f dca=%8.3f",t->chi2,t->dca)<<endl;
413  if(t->dca>3 || t->chi2>0.01) continue;
414  vector<AVPoint>* points=t->points;
415  int quad=-1,disc=-1,nhit=0;
416  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
417  nhit++;
418  disc=p->dID;
419  quad = p->quadID;
420  histTrk[quad][5]->Fill(float(disc+1));
421  hist1[disc][quad][8]->Fill(p->rCharge);
422  hist1[disc][quad][8]->Fill(p->phiCharge);
423  hist1[disc][quad][9]->Fill(p->fgtHitR->getMaxAdc());
424  hist1[disc][quad][9]->Fill(p->fgtHitPhi->getMaxAdc());
425  hist1[disc][quad][10]->Fill((p->phiCharge-p->rCharge)/(p->phiCharge+p->rCharge));
426  hist1[disc][quad][11]->Fill(p->fgtHitPhi->getLandauNorm());
427  hist1[disc][quad][11]->Fill(p->fgtHitR->getLandauNorm());
428  hist2[disc][3]->Fill(p->x,p->y);
429  StFgtHit *phi=p->fgtHitPhi;
430  StFgtHit *r=p->fgtHitR;
431  hist0[6]->Fill(phi->getLandauChi2()); hist0[6]->Fill(r->getLandauChi2());
432  hist0[7]->Fill(phi->getLandauSigma()); hist0[7]->Fill(r->getLandauSigma());
433  hist0[8]->Fill(phi->getLandauMpv()); hist0[8]->Fill(r->getLandauMpv());
434  hist0[9]->Fill(phi->getLandauMpv()- phi->getLandauSigma()*3);
435  hist0[9]->Fill(r->getLandauMpv() - r->getLandauSigma()*3);
436  pulseFit(p->fgtHitPhi); pulseFit(p->fgtHitR);
437  }
438  histTrk[quad][1]->Fill(nhit);
439  histTrk[quad][2]->Fill(t->dca);
440  histTrk[quad][3]->Fill(t->trkZ);
441  histTrk[quad][4]->Fill(t->chi2);
442  ntrk[quad]++;
443  }
444  for(int quad=0; quad<kFgtNumQuads; quad++) {histTrk[quad][0]->Fill(float(ntrk[quad]));}
445  }
446  textDump();
447 }
448 
449 void StFgtQAMaker::pulseFit(StFgtHit* cluster){
450  static TCanvas *C1=0;
451  static int NPLOT=0;
452  static int NPAGE=0;
453  static int MAXPAGE=10;
454  if(NPAGE>MAXPAGE) return;
455 
456  if(C1==0) C1=new TCanvas("Pulse","Pulse",50,50,750,800);
457  if(NPLOT==0){
458  C1->Clear();
459  C1->Divide(3,3);
460  }
461  TVirtualPad *pad = C1->cd(NPLOT+1);
462 
463  int geoid=cluster->getCentralStripGeoId();
464  char layer;
465  short disc,quad,str,pr;
466  StFgtGeom::decodeGeoId(geoid,disc,quad,layer,str);
467  if(layer=='P') {pr=0;}
468  if(layer=='R') {pr=1;}
469 
470  int nstrip=0, center=0;
471  StFgtStrip* strips[20];
472  StFgtStripCollection *cstrip = fgtCollectionPtr->getStripCollection(disc);
473  StSPtrVecFgtStrip &strip=cstrip->getStripVec();
474  for(StSPtrVecFgtStripIterator it=strip.begin(); it!=strip.end(); it++) {
475  int igeoid =(*it)->getGeoId();
476  char ilayer;
477  short idisc,iquad,ipr,istr;
478  StFgtGeom::decodeGeoId(igeoid,idisc,iquad,ilayer,istr);
479  if(ilayer=='P') {ipr=0;}
480  if(ilayer=='R') {ipr=1;}
481  if(geoid==igeoid) {center=nstrip;}
482  if(quad==iquad && pr==ipr && abs(str-istr)<10){ strips[nstrip]=(*it); nstrip++; }
483  }
484  float a[NTB]; memset(a,0,sizeof(a));
485  float e[NTB]; memset(e,0,sizeof(e));
486  int f[NTB]; memset(f,0,sizeof(f));
487  int nstr=0, nskip=0;
488  for(int i=center; i<nstrip; i++){
489  short type = strips[i]->getClusterSeedType();
490  if(type<kFgtSeedType1 || type>=kFgtClusterTooBig) nskip++;
491  if(nskip==2) break;
492  if(nskip==1) continue;
493  nstr++;
494  float ped = strips[i]->getPed();
495  float rms = strips[i]->getPedErr();
496  for(int t=0; t<NTB; t++){
497  int adc = strips[i]->getAdc(t);
498  a[t]+=adc;
499  e[t]+=rms*rms;
500  if(adc+ped>3200){f[t]=1;}// printf("Saturation??? %f\n",adc+ped);}
501  }
502  }
503  nskip=0;
504  for(int i=center-1; i>=0; i--){
505  short type = strips[i]->getClusterSeedType();
506  if(type<kFgtSeedType1 || type>=kFgtClusterTooBig) nskip++;
507  if(nskip==2) break;
508  if(nskip==1) continue;
509  nstr++;
510  float ped = strips[i]->getPed();
511  float rms = strips[i]->getPedErr();
512  for(int t=0; t<NTB; t++){
513  float adc = strips[i]->getAdc(t);
514  a[t]+=adc;
515  e[t]+=rms*rms;
516  if(adc+ped>3400){f[t]=1;}
517  }
518  }
519  TGraphAsymmErrors *g1 = new TGraphAsymmErrors(mNTimeBin);
520  TGraphAsymmErrors *g2 = new TGraphAsymmErrors(7);
521  TGraphAsymmErrors *g3 = new TGraphAsymmErrors(6);
522  //printf("%1d nstrip=%d adc=",NPLOT,nstrip);
523  for(int t=0; t<mNTimeBin; t++){
524  //printf("%4.0f(%4.0f) ",a[t],e[t]);
525  e[t]=sqrt(e[t]);
526  g1->SetPoint(t,float(t),a[t]);
527  if(f[t]==0) {g1->SetPointError(t,0,0,e[t],e[t]);}
528  else {g1->SetPointError(t,0,0,e[t],2000);}
529  }
530  //printf("\n");
531  for(int t=0; t<7; t++){
532  g2->SetPoint(t,float(t+2),a[t+2]);
533  if(f[t]==0) {g2->SetPointError(t,0,0,e[t+2],e[t+2]);}
534  else {g2->SetPointError(t,0,0,e[t+2],2000);}
535  }
536  for(int t=0; t<6; t++){
537  g3->SetPoint(t,float(t+2),a[t+2]);
538  if(f[t]==0) {g3->SetPointError(t,0,0,e[t+2],e[t+2]);}
539  else {g3->SetPointError(t,0,0,e[t+2],2000);}
540  }
541  g1->SetMarkerStyle(20); g1->SetMarkerSize(1);
542  g1->Draw("AP");
543 
544  gStyle->SetStatW(0.4);
545  int res1=g1->Fit("landau","Q0");
546  TF1 *f1 = g1->GetFunction("landau"); f1->SetLineColor(2); f1->SetLineWidth(2); f1->Draw("same");
547 
548  int res2=g2->Fit("landau","Q0");
549  TF1 *f2 = g2->GetFunction("landau"); f2->SetLineColor(4); f2->SetLineWidth(2); f2->Draw("same");
550 
551  int res3=g3->Fit("landau","Q0");
552  TF1 *f3 = g3->GetFunction("landau"); f3->SetLineColor(7); f3->SetLineWidth(2); f3->Draw("same");
553 
554  NPLOT++;
555  if(NPLOT==9){
556  C1->Update();
557  //char tmp[100];
558  char file[100];
559  if(NPAGE==0){
560  sprintf(file,"%d/%d.pulsefit.pdf(",int(mRunNumber)/1000,mRunNumber);
561  }else if(NPAGE==MAXPAGE){
562  sprintf(file,"%d/%d.pulsefit.pdf)",int(mRunNumber)/1000,mRunNumber);
563  }else{
564  sprintf(file,"%d/%d.pulsefit.pdf",int(mRunNumber)/1000,mRunNumber);
565  }
566  C1->Print(file,"pdf");
567  //printf("Press key!!!! -->"); cin >> tmp;
568  NPLOT=0;
569  NPAGE++;
570  }
571 }
572 
573 void StFgtQAMaker::dip(){
574 }
575 
576 void StFgtQAMaker::textDump(){
577  static int MAXEVT=100;
578  static int NEVT=0;
579  static FILE *mTextFile=0;
580  if(NEVT==0){
581  char filename[100];
582  sprintf(filename,"%d/%d.evtdump.txt",int(mRunNumber)/1000,mRunNumber);
583  printf("Writing %s\n",filename);
584  mTextFile = fopen(filename,"w+");
585  if(!mTextFile) {printf("Couldn't open %s\n",filename); NEVT=MAXEVT+1;}
586  }
587  //if(NEVT==MAXEVT) {fclose(mTextFile);}
588  if(NEVT>=MAXEVT) {return;}
589  fprintf(mTextFile,"===== Event %d =====\n",NEVT);
590  char layer[2]={'P','R'};
591  StEvent* eventPtr = (StEvent*)GetInputDS("StEvent");
592  if( !eventPtr ) { return; }
593  StFgtCollection* fgtCollectionPtr = eventPtr->fgtCollection();
594  if( !fgtCollectionPtr) { return; }
595  StFgtStraightTrackMaker *fgtSTracker = static_cast<StFgtStraightTrackMaker * >( GetMaker("fgtStraightTracker"));
596  if (! fgtSTracker){ return; }
597  vector<AVTrack>& tracks=fgtSTracker->getTracks();
598  int missedhit[kFgtNumGeoIds]; memset(missedhit,0,sizeof(missedhit));
599 
600  for(int quad=0; quad<kFgtNumQuads; quad++){
601  //print tracks
602  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
603  vector<AVPoint>* points=t->points;
604  int nhit=0, quad2=-1;
605  int dischit[kFgtNumDiscs]; memset(dischit,0,sizeof(dischit));
606  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
607  quad2=p->quadID;
608  if(quad2 != quad) break;
609  int disc=p->dID;
610  dischit[disc]=1;
611  if(nhit==0) fprintf(mTextFile,"TRK quad=%1s chi2=%6.3f dca=%6.2f vtxz=%6.2f\n",
612  cquad[quad],t->chi2,t->dca,t->trkZ);
613  fprintf(mTextFile,"TRK quad=%1s hit=%1d disc=%1d Pchg=%6.0f Rchg=%6.0f x=%6.2f y=%6.2f\n",
614  cquad[quad], nhit, disc+1, p->fgtHitPhi->getLandauNorm(),p->fgtHitR->getLandauNorm(),p->x,p->y);
615  //printf("TRK quad=%1s hit=%1d disc=%1d Pchg=%6.0f Rchg=%6.0f x=%6.2f y=%6.2f\n",
616  // cquad[quad], nhit, disc+1, p->fgtHitPhi->getLandauNorm(),p->fgtHitR->getLandauNorm(),p->x,p->y);
617  nhit++;
618  }
619  if(quad2 != quad) break;
620  for(int disc=0; disc<kFgtNumDiscs; disc++){
621  if(dischit[disc]==0){ //this track is in this quad, but no hit in this disc
622  double binFrac[1];
623  double z=StFgtGeom::getDiscZ(disc);
624  double x=t->mx*z+t->ax;
625  double y=t->my*z+t->ay;
626  double r=sqrt(x*x+y*y);
627  double p=atan2(y,x);
628  double localp=localphi(quad,p);
629  int ip = StFgtGeom::phi2LocalStripId(r,localp,binFrac);
630  int ir = StFgtGeom::rad2LocalStripId(r,localp,binFrac);
631  int gp= StFgtGeom::encodeGeoId(disc,quad,'P',ip); missedhit[gp]=1;
632  int gr= StFgtGeom::encodeGeoId(disc,quad,'R',ir); missedhit[gr]=1;
633  //printf("%1d%1s x=%6.2f y=%6.2f r=%6.2f ir=%3d phi=%6.2f local=%6.2f ip=%3d\n",
634  // disc+1,cquad[quad],x,y,r,ir,p,localp,ip);
635  }
636  }
637  }
638  //Print strips & clusters (with track if associated)
639  for (int disc=0; disc<kFgtNumDiscs; disc++){
640  StFgtStripCollection *cstrip = fgtCollectionPtr->getStripCollection(disc);
641  StSPtrVecFgtStrip &strips=cstrip->getStripVec();
642  StFgtHitCollection *chit = fgtCollectionPtr->getHitCollection(disc);
643  StSPtrVecFgtHit &hits=chit->getHitVec();
644  for(int ipr=0; ipr<2; ipr++){
645  for(int strip=0; strip<kFgtNumStrips; strip++){
646  int geoid = StFgtGeom::encodeGeoId(disc,quad,layer[ipr],strip);
647  if(geoid<0) continue;
648  StFgtStrip* str = 0;
649  StFgtHit* hit = 0;
650  AVTrack* trk = 0;
651  int flags=0, flagc=0, flagt=0;
652  for(StSPtrVecFgtStripIterator its=strips.begin(); its!=strips.end(); its++) {
653  if((*its)->getGeoId() != geoid) continue;
654  str = (*its);
655  if(str->getClusterSeedType()>0) flags=1;
656  }
657  for(StSPtrVecFgtHitIterator ith=hits.begin(); ith!=hits.end(); ith++) {
658  if((*ith)->getCentralStripGeoId() != geoid) continue;
659  hit = (*ith);
660  flagc=1;
661  }
662  for(vector<AVTrack>::iterator t=tracks.begin();t!=tracks.end();t++){
663  vector<AVPoint>* points=t->points;
664  for(vector<AVPoint>::iterator p=points->begin(); p!=points->end();p++){
665  if(p->quadID != quad) break;
666  if(p->dID != disc) continue;
667  StFgtHit* thit;
668  if(ipr==0) {thit=p->fgtHitPhi;}
669  else {thit=p->fgtHitR;}
670  if(thit->getCentralStripGeoId() != geoid) continue; //hit is in same disc/quad but somewhere else
671  flagt=1;
672  }
673  }
674  if(missedhit[geoid]>0) flagt=2;
675  if(flags>0 || flagc>0 || flagt>0){
676  int rdo, arm, apv, ch;
677  mDb->getElecCoordFromGeoId(geoid,rdo,arm,apv,ch);
678  Int_t elecId = StFgtGeom::encodeElectronicId( rdo, arm, apv, ch );
679  int status=mDb->getStatusFromElecId(elecId);
680  char cl;
681  short id,iq;
682  double rphi,l,h;
683  StFgtGeom::getGlobalPhysicalCoordinate(geoid,id,iq,cl,rphi,l,h);
684  fprintf(mTextFile,"%05d=%1d%1s%1c%03d|%05d=R%1dA%1dAPV%02dC%03d|",
685  geoid,disc+1,cquad[quad],layer[ipr],strip,elecId,rdo,arm,apv,ch);
686  if(str){
687  fprintf(mTextFile,"s=%1d psig=%3.0f|adc=", status,str->getPedErr());
688  for(int i=0; i<mNTimeBin; i++){ fprintf(mTextFile,"%4d ", str->getAdc(i)); }
689  fprintf(mTextFile," C=%7.0f sd=%02d",str->getCharge(), str->getClusterSeedType());
690  }
691  else {
692  fprintf(mTextFile,"s=%1d no strip info ", status);
693  for(int i=0; i<mNTimeBin; i++){ fprintf(mTextFile," ");}
694  }
695  fprintf(mTextFile," %1c=%8.4f ",layer[ipr],rphi);
696  if(flagc==1){
697  if(ipr==0){
698  fprintf(mTextFile,"C=%7.0f P=%8.3f N=%2d Sd=%2d",
699  hit->charge(),hit->getPositionPhi(),hit->getNstrip(),hit->getSeedType());
700  }else if(ipr==1){
701  fprintf(mTextFile,"C=%7.0f R=%8.3f N=%2d Sd=%2d",
702  hit->charge(),hit->getPositionR(),hit->getNstrip(),hit->getSeedType());
703  }
704  }
705  if(flagt==1) {fprintf(mTextFile," Track!");}
706  if(flagt==2) {fprintf(mTextFile," Miss (track should have left hit here)");}
707  fprintf(mTextFile,"\n");
708  }
709  }//strip
710  }//phir
711  }//disc
712  }//quad
713  NEVT++;
714 }
StFgtQAMaker::Finish
Int_t Finish()
Definition: StFgtQAMaker.cxx:98
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
kStOK
Definition: Stypes.h:40
StFgtQAMaker
A typical Analysis Class.
Definition: StFgtQAMaker.h:68
TObjectSet::GetObject
virtual TObject * GetObject() const
The depricated method (left here for the sake of the backward compatibility)
Definition: TObjectSet.h:56
Generated by   doxygen 1.8.5