StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFgtGenAVEMaker.cxx
1 #include "StFgtGenAVEMaker.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 
11 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
12 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
13 #include "StMuDSTMaker/COMMON/StMuDst.h"
14 #include "StMuDSTMaker/COMMON/StMuEvent.h"
15 #include "StarClassLibrary/StThreeVectorF.hh"
16 
17 
18 
19 #include <TH2D.h>
20 #include <TROOT.h>
21 #include <TStyle.h>
22 #include <TCanvas.h>
23 #include <utility>
24 #include <TArc.h>
25 #include <TLine.h>
26 #include <set>
27 
28 //#define PRINT_1D
29 //max num clusters any disk is allowed to have
30 
31 //#define COSMIC
32 #ifdef COSMIC
33 #include "StFgtCosmicAlignment.h"
34 #endif
35 #define MAX_PHI_DIFF 0.5//the maximal difference in phi a track is allowed
36 #define MAX_CLUSTERS 20
37 #define CHARGE_MEASURE clusterCharge
38 #define MAX_DIST_STRIP_R 0.7
39 #define MAX_DIST_STRIP_PHI 0.03
40 #include "StRoot/StFgtUtil/geometry/StFgtGeom.h"
41 //#define LEN_CONDITION
42 #define PULSE_CONDITION
43 #define DO_PRINT
44 #ifndef COSMIC
45 #define VERTEX_CUT 50
46 #else
47 #define VERTEX_CUT 100000000
48 #endif
49 
50 #include "StRoot/StEvent/StEvent.h"
51 #include "StRoot/StEvent/StFgtCollection.h"
52 //disk for which I want to calculate efficieny
53 
54 
55 #define MAX_CHARGE_RATIO
56 #define MIN_CHARGE_RATIO
57 
58 //#define DISK_EFF 2
59 #define QUAD_EFF 1
60 #define MY_PI 3.14159265359
61 //#define REFIT_WITH_VERTEX
62 //#define FIT_WITH_VERTEX
63 
64 
65 #define MAX_DIST_CHI 1.0
66 #define MAX_DIST2_EFF 1.0
67 #define MAX_DIST2 1.0
68 
69 #define MIN_NUM_POINTS 4
70 #define DISK_DIM 40
71 #ifndef COSMIC
72 #define NUM_EFF_BIN 30
73 #else
74 #define NUM_EFF_BIN 100
75 #endif
76 
77 Bool_t StFgtGenAVEMaker::fitTheStrip(generalStrip* pStrip, generalStrip* pStripOtherLayer,float* amp, float* t0, float* chi2Ndf, int iD, int iq, int apvBin, Char_t layer)
78 {
79  if(fitCounter>2000)
80  return true;
81  fitCounter++;
82  char buffer[100];
83  sprintf(buffer,"d%d_quad%d",iD,iq);
84  pulsePictureFile->cd();
85  pulsePictureFile->cd(buffer);
86  if(layer=='R')
87  sprintf(buffer,"apv%d_R",apvBin);
88  else
89  sprintf(buffer,"apv%d_P",apvBin);
90  // cout <<"changing to " << buffer <<endl;
91  gDirectory->cd(buffer);
92  Int_t minAdcCount=100000;
93  Int_t maxAdcCount=0;
94  for(Int_t tb=0;tb<7;tb++)
95  {
96  mHistPtr->SetBinContent(tb+1,0);
97  mHistPtr->SetBinError(tb+1,10000);
98  mHistPtr->SetBinContent(tb+1, pStrip->adc[tb]);
99  if(pStrip->adc[tb]<minAdcCount)
100  minAdcCount=pStrip->adc[tb]-pStrip->pedErr;
101  if(pStrip->adc[tb]>maxAdcCount)
102  maxAdcCount=pStrip->adc[tb]+pStrip->pedErr;
103  mHistPtr->SetBinError(tb+1,pStrip->pedErr);
104  }
105 
106  // mHistPtr->Fit(mPulseShapePtr);
107  (*amp)=mPulseShapePtr->GetParameter(0);
108  (*t0)=mPulseShapePtr->GetParameter(4);
109  (*chi2Ndf)=mPulseShapePtr->GetChisquare()/mPulseShapePtr->GetNDF();
110  sprintf(buffer,"pulse histo_D%d_Q%d_APV%d_ev%d",iD,iq,apvBin,evtNr);
111  // cout <<"histo name: "<< buffer <<endl;
112  TH1F* tmpPulseHisto=(TH1F*)mHistPtr->Clone(buffer);
113  tmpPulseHisto->Write();
114  if(pStripOtherLayer!=0)
115  {
116  mCanvas->cd();
117  sprintf(buffer,"tmpCnvsD%d_Q%d_APV%d_ev%d",iD,iq,apvBin,evtNr);
118  TCanvas* tmpCnv=(TCanvas*)mCanvas->Clone(buffer);
119  tmpCnv->SetTitle(buffer);
120  tmpCnv->SetName(buffer);
121  tmpCnv->cd();
122 
123  for(Int_t tb=0;tb<7;tb++)
124  {
125  mHistPtr2->SetBinContent(tb+1,0);
126  mHistPtr2->SetBinError(tb+1,10000);
127  mHistPtr2->SetBinContent(tb+1, pStripOtherLayer->adc[tb]);
128  mHistPtr2->SetBinError(tb+1,pStripOtherLayer->pedErr);
129  if(pStripOtherLayer->adc[tb]<minAdcCount)
130  minAdcCount=pStripOtherLayer->adc[tb]-pStripOtherLayer->pedErr;
131  if(pStripOtherLayer->adc[tb]>maxAdcCount)
132  maxAdcCount=pStripOtherLayer->adc[tb]+pStripOtherLayer->pedErr;
133  }
134  tmpPulseHisto->GetYaxis()->SetRangeUser(minAdcCount-10,maxAdcCount+10);
135  mHistPtr2->GetYaxis()->SetRangeUser(minAdcCount-10,maxAdcCount+10);
136  tmpPulseHisto->Draw();
137  mHistPtr2->Draw("SAME");
138  tmpCnv->Write();
139  }
140  pulsePictureFile->cd();
141  myRootFile->cd();
142  return true;
143 }
144 
145 
146 template<class T> void createPlots(T*** pH, int numH, const char* nameBase, int numBin, int first, int last)
147 {
148  char buffer[200];
149  (*pH)=new T*[numH];
150 
151  if (numH==22)
152  {
153  for(int iD=0;iD<numH;iD++)
154  {
155  sprintf(buffer, "%s_APV%d", nameBase, iD);
156  (*pH)[iD]=new T(buffer,buffer,numBin, first, last);
157  }
158  }
159 
160  if ((numH == kFgtNumDiscs*4)&&(numH!=22))
161  //if ((numH == kFgtNumDiscs*4))
162  {
163  for(int iD=0;iD<kFgtNumDiscs;iD++)
164  {
165  for(int iQ=0;iQ<4;iQ++)
166  {
167  sprintf(buffer,"%s_disc%d_quad%d",nameBase,iD+1,iQ);
168  (*pH)[iD*4+iQ]=new T(buffer,buffer,numBin,first,last);
169  }
170  }
171  }
172  else
173  {
174  for(int iD=0;iD<kFgtNumDiscs;iD++)
175  {
176 
177  if (numH==kFgtNumDiscs)
178  {
179  sprintf(buffer, "%s_disc%d", nameBase, iD+1);
180  (*pH)[iD]=new T(buffer, buffer,numBin, first, last);
181  }
182  else
183  {
184  for(int binAPVi=0;binAPVi<40;binAPVi++)
185  {
186  int iQ = -1;
187  if((binAPVi>= 0) && (binAPVi<= 9)) iQ=0;
188  if((binAPVi>=10) && (binAPVi<=19)) iQ=1;
189  if((binAPVi>=20) && (binAPVi<=29)) iQ=2;
190  if((binAPVi>=30) && (binAPVi<=39)) iQ=3;
191  sprintf(buffer,"%s_disc%d_quad%d_apvBIN%d",nameBase,iD+1,iQ,binAPVi);
192  (*pH)[iD*40+binAPVi]=new T(buffer,buffer,numBin,first,last);
193  }
194  }
195  }
196  }
197 }
198 
199 
200 template void createPlots(TH1I*** pH, int numH,const char* nameBase, int numBin, int first, int last);
201 template void createPlots(TH1F*** pH, int numH, const char* nameBase, int numBin, int first, int last);
202 
203 void StFgtGenAVEMaker::doNormalize(TH2D** hEff, TH2D** hNonEff)
204 {
205  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
206  {
207  TH2D* tmpAllCounts=(TH2D*)hEff[iD]->Clone("tmp");
208  hEff[iD]->Add(hNonEff[iD]);//all counts
209  for(int nx=1;nx<hEff[iD]->GetNbinsX()+1;nx++)
210  {
211  for(int ny=1;ny<hEff[iD]->GetNbinsY()+1;ny++)
212  {
213  Double_t denom=hEff[iD]->GetBinContent(nx,ny);
214  if(denom>0 && (tmpAllCounts->GetBinContent(nx,ny)/denom)<=1.0)
215  {
216  hEff[iD]->SetBinContent(nx,ny,tmpAllCounts->GetBinContent(nx,ny)/denom);
217  }
218  else
219  {
220  hEff[iD]->SetBinContent(nx,ny,0.0);
221  }
222  }
223  }
224  delete tmpAllCounts;
225  }
226 
227 }
228 
229 
230 
231 void StFgtGenAVEMaker::saveSigs(Double_t* sigR, Double_t* sigP, Double_t r, Double_t phi, Int_t maxR, Int_t maxPhi, Int_t discId, Int_t quad)
232 {
233  Char_t buffer[200];
234  sprintf(buffer,"Sig_Disc%d_quad%d_Phi_Evt_%d_R_%f_Phi_%f",discId,quad,evtNr,r,phi);
235  TH2D* histoP=new TH2D(buffer,buffer,7,0,6,maxPhi,0,maxPhi-1);
236  for(int i=0;i<maxPhi;i++)
237  {
238  for(int j=0;j<7;j++)
239  {
240  histoP->SetBinContent(j+1,i+1,sigP[i*7+j]);
241  // cout <<"P: setting bin : i: " << i << " j: " << j << " index: "<< i*7+j << " sig: " << sigP[i*7+j]<<endl;
242  }
243  }
244  v_hClusP.push_back(histoP);
245  sprintf(buffer,"Sig_Disc%d_quad%d_R_Evt_%d_R_%f_Phi_%f",discId,quad,evtNr,r,phi);
246  TH2D* histoR=new TH2D(buffer,buffer,7,0,6,maxR,0,maxR-1);
247  for(int i=0;i<maxR;i++)
248  {
249  for(int j=0;j<7;j++)
250  {
251  histoR->SetBinContent(j+1,i+1,sigR[i*7+j]);
252  // cout <<"R: setting bin : i: " << i << " j: " << j << " index: "<< i*7+j << " sig: " << sigR[i*7+j]<<endl;
253  }
254  }
255  v_hClusR.push_back(histoR);
256 }
257 
258 pair<double,double> StFgtGenAVEMaker::getDca( vector<AVTrack>::iterator it)
259 {
260 
261  float oldDist=100000;
262  float zStep=1.0;
263  float optZ=-200;
264  float dist=0;
265  for(float z=-100;z<100;z=z+zStep)
266  {
267  float x=it->mx*z+it->ax;
268  float y=it->my*z+it->ay;
269  dist=x*x+y*y;
270  if(dist>oldDist)
271  {
273  dist=oldDist;
274  optZ=z-zStep;//last z was better
275  break;
276  }
277  else
278  oldDist=dist;
279  }
280  return pair<double,double>(optZ,sqrt(dist));
281 }
282 
283 
284 Double_t StFgtGenAVEMaker::findClosestStrip(Char_t layer, double ord, Int_t iD, Int_t iQ)
285 {
286  if(iD<0 || iD >5)
287  {
288  return -99999;
289  }
290  vector<generalCluster> &hitVec=*(pClusters[iD]);
291  Double_t dist=99999;
292  for(vector<generalCluster>::iterator it=hitVec.begin();it!=hitVec.end();it++)
293  {
294  if(useChargeMatch && !it->hasMatch)
295  continue;
296  if(it->layer!=layer)
297  continue;
298  Double_t mDist=9999;
299  if(it->layer!='R')
300  {
301  mDist=fabs(it->posPhi-ord);
302  }
303  else
304  mDist=fabs(it->posR-ord);
305 
306  if(mDist<dist)
307  dist=mDist;
308  }
309  return dist;
310 }
311 
312 
315 void StFgtGenAVEMaker::fillStripHistos(Float_t r, Float_t phi, Int_t iD, Int_t iq)
316 {
317 
318  // cout <<"filling strip histos " <<endl;
319  bool partOfClusterP=false;
320  bool partOfClusterR=false;
321  Double_t clusterChargeR=-9999;
322  Double_t clusterChargeP=-9999;
323 
324  Double_t mClusterSizeP=-9999;
325  Double_t mClusterSizeR=-9999;
326 
327  Double_t maxRCharge=-9999;
328  Double_t maxPhiCharge=-9999;
329  Double_t maxRChargeUncert=-9999;
330  Double_t maxPhiChargeUncert=-9999;
331  Int_t maxRInd=-1;
332  Int_t maxPInd=-1;
333  Int_t maxRTb=-1;
334  Int_t maxPTb=-1;
335  Int_t maxRAdc=-1;
336  Int_t maxPAdc=-1;
337  Double_t maxSigAdcR=-1;
338  Double_t maxSigAdcP=-1;
339  Int_t numFSigP=0;
340  Int_t numFSigR=0;
341  Int_t firstFSigTbP=0;
342  Int_t firstFSigTbR=0;
343 
344  Float_t secondToLastRatioP=0;
345  Float_t secondToLastRatioR=0;
346 
347  Float_t firstTbSigR=-1;
348  Float_t firstTbSigP=-1;
349 
350 
351  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
352  Short_t disc, quadrant,strip;
353  Char_t layer;
354  //These variables will store the APV number at which the relevant stuff happens for the overall quad.
355  Int_t APVmaxRCharge=-9999;
356  Int_t APVmaxPhiCharge=-9999;
357  //Int_t APVmaxRChargeUncert=-9999;
358  //Int_t APVmaxPhiChargeUncert=-9999;
359  Int_t APVmaxRInd=-1;
360  Int_t APVmaxPInd=-1;
361  Int_t APVmaxRTb=-1;
362  Int_t APVmaxPTb=-1;
363  Int_t APVmaxRAdc=-1;
364  Int_t APVmaxPAdc=-1;
365  Int_t APVmaxSigAdcR=-1;
366  Int_t APVmaxSigAdcP=-1;
367  Int_t APVnumFSigP=0;
368  Int_t APVnumFSigR=0;
369  Int_t APVfirstFSigTbP=0;
370  Int_t APVfirstFSigTbR=0;
371 
372  Int_t APVsecondToLastRatioP=0;
373  Int_t APVsecondToLastRatioR=0;
374 
375  Int_t APVfirstTbSigR=-1;
376  Int_t APVfirstTbSigP=-1;
377 
378 
379 
380  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
381  {
382  Int_t geoId=pStrips[iD*4+iq][i].geoId;
383  generalStrip& pStrip=pStrips[iD*4+iq][i];
384 
385  Int_t rdo, arm, apv, chan;
386  mDb->getElecCoordFromGeoId(geoId, rdo,arm,apv,chan);
387  //using binAPV for simplicity
388  Int_t binAPV = (iq*10)+(apv%12);
389 
390  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
391  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
392  //char buffer[100];
393  // if(layer=='P' && disc==iD && iq==quadrant)
394  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
395  if(disc==iD && iq==quadrant && ((layer =='R' && fabs(ordinate-r)<0.7) || (layer=='P' && fabs(ordinate-phi)<0.03) || (layer=='P' && fabs(ordinate-phi+2*MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-2*MY_PI)<0.03)|| (layer=='P' && fabs(ordinate-phi+MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-MY_PI)<0.03)))
396  {
397  if(layer=='P')
398  {
399  if(pStrip.charge>maxPhiCharge)
400  {
401  if(pStrip.seedType==kFgtSeedType1|| pStrip.seedType==kFgtSeedType2 || pStrip.seedType==kFgtSeedType3 || pStrip.seedType==kFgtClusterPart || pStrip.seedType==kFgtClusterEndUp ||pStrip.seedType==kFgtClusterEndDown)
402  {
403  partOfClusterP=true;
404  pair<Double_t, Double_t> cluSize=findCluChargeSize(iD,'P',ordinate);
405  clusterChargeP=cluSize.first;
406  mClusterSizeP=cluSize.second;
407  }
408  else
409  {
410  partOfClusterP=false;
411  clusterChargeP=-9999;
412  }
413  maxPhiCharge=pStrip.charge;
414  maxPhiChargeUncert=pStrip.chargeUncert;
415  maxPInd=i;
416  maxPAdc=-9999;
417  maxSigAdcP=-9999;
418  numFSigP=0;
419  firstFSigTbP=-1;
420  firstTbSigP=-1;
421  APVmaxPhiCharge=-1;
422  APVmaxPInd=-1;
423  APVmaxPAdc=-1;
424  APVmaxSigAdcP=-1;
425  APVnumFSigP=-1;
426  APVfirstFSigTbP=-1;
427  APVfirstTbSigP=-1;
428 
429  if(pStrip.adc[6]>0)
430  {
431  secondToLastRatioP=pStrip.adc[5]/(float)pStrip.adc[6];
432  APVsecondToLastRatioP=binAPV;
433  }
434  if(pStrip.pedErr>0)
435  {
436  firstTbSigP=pStrip.adc[0]/(float)pStrip.pedErr;
437  APVfirstTbSigP=binAPV;
438  }
439  for(int iAdc=0;iAdc<7;iAdc++)
440  {
441  if(pStrip.adc[iAdc]>5*pStrip.pedErr)
442  {
443  numFSigP++;
444  APVnumFSigP=binAPV;
445  if(firstFSigTbP<0)
446  {
447  firstFSigTbP=iAdc;
448  APVfirstFSigTbP=binAPV;
449  }
450  }
451  if(pStrip.adc[iAdc]>maxPAdc)
452  {
453  maxPAdc=pStrip.adc[iAdc];
454  APVmaxPAdc=binAPV;
455  maxPTb=iAdc;
456  APVmaxPTb=binAPV;
457  if(pStrip.pedErr>0)
458  {
459  maxSigAdcP=(Double_t)maxPAdc/pStrip.pedErr;
460  APVmaxSigAdcP=binAPV;
461  }
462  }
463  }
464  }
465  }
466  else
467  {
468  if(pStrip.charge>maxRCharge)
469  {
470  if(pStrip.seedType==kFgtSeedType1|| pStrip.seedType==kFgtSeedType2 || pStrip.seedType==kFgtSeedType3 || pStrip.seedType==kFgtClusterPart || pStrip.seedType==kFgtClusterEndUp ||pStrip.seedType==kFgtClusterEndDown)
471  {
472  partOfClusterR=true;
473  pair<Double_t, Double_t> cluSize=findCluChargeSize(iD,'R',ordinate);
474  clusterChargeR=cluSize.first;
475  mClusterSizeR=cluSize.second;
476  }
477  else
478  {
479  clusterChargeR=-9999;
480  partOfClusterR=false;
481  }
482 
483  maxRCharge=pStrip.charge;
484  maxRInd=i;
485  maxRChargeUncert=pStrip.chargeUncert;
486  maxRAdc=-9999;
487  maxSigAdcR=-9999;
488  numFSigR=0;
489  firstFSigTbR=-1;
490  firstTbSigR=-1;
491  APVmaxRCharge=-1;
492  APVmaxRInd=-1;
493  APVmaxRAdc=-1;
494  APVmaxSigAdcR=-1;
495  APVnumFSigR=-1;
496  APVfirstFSigTbR=-1;
497  APVfirstTbSigR=-1;
498 
499  if(pStrip.adc[6]>0)
500  {
501  secondToLastRatioR=pStrip.adc[5]/(float)pStrip.adc[6];
502  APVsecondToLastRatioR=binAPV;
503  }
504  if(pStrip.pedErr>0)
505  {
506  firstTbSigR=pStrip.adc[0]/(float)pStrip.pedErr;
507  APVfirstTbSigR=binAPV;
508  }
509  for(int iAdc=0;iAdc<7;iAdc++)
510  {
511  // cout <<"adc: "<< pStrip.adc[iAdc] <<endl;
512  if(pStrip.adc[iAdc]>5*pStrip.pedErr)
513  {
514  numFSigR++;
515  APVnumFSigR=binAPV;
516  if(firstFSigTbR<0)
517  {
518  firstFSigTbR=iAdc;
519  APVfirstFSigTbR=binAPV;
520  }
521  }
522  if(pStrip.adc[iAdc]>maxRAdc)
523  {
524  maxRAdc=pStrip.adc[iAdc];
525  APVmaxRAdc=binAPV;
526  maxRTb=iAdc;
527  APVmaxRTb=binAPV;
528  if(pStrip.pedErr>0)
529  {
530  maxSigAdcR=(Double_t)maxRAdc/pStrip.pedErr;
531  APVmaxSigAdcR=binAPV;
532  }
533  }
534  }
535  // cout <<"numfSigma: " << numFSigR <<endl;
536  }
537  }
538  }
539  }
540  //fill histos
541  if(maxPhiCharge>200)
542  {
543  firstTbSigCloseClusterP[iD*4+iq]->Fill(firstTbSigP);
544  maxAdcCloseClusterP[iD*4+iq]->Fill(maxPAdc);
545 
546  maxTbCloseClusterP[iD*4+iq]->Fill(maxPTb);
547  numFSigCloseClusterP[iD*4+iq]->Fill(numFSigP);
548  numFirstHighCloseClusterP[iD*4+iq]->Fill(firstFSigTbP);
549  // maxAdcCloseClusterP[iD*4+iq]->Fill(maxPAdc);
550  maxSigCloseClusterP[iD*4+iq]->Fill(maxSigAdcP);
551  // cout <<"firstFSigtbP: " << firstFSigTbP <<endl;
552  // cout <<"filling ratio with : " << secondToLastRatioP <<endl;
553  secondToLastRatioCloseClusterP[iD*4+iq]->Fill(secondToLastRatioP);
554 
555  if(APVfirstTbSigP>-1 && APVmaxPAdc>-1 && APVmaxPTb>-1 && APVnumFSigP>-1 && APVfirstFSigTbP>-1 && APVmaxSigAdcP>-1 && APVsecondToLastRatioP>-1)
556  {
557  APVfirstTbSigCloseClusterP[iD*40+APVfirstTbSigP]->Fill(firstTbSigP);
558  APVmaxAdcCloseClusterP[iD*40+APVmaxPAdc]->Fill(maxPAdc);
559  APVmaxTbCloseClusterP[iD*40+APVmaxPTb]->Fill(maxPTb);
560  APVnumFSigCloseClusterP[iD*40+APVnumFSigP]->Fill(numFSigP);
561  APVnumFirstHighCloseClusterP[iD*40+APVfirstFSigTbP]->Fill(firstFSigTbP);
562  APVmaxSigCloseClusterP[iD*40+APVmaxSigAdcP]->Fill(maxSigAdcP);
563  APVsecondToLastRatioCloseClusterP[iD*40+APVsecondToLastRatioP]->Fill(secondToLastRatioP);
564  }
565 
566  if(iD==2 && iq==1 && (float)pStrips[iD*4+iq][maxPInd].pedErr>0)
567  {
568  pulseCounterP++;
569  for(int iB=0;iB<7;iB++)
570  {
571  exPulseMaxAdcNormP->SetBinContent(iB+1,exPulseMaxAdcNormP->GetBinContent(iB+1)+pStrips[iD*4+iq][maxPInd].adc[iB]/(float)maxPAdc);
572  exPulseSigP->SetBinContent(iB+1,exPulseSigP->GetBinContent(iB+1)+pStrips[iD*4+iq][maxPInd].adc[iB]/(float)pStrips[iD*4+iq][maxPInd].pedErr);
573  }
574  }
575  }
576 
577  if(maxRCharge>200)
578  {
579  firstTbSigCloseClusterR[iD*4+iq]->Fill(firstTbSigR);
580  maxTbCloseClusterR[iD*4+iq]->Fill(maxRTb);
581  maxAdcCloseClusterR[iD*4+iq]->Fill(maxRAdc);
582  numFSigCloseClusterR[iD*4+iq]->Fill(numFSigR);
583  numFirstHighCloseClusterR[iD*4+iq]->Fill(firstFSigTbR);
584  // cout <<"firstFSigtbR: " << firstFSigTbR <<endl;
585  maxSigCloseClusterR[iD*4+iq]->Fill(maxSigAdcR);
586  // cout <<"filling ratio R with : " << secondToLastRatioR <<endl;
587  secondToLastRatioCloseClusterR[iD*4+iq]->Fill(secondToLastRatioR);
588 
589  if(APVfirstTbSigR>-1 && APVmaxRAdc>-1 && APVmaxRTb>-1 && APVnumFSigR>-1 && APVfirstFSigTbR>-1 && APVmaxSigAdcR>-1 && APVsecondToLastRatioR>-1)
590  {
591  APVfirstTbSigCloseClusterR[iD*40+APVfirstTbSigR]->Fill(firstTbSigR);
592  APVmaxTbCloseClusterR[iD*40+APVmaxRTb]->Fill(maxRTb);
593  APVmaxAdcCloseClusterR[iD*40+APVmaxRAdc]->Fill(maxRAdc);
594  APVnumFSigCloseClusterR[iD*40+APVnumFSigR]->Fill(numFSigR);
595  APVnumFirstHighCloseClusterR[iD*40+APVfirstFSigTbR]->Fill(firstFSigTbR);
596  APVmaxSigCloseClusterR[iD*40+APVmaxSigAdcR]->Fill(maxSigAdcR);
597  APVsecondToLastRatioCloseClusterR[iD*40+APVsecondToLastRatioR]->Fill(secondToLastRatioR);
598  }
599 
600  if(iD==2 && iq==1 && (float)pStrips[iD*4+iq][maxRInd].pedErr>0)
601  {
602  pulseCounterR++;
603  for(int iB=0;iB<7;iB++)
604  {
605  exPulseMaxAdcNormR->SetBinContent(iB+1,exPulseMaxAdcNormR->GetBinContent(iB+1)+pStrips[iD*4+iq][maxRInd].adc[iB]/(float)maxRAdc);
606  exPulseSigR->SetBinContent(iB+1,exPulseSigR->GetBinContent(iB+1)+pStrips[iD*4+iq][maxRInd].adc[iB]/(float)pStrips[iD*4+iq][maxRInd].pedErr);
607  }
608  }
609  }
610 
611 
612  if(maxRCharge> 200 && maxPhiCharge>200 && iD==m_effDisk)// && (float)pStrips[iD*4+iq][maxRInd].charge)
613  {
614  StFgtGeom::getPhysicalCoordinate((float)pStrips[iD*4+iq][maxRInd].geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
615  // if(ordinate>20)
616  {
617  chargeCorrMaxStrip->Fill(maxRCharge,maxPhiCharge);
618  chargeCorrMaxAdc->Fill(maxRAdc,maxPAdc);
619 
620 
621  }
622  }
623  // if(iD==2 && iq==1 && (float)pStrips[iD*4+iq][maxRInd].pedErr>0)
624 
625 
626 
627  if(partOfClusterR&& partOfClusterP && iD==m_effDisk)
628  {
629  chargeCorrInEffDisk->Fill(clusterChargeR,clusterChargeP);
630  StFgtGeom::getPhysicalCoordinate((float)pStrips[iD*4+iq][maxRInd].geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
631  // if(ordinate>20)
632  {
633  chargeCorr[iD*4+iq]->Fill(clusterChargeR,clusterChargeP);
634  }
635  //basically only here we fill with the estimated cluster, for the other disks we fill with the cluster on the track
636  // if(r>20)
637  {
638  clusterSizeR[iD*4+iq]->Fill(mClusterSizeR);
639  clusterSizeP[iD*4+iq]->Fill(mClusterSizeP);
640  }
641  }
642 
643 
644  if(partOfClusterP)
645  {
646 
647  // cout <<" part of p cluster " << endl;
648  float intPCharge=(float)pStrips[iD*4+iq][maxPInd].charge+(float)pStrips[iD*4+iq][maxPInd-1].charge+(float)pStrips[iD*4+iq][maxPInd+1].charge;
649  float chi2Ndf;
650  float amp;
651  float t0;
653  if(partOfClusterR)
654  fitTheStrip(&(pStrips[iD*4+iq][maxPInd]),&(pStrips[iD*4+iq][maxRInd]),&amp,&t0,&chi2Ndf,iD,iq,APVmaxPAdc,'P');
655  else
656  fitTheStrip(&(pStrips[iD*4+iq][maxPInd]),0,&amp,&t0,&chi2Ndf,iD,iq,APVmaxPAdc,'P');
657 
658  APVfitChi2P[iD*40+APVmaxPAdc]->Fill(chi2Ndf);
659 
660 
661  firstTbSigTrackClusterP[iD*4+iq]->Fill(firstTbSigP);
662  maxAdcTrackClusterP[iD*4+iq]->Fill(maxPAdc);
663  maxSigTrackClusterP[iD*4+iq]->Fill(maxSigAdcP);
664  maxTbTrackClusterP[iD*4+iq]->Fill(maxPTb);
665  numFSigTrackClusterP[iD*4+iq]->Fill(numFSigP);
666  numFirstHighTrackClusterP[iD*4+iq]->Fill(firstFSigTbP);
667  secondToLastRatioTrackClusterP[iD*4+iq]->Fill(secondToLastRatioP);
668 
669  if(iD==2 && iq==1 && (float)pStrips[iD*4+iq][maxPInd].pedErr>0)
670  {
671  pulseCounterTP++;
672  for(int iB=0;iB<7;iB++)
673  {
674  exPulseMaxAdcNormTrackP->SetBinContent(iB+1,exPulseMaxAdcNormTrackP->GetBinContent(iB+1)+pStrips[iD*4+iq][maxPInd].adc[iB]/(float)maxPAdc);
675 
676  exPulseSigTrackP->SetBinContent(iB+1,exPulseSigTrackP->GetBinContent(iB+1)+pStrips[iD*4+iq][maxPInd].adc[iB]/(float)pStrips[iD*4+iq][maxPInd].pedErr);
677 
678  }
679  }
680 
681  }
682  if(partOfClusterR)
683  {
684 
685  // cout <<" part of R cluster " << endl;
686  float intRCharge=(float)pStrips[iD*4+iq][maxRInd].charge+(float)pStrips[iD*4+iq][maxRInd-1].charge+(float)pStrips[iD*4+iq][maxRInd+1].charge;
687  float chi2Ndf;
688  float amp;
689  float t0;
691  if(partOfClusterP)
692  fitTheStrip(&(pStrips[iD*4+iq][maxRInd]),&pStrips[iD*4+iq][maxPInd],&amp,&t0,&chi2Ndf, iD, iq, APVmaxRAdc,'R');
693  else
694  fitTheStrip(&(pStrips[iD*4+iq][maxRInd]),0,&amp,&t0,&chi2Ndf, iD, iq, APVmaxRAdc,'R');
695  APVfitChi2R[iD*40+APVmaxRAdc]->Fill(chi2Ndf);
696 
697 
698  firstTbSigTrackClusterR[iD*4+iq]->Fill(firstTbSigR);
699  maxTbTrackClusterR[iD*4+iq]->Fill(maxRTb);
700  maxAdcTrackClusterR[iD*4+iq]->Fill(maxRAdc);
701  maxSigTrackClusterR[iD*4+iq]->Fill(maxSigAdcR);
702  numFSigTrackClusterR[iD*4+iq]->Fill(numFSigR);
703  numFirstHighTrackClusterR[iD*4+iq]->Fill(firstFSigTbR);
704  secondToLastRatioTrackClusterR[iD*4+iq]->Fill(secondToLastRatioR);
705 
706  if(iD==2 && iq==1 && (float)pStrips[iD*4+iq][maxRInd].pedErr>0)
707  {
708  pulseCounterTR++;
709  for(int iB=0;iB<7;iB++)
710  {
711  exPulseMaxAdcNormTrackR->SetBinContent(iB+1,exPulseMaxAdcNormTrackR->GetBinContent(iB+1)+pStrips[iD*4+iq][maxRInd].adc[iB]/(float)maxRAdc);
712  exPulseSigTrackR->SetBinContent(iB+1,exPulseSigTrackR->GetBinContent(iB+1)+pStrips[iD*4+iq][maxRInd].adc[iB]/(float)pStrips[iD*4+iq][maxRInd].pedErr);
713  }
714  }
715  }
716 
717 }
718 //is this a good hit in the loose term. We can implement different criteria, so either pulse or some charge sum
719 //Bool_t StFgtGenAVEMaker::isGoodHit(
720 
721 
722 //is there something where we expect it? Again, here phi is relative to the quadrant!!
723 Bool_t StFgtGenAVEMaker::isSomewhatEff(Float_t r, Float_t phi, Int_t iD, Int_t iq)
724 {
725  Double_t maxRCharge=-9999;
726  Double_t maxPhiCharge=-9999;
727  Double_t maxRChargeUncert=-9999;
728  Double_t maxPhiChargeUncert=-9999;
729  Int_t maxRInd=-1;
730  Int_t maxPInd=-1;
731  Bool_t validPhiPulse=false;
732  Bool_t validRPulse=false;
733  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
734  {
735  Int_t geoId=pStrips[iD*4+iq][i].geoId;
736  generalStrip& pStrip=pStrips[iD*4+iq][i];
737  Short_t disc, quadrant,strip;
738  Char_t layer;
739  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
740  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
741  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
742  // if(layer=='P' && disc==iD && iq==quadrant)
743  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
744  if(disc==iD && iq==quadrant && ((layer =='R' && fabs(ordinate-r)<0.7) || (layer=='P' && fabs(ordinate-phi)<0.04) || (layer=='P' && fabs(ordinate-phi+2*MY_PI)<0.04 ) || (layer=='P' && fabs(ordinate-phi-2*MY_PI)<0.04)|| (layer=='P' && fabs(ordinate-phi+MY_PI)<0.04 ) || (layer=='P' && fabs(ordinate-phi-MY_PI)<0.04)))
745  {
746  if(layer=='P')
747  {
748  if(validPulse(pStrip))
749  validPhiPulse=true;
750 
751  if(pStrip.charge>maxPhiCharge)
752  {
753  maxPhiCharge=pStrip.charge;
754  maxPhiChargeUncert=pStrip.chargeUncert;
755  maxPInd=i;
756  }
757  }
758  else
759  {
760  if(validPulse(pStrip))
761  validRPulse=true;
762 
763  if(pStrip.charge>maxRCharge)
764  {
765  maxRCharge=pStrip.charge;
766  maxRInd=i;
767  maxRChargeUncert=pStrip.chargeUncert;
768  }
769  }
770  }
771  }
772 #ifdef PULSE_CONDITION
773  if(validPhiPulse && validRPulse)
774  return true;
775 #endif
776 
777 
778  if(maxRInd>=0 && maxPInd>=0)
779  {
780 #ifdef LEN_CONDITION
781  if(maxRCharge>1000 && maxPhiCharge>1000)
782  return true;
783 #endif
784 
786  if(maxRInd>0)
787  maxRCharge+=pStrips[iD*4+iq][maxRInd-1].charge;
788  if(maxRInd< (int)(pStrips[iD*4+iq].size()-1))
789  maxRCharge+=pStrips[iD*4+iq][maxRInd+1].charge;
790  if(maxPInd>0)
791  maxPhiCharge+=pStrips[iD*4+iq][maxPInd-1].charge;
792  if(maxPInd< (int)(pStrips[iD*4+iq].size()-1))
793  maxPhiCharge+=pStrips[iD*4+iq][maxPInd+1].charge;
795 
796  // cout <<"charge: "<< maxRCharge<< " 3* uncert " << 3*maxRChargeUncert <<" maxPhiCharge " << maxPhiCharge<<" 3*phi "<< 3*maxPhiChargeUncert <<endl;
797  if(maxRCharge>3*maxRChargeUncert && maxPhiCharge>3*maxPhiChargeUncert)
798  {
799  //if we don't care about charge ratio
800  // return true;
801  if(maxRCharge>maxPhiCharge)
802  {
803  if(maxRCharge/maxPhiCharge<4)
804  return true;
805  else
806  return false;
807  }
808  else
809  {
810  if(maxPhiCharge/maxPhiCharge<4)
811  return true;
812  else
813  return false;
814  }
815  }
816 
817  }
818  return false;
819 }
820 
821 //if required only return clusters with matches...
822 pair<Double_t,Double_t> StFgtGenAVEMaker::findCluChargeSize(Int_t iD,Char_t layer, Double_t ordinate)
823 {
824  if(iD<0 || iD >5)
825  {
826  return pair<Double_t,Double_t>(-99999,-99999);
827  }
828  vector<generalCluster> &hitVec=*(pClusters[iD]);
829  Double_t charge=-99999;
830  Double_t cluSize=-9999;
831  Double_t minDist=99999;
832  for(vector<generalCluster>::iterator it=hitVec.begin();it!=hitVec.end();it++)
833  {
834  if(it->layer!=layer)
835  continue;
836  Float_t ord=-9999;
837  if(layer=='R')
838  {
839  ord=it->posR;
840  }
841  else
842  {
843  ord=it->posPhi;
844  if(fabs(ordinate-(ord+MY_PI))<fabs(ord-ordinate))
845  ord=ord+MY_PI;
846  if(fabs(ordinate-(ord-MY_PI))<fabs(ord-ordinate))
847  ord=ord-MY_PI;
848  }
849 
850  if((fabs(ordinate-ord)<minDist) && (!useChargeMatch || it->hasMatch))
851  {
852  charge=it->clusterCharge;
853  cluSize=it->clusterSize;
854  minDist=fabs(ordinate-ord);
855  }
856  }
857  return pair<Double_t,Double_t>(charge,cluSize);
858 }
859 
860 Double_t StFgtGenAVEMaker::findClosestPoint(float mx, float bx, float my, float by, double xE, double yE, Int_t iD)
861 {
862  // cout <<"expecting point at " << xE <<", " <<yE <<endl;
863  if(iD<0 || iD >5)
864  {
865  return 99999;
866  }
867  vector<generalCluster> &hitVec=*(pClusters[iD]);
868  Double_t dist2=99999;
869  for(vector<generalCluster>::iterator it=hitVec.begin();it!=hitVec.end();it++)
870  {
871  for(vector<generalCluster>::iterator it2=hitVec.begin();it2!=hitVec.end();it2++)
872  {
873  //
874  if(useChargeMatch && !arePointsMatched(it,it2))
875  continue;
876 
877  if(it->layer==it2->layer)
878  continue;
879  Float_t r=it->posR;
880  Float_t phi=it2->posPhi;
881  if(it->layer!='R')
882  {
883  phi=it->posPhi;
884  r=it2->posR;
885  }
886 
887  Float_t x=r*cos(phi);
888  Float_t y=r*sin(phi);
890  Double_t mDist=(x-xE)*(x-xE)+(y-yE)*(y-yE);
891 
892  if(mDist<dist2)
893  {
894  dist2=mDist;
895  //recalculate distance with proper alignment
896  float tmpX, tmpY,tmpZ,tmpP,tmpR;
897 #ifdef COSMIC
898  getAlign(iD,phi,r,tmpX,tmpY,tmpZ,tmpP,tmpR);
899  Double_t xExpUpdate=mx*tmpZ+bx;
900  Double_t yExpUpdate=my*tmpZ+by;
901  // cout<<"tmpx: " << tmpX <<" old: " << x <<" xE old: " << xE << " updated: " << xExpUpdate;
902  // cout<<"tmpy: " << tmpY <<" old: " << y <<" yE old: " << yE << " updated: " << yExpUpdate<<endl;
903 
904  mDist=(tmpX-xExpUpdate)*(tmpX-xExpUpdate)+(tmpY-yExpUpdate)*(tmpY-yExpUpdate);
905 #endif
906  dist2=mDist;
908  // Double_t yExp=my*StFgtGeom::getDiscZ(i)+by;
909 
910 
911  // (*outTxtFile) <<"point found, x: " << x <<" y: " << y << " dist: " << dist2 <<endl;
912  }
913  }
914  }
915  // (*outTxtFile) <<"returning : " << dist2<<endl;
916  return dist2;
917 }
918 
919 
921 Short_t StFgtGenAVEMaker::getQuadFromCoo(Double_t x, Double_t y)
922 {
923  cout <<"do not use this function!!!" <<endl;
924  if(x>0 && y>0)
925  return 0;
926  if(x>0 && y<0)
927  return 1;
928  if(x<0 && y<0)
929  return 2;
930  if(x<0 && y>0)
931  return 3;
932 
933  return -9999;
934 }
935 
936 
937 Bool_t StFgtGenAVEMaker::printArea1D(Int_t iD,Int_t iq, Int_t centerGeoId)
938 {
939 
940  for(unsigned int i=0;i<pStrips[iD*4+iq].size();i++)
941  {
942  Int_t geoId=pStrips[iD*4+iq][i].geoId;
943  generalStrip& pStrip=pStrips[iD*4+iq][i];
944  Short_t disc, quadrant,strip;
945  Char_t layer;
946  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
947  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
948  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
949  char buffer[100];
950  switch(pStrip.seedType)
951  {
952  case kFgtSeedTypeNo:
953  sprintf(buffer,"No Seed");
954  break;
955  case kFgtSeedType1:
956  sprintf(buffer,"Seed1");
957  break;
958  case kFgtSeedType2:
959  sprintf(buffer,"Seed2");
960  break;
961  case kFgtSeedType3:
962  sprintf(buffer,"Seed3");
963  break;
964  case kFgtClusterPart:
965  sprintf(buffer,"PartOfCluster");
966  break;
967  case kFgtClusterEndUp:
968  sprintf(buffer,"EoC");
969  break;
970  case kFgtClusterEndDown:
971  sprintf(buffer,"BoC");
972  break;
973  case kFgtDeadStrip:
974  sprintf(buffer,"DeadStrip");
975  break;
976  case kFgtClusterTooBig:
977  sprintf(buffer,"cluster too big");
978  break;
979  case kFgtClusterSeedInSeaOfNoise:
980  sprintf(buffer,"seed in noise");
981  break;
982  default:
983  sprintf(buffer,"somethingWrong: %d", pStrip.seedType);
984  }
985 
986  // if(layer=='P' && disc==iD && iq==quadrant)
987  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
988  if(abs(geoId-centerGeoId)<8)
989  {
990  // cout <<" found!!!" << endl;
991  (*outTxtFile) <<StFgtGeom::encodeGeoName(iD,iq,layer,strip)<<"geo: " << geoId<< " ord: " << ordinate <<" layer: " <<layer<<" ped: " << pStrip.ped <<" pedErr: " << pStrip.pedErr <<" seedType: " <<buffer<<" ";
992  for(int iT=0;iT<7;iT++)
993  {
994  if(pStrip.adc[iT]<pStrip.pedErr)
995  (*outTxtFile) << setw(4) << " . "<< " ";
996  else
997  (*outTxtFile) << setw(4) <<pStrip.adc[iT] <<" ";
998  }
999  (*outTxtFile) <<endl;
1000  }
1001  }
1002  return kStOk;
1003 
1004 }
1005 
1006 //print the strips around the place where we expect hit
1007 Bool_t StFgtGenAVEMaker::printArea(Float_t r, Float_t phi, Int_t iD, Int_t iq)
1008 {
1009  //just print the first 1000 clusters
1010 
1011  // cout <<" looking for r: " << r <<" phi: " << phi<< " to print " <<endl;
1012  if(printCounter>1000)
1013  return true;
1014  printCounter++;
1015  //first r:
1016  Double_t signalsP[900];
1017  Double_t signalsR[900];
1018 
1019  Int_t counterR=0;
1020  Int_t counterP=0;
1021  // cout <<" looking for strips.., we have: " << pStrips[iD*4+iq].size() << " in iD: " << iD << " iq: " << iq <<endl;
1022  for(unsigned int i=0;i<pStrips[iD*4+iq].size();i++)
1023  {
1024  Int_t geoId=pStrips[iD*4+iq][i].geoId;
1025  generalStrip& pStrip=pStrips[iD*4+iq][i];
1026  Short_t disc, quadrant,strip;
1027  Char_t layer;
1028  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
1029  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
1030  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
1031  char buffer[100];
1032  switch(pStrip.seedType)
1033  {
1034  case kFgtSeedTypeNo:
1035  sprintf(buffer,"No Seed");
1036  break;
1037  case kFgtSeedType1:
1038  sprintf(buffer,"Seed1");
1039  break;
1040  case kFgtSeedType2:
1041  sprintf(buffer,"Seed2");
1042  break;
1043  case kFgtSeedType3:
1044  sprintf(buffer,"Seed3");
1045  break;
1046  case kFgtClusterPart:
1047  sprintf(buffer,"PartOfCluster");
1048  break;
1049  case kFgtClusterEndUp:
1050  sprintf(buffer,"EoC");
1051  break;
1052  case kFgtClusterEndDown:
1053  sprintf(buffer,"BoC");
1054  break;
1055  case kFgtDeadStrip:
1056  sprintf(buffer,"DeadStrip");
1057  break;
1058  case kFgtClusterTooBig:
1059  sprintf(buffer,"cluster too big");
1060  break;
1061  case kFgtClusterSeedInSeaOfNoise:
1062  sprintf(buffer,"seed in noise");
1063  break;
1064  default:
1065  sprintf(buffer,"somethingWrong: %d", pStrip.seedType);
1066  }
1067 
1068  // if(layer=='P' && disc==iD && iq==quadrant)
1069  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
1070  // if(layer=='R' && disc==iD && iq==quadrant)
1071  // cout <<"looking for r=" << r << " have: " << ordinate <<" diff: " << fabs(ordinate-r) <<endl;
1072  if(disc==iD && iq==quadrant && ((layer =='R' && fabs(ordinate-r)<1.0) || (layer=='P' && fabs(ordinate-phi)<0.04) || (layer=='P' && fabs(ordinate-phi+2*MY_PI)<0.04 ) || (layer=='P' && fabs(ordinate-phi-2*MY_PI)<0.04)|| (layer=='P' && fabs(ordinate-phi+MY_PI)<0.04 ) || (layer=='P' && fabs(ordinate-phi-MY_PI)<0.04)))
1073  {
1074  // cout <<" found!!!" << endl;
1075  (*outTxtFile) <<"Id: " << geoId<<" "<<StFgtGeom::encodeGeoName(iD,iq,layer,strip)<<" ord: " << ordinate <<" layer: " <<layer<<" ped: " << pStrip.ped <<" pedErr: " << pStrip.pedErr <<" seedType: " <<buffer<<" ";
1076  for(int iT=0;iT<7;iT++)
1077  {
1078  if(pStrip.adc[iT]<pStrip.pedErr)
1079  (*outTxtFile) << setw(4) << " . "<< " ";
1080  else
1081  (*outTxtFile) << setw(4) <<pStrip.adc[iT] <<" ";
1082 
1083 
1084  if(layer=='P'&& counterP<40)
1085  {
1086  signalsP[counterP*7+iT]=pStrip.adc[iT];
1087  // cout <<"first p: setting bin : i: " << counterP << " j: " << iT << " index: "<< counterP*7+iT << " sig: " << signalsP[counterP*7+iT]<<endl;
1088  if(iT==6)
1089  counterP++;
1090  }
1091  if(layer=='R'&& counterR<40)
1092  {
1093  signalsR[counterR*7+iT]=pStrip.adc[iT];
1094  // cout <<"first R: setting bin : i: " << counterR << " j: " << iT << " index: "<< counterR*7+iT << " sig: " << signalsR[counterP*7+iT]<<endl;
1095  if(iT==6)
1096  counterR++;
1097  }
1098  }
1099  (*outTxtFile) <<endl;
1100  }
1101  }
1102 
1103  saveSigs(signalsR,signalsP,r,phi,counterR,counterP, iD, iq);
1104  return kStOk;
1105 }
1106 
1107 
1108 
1109 //print the strips around the place where we expect hit
1110 pair<Double_t,Double_t> StFgtGenAVEMaker::getChargeRatio(Float_t r, Float_t phi, Int_t iD, Int_t iq)
1111 {
1112  //first r:
1113  Double_t maxRCharge=-9999;
1114  Double_t maxPhiCharge=-9999;
1115  Int_t maxRInd=-1;
1116  Int_t maxPInd=-1;
1117  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
1118  {
1119  Int_t geoId=pStrips[iD*4+iq][i].geoId;
1120  generalStrip& pStrip=pStrips[iD*4+iq][i];
1121  Short_t disc, quadrant,strip;
1122  Char_t layer;
1123  Double_t ordinate, lowerSpan, upperSpan;//, prvOrdinate;
1124  StFgtGeom::getPhysicalCoordinate(geoId,disc,quadrant,layer,ordinate,lowerSpan,upperSpan);
1125  StFgtGeom::decodeGeoId(geoId,disc, quadrant, layer, strip);
1126  // if(layer=='P' && disc==iD && iq==quadrant)
1127  // cout <<"looking for " << phi << " have: " << ordinate <<" diff: " << fabs(ordinate-phi) <<endl;
1128  if(disc==iD && iq==quadrant && ((layer =='R' && fabs(ordinate-r)<0.7) || (layer=='P' && fabs(ordinate-phi)<0.03) || (layer=='P' && fabs(ordinate-phi+2*MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-2*MY_PI)<0.03)|| (layer=='P' && fabs(ordinate-phi+MY_PI)<0.03 ) || (layer=='P' && fabs(ordinate-phi-MY_PI)<0.03)))
1129  {
1130  if(layer=='P')
1131  {
1132  if(pStrip.charge>maxPhiCharge)
1133  {
1134  maxPhiCharge=pStrip.charge;
1135  maxPInd=i;
1136  }
1137  }
1138  else
1139  {
1140  if(pStrip.charge>maxRCharge)
1141  {
1142  maxRCharge=pStrip.charge;
1143  maxRInd=i;
1144  }
1145  }
1146  }
1147  }
1148  if(maxRInd>=0 && maxPInd>=0)
1149  {
1150  if(maxRCharge>1000 && maxPhiCharge>1000)
1151  {
1153  if(maxRInd>0)
1154  maxRCharge+=pStrips[iD*4+iq][maxRInd-1].charge;
1155  if(maxRInd< (int)(pStrips[iD*4+iq].size()-1))
1156  maxRCharge+=pStrips[iD*4+iq][maxRInd+1].charge;
1157  if(maxPInd>0)
1158  maxPhiCharge+=pStrips[iD*4+iq][maxPInd-1].charge;
1159  if(maxPInd< (int)(pStrips[iD*4+iq].size()-1))
1160  maxPhiCharge+=pStrips[iD*4+iq][maxPInd+1].charge;
1161 
1162  return pair<Double_t,Double_t>(maxRCharge,maxPhiCharge);
1163  }
1164 
1165  }
1166  return pair<Double_t,Double_t>(-9999,-9999);
1167 }
1168 
1169 
1170 
1171 
1172 
1173 
1174 Bool_t StFgtGenAVEMaker::getTrack(vector<AVPoint>& points, Double_t ipZ)
1175 {
1176  // (*outTxtFile) <<"getTrack" <<endl;
1177  // cout <<"get track" <<endl;
1178  ipZ=-9999; //get ourselves
1179  ipZ=vtxZ;
1180  vector<AVPoint>::iterator iter=points.begin();
1181  Double_t A = 0, B = 0, Cx = 0, Cy = 0, D = 0, Ex = 0, Ey = 0;
1182  Double_t dist = 0;
1183 
1184  //LOG_INFO << " --------------------------> calling makeLine with " << points.size() << " 0x" << std::hex << discBitArray << std::dec << endm;
1185  // cout <<"get track2" <<endl;
1186  //do the alignment
1187 
1188 #ifdef COSMIC
1189  float tmpX, tmpY,tmpZ,tmpP,tmpR;
1190 
1191  for( ; iter != points.end(); ++iter ){
1192  getAlign(iter->dID,iter->phi,iter->r,tmpX,tmpY,tmpZ,tmpP,tmpR);
1193  // cout <<"before: " << iter->phi << ", " << iter->r <<" " << iter->x <<" " << iter->y << ", " << iter->z <<endl;
1194  iter->phi=tmpP;
1195  iter->r=tmpR;
1196  iter->x=tmpX;
1197  iter->y=tmpY;
1198  iter->z=tmpZ;
1199  // cout <<"after: " << iter->phi << ", " << iter->r <<" " << iter->x <<" " << iter->y << ", " << iter->z <<endl;
1200  }
1201 #endif
1202 
1203  // cout <<" we have " << points.size() << " points " << endl;
1204  for( iter=points.begin(); iter != points.end(); ++iter ){
1205 
1206  Double_t x = iter->x;
1207  Double_t y = iter->y;
1208  Double_t z = iter->z;
1209  // cout <<"x: " << x << " y: " << y << " z: " << z <<endl;
1210  A += z*z;
1211  B += z;
1212  Cx += x*z;
1213  Cy += y*z;
1214  Ex += x;
1215  Ey += y;
1216 
1217  // cout << "*** Point located at " << x << ' ' << y << ' ' << z << " Disk: " << iter->dID <<endl;
1218  }
1219  // cout <<"ipZ: " << ipZ <<endl;
1220  // cout <<"get track3" <<endl;
1221  D = points.size();
1222  //invalid is -999
1223 #ifdef FIT_WITH_VERTEX
1224  if(muDst && ipZ>-100 && ipZ<100)
1225  {
1226  A+=ipZ*ipZ;
1227  B+=ipZ;
1228  D++;
1229  }
1230 #endif
1231  // cout << "*** Consts " << A << ' ' << B << ' ' << Cx << ' ' << Cy << ' ' << D << ' ' << Ex << ' ' << Ey << endl;
1232  Double_t denom = D*A - B*B;
1233  if( denom ){
1234  Double_t bx = (-B*Cx + A*Ex)/denom;
1235  Double_t by = (-B*Cy + A*Ey)/denom;
1236  Double_t mx = ( D*Cx - B*Ex)/denom;
1237  Double_t my = ( D*Cy - B*Ey)/denom;
1238  // cout <<"bx: " << bx <<" by: " << by <<" mx: " << mx << " my: " << my <<endl;
1239  for( iter = points.begin(); iter != points.end(); ++iter ){
1240  // cout << "--- Location at each disc, " << iter->dID <<" "
1241  // << "X: " << mx*iter->z+bx << " vs " << iter->x << ' '
1242  // << "Y: " << my*iter->z+by << " vs " << iter->y << " " <<" charge r: " << iter->rCharge <<" phi: " << iter->phiCharge <<
1243  // " r: " << iter->r <<" phi: " << iter->phi <<endl;
1244  (*outTxtFile) << "--- Location at each disc, " << iter->dID <<" "
1245  << "X: " << mx*iter->z+bx << " vs " << iter->x << ' '
1246  << "Y: " << my*iter->z+by << " vs " << iter->y << " " <<" charge r: " << iter->rCharge <<" phi: " << iter->phiCharge <<
1247  " r: " << iter->r <<" phi: " << iter->phi <<endl;
1248  };
1249  dist = 0;
1250  // cout <<"get track4" <<endl;
1251  vector<AVPoint> redPoints;
1252  // find closest points...
1253  for(int iDx=0;iDx<6;iDx++)
1254  {
1255  Double_t minDistance=9999;
1256  Int_t pointIdx=-1;
1257  Int_t cnt=0;
1258  for( iter = points.begin(); iter != points.end(); ++iter ){
1259  Int_t dId=iter->dID;
1260  Double_t distX=fabs((iter->z*mx+bx)-(iter->x));
1261  Double_t distY=fabs((iter->z*my+by)-(iter->y));
1262  // cout <<"distX: " << distX <<" distY: " << distY <<endl;
1263  Double_t distance=distX*distX+distY*distY;
1264  // cout << " got distance " << distance << endl;
1265  if((iDx==dId)&&(distance<minDistance))
1266  {
1267  minDistance=distance;
1268  pointIdx=cnt;
1269  // cout <<"min dist now:" << minDistance<<" for this dik " << iDx <<" pointIdx: " << pointIdx <<endl;
1270  }
1271  cnt++;
1272  }
1273  if(pointIdx>=0)
1274  {
1275  // cout <<"pushing back " << pointIdx <<endl;
1276  redPoints.push_back(points[pointIdx]);
1277  }
1278  }//end of looping over discs
1279 
1282  // cout <<"doing refit... " <<endl;
1283 // cout <<"get track5" <<endl;
1284  A=0.0;
1285  B=0.0;
1286  Cx=0.0;
1287  Cy=0.0;
1288  Ex=0.0;
1289  Ey=0.0;
1290 
1291  for(vector<AVPoint>::iterator iterR=redPoints.begin() ; iterR != redPoints.end(); ++iterR )
1292  {
1293  Double_t x = iterR->x;
1294  Double_t y = iterR->y;
1295  Double_t z = iterR->z;
1296 
1297  A += z*z;
1298  B += z;
1299  Cx += x*z;
1300  Cy += y*z;
1301  Ex += x;
1302  Ey += y;
1303  }
1304  D = redPoints.size();
1305 #ifdef REFIT_WITH_VERTEX
1306  if(muDst && ipZ>-100 && ipZ<100)
1307  {
1308  A+=ipZ*ipZ;
1309  B+=ipZ;
1310  D++;
1311  }
1312 #endif
1313  Double_t denom = D*A - B*B;
1314  // cout <<"get track6" <<endl;
1315  if( denom ){
1316  bx = (-B*Cx + A*Ex)/denom;
1317  by = (-B*Cy + A*Ey)/denom;
1318  mx = ( D*Cx - B*Ex)/denom;
1319  my = ( D*Cy - B*Ey)/denom;
1320  }
1321  // cout <<"after refit: bx: " << bx <<" by: " << by <<" mx: " << mx << " my: " << my <<endl;
1322  // cout <<"we have refit line: "<< bx << " by: " << by <<" mx: " << mx << " my: " << my <<endl;
1324 
1325  for(vector<AVPoint>::iterator iterR = redPoints.begin(); iterR != redPoints.end(); ++iterR )
1326  {
1327  Double_t distX, distY;
1328  distX=fabs((iterR->z*mx+bx)-(iterR->x));
1329  distY=fabs((iterR->z*my+by)-(iterR->y));
1330  // cout <<"distX: " << distX <<" distY: " << distY <<endl;
1331  dist += (distX*distX+distY*distY);
1332  // cout <<"adding " << (distX*distX+distY*distY) <<" to chi2: " << endl;
1333  // cout << "*** DistSq " << dist << endl;
1334  }
1335  // cout <<"get track8" <<endl;
1336  // cout <<"dist : " << dist <<" D: "<< D <<endl;
1337  dist/=D;
1338  // cout <<" end chi2: " <<dist <<endl;
1339  m_tracks.push_back(AVTrack(mx,my,bx,by,ipZ,dist));
1340  // cout <<" we have " <<m_tracks.size() <<" track now " <<endl;
1341 
1342  points.clear();
1343  for(vector<AVPoint>::iterator it=redPoints.begin();it!=redPoints.end();it++)
1344  {
1345  points.push_back(*it);
1346  }
1347  // copy to pointer, if one provided
1350 
1351  // #ifdef DEBUG
1352  // cout << "*** Event " << GetEventNumber() << " final distSQ " << dist << endl;
1353  // #endif
1354 
1355  // add to vector, if thes OK
1358 
1359  // cout <<"get track9" <<endl;
1360  //Double_t zIpExpX0=0;//x(D6Pos-(xD6/xD1)*D1Pos)*1/(1-xD6/xD1);
1362  //Double_t zIpExpY0=0;//(D6Pos-yD6/yD1*D1Pos)*1/(1-yD6/yD1);
1363 
1364  for(vector<AVPoint>::iterator iter = points.begin(); iter != points.end(); ++iter ){
1365  // cout << "--- Location at each disc at z: " << iter->z << " "
1366  // << "X: " << mx*iter->z+bx << " vs " << iter->x << ' '
1367  // << "Y: " << my*iter->z+by << " vs " << iter->y << " "
1368  // << " charge phi: " << iter->phiCharge <<" rcharge: "<< iter->rCharge <<endl;
1369 
1370  };
1371  // cout <<endl<<endl;
1372  // cout <<"dist again: " <<dist <<endl;
1373  vector<AVTrack>::iterator it_lastTrack=m_tracks.end();
1374  it_lastTrack--;
1375  pair<double,double> dca=getDca(it_lastTrack);
1376  // cout <<"mx: " << it_lastTrack->mx <<" ax: " << it_lastTrack->ax <<" my: " << it_lastTrack->my <<" ay: " << it_lastTrack->ay <<endl;
1377  Double_t vertZ = ( -( it_lastTrack->mx*it_lastTrack->ax + it_lastTrack->my*it_lastTrack->ay )/(it_lastTrack->mx*it_lastTrack->mx+it_lastTrack->my*it_lastTrack->my));
1378  (it_lastTrack)->trkZ=vertZ;
1379  it_lastTrack->dca=dca.second;
1380  it_lastTrack->ipZ=dca.first;
1381  // cout <<"dca: " << dca.first <<" vertZ: " << vertZ <<endl;
1382  // cout <<"get track10" <<endl;
1383  // cout <<"dist: "<< dist <<" vertZ: " << vertZ <<endl;
1384  if(dist< MAX_DIST_CHI && fabs(vertZ)<VERTEX_CUT)// && fabs(bx)<40 && fabs(by)<40)
1385  {
1386  // cout <<" track accepted " <<endl;
1387  // cout <<"get track10-10" <<endl;
1388  // cout <<"track cuts passed " <<endl;
1389  set<Short_t> disksHit;
1390  // cout <<"track has " << points.size() <<" points " <<endl;
1391  for(vector<AVPoint>::iterator iterP = points.begin(); iterP != points.end(); ++iterP ){
1392  // cout <<"testing" << endl;
1393  // cout <<"get track10-3, " << iterP->dID<<" quad: " << iterP->quadID << endl;
1394  // cout <<" filling disk " << iterP->dID <<" quad " << iterP->quadID <<" with r charge: " << iterP->rCharge <<" phic " << iterP->phiCharge<<endl;
1395  // if(iterP->r>20)
1396  {
1397  chargeCorr[iterP->dID*4+iterP->quadID]->Fill(iterP->rCharge,iterP->phiCharge);
1398  clusterSizeP[iterP->dID*4+iterP->quadID]->Fill(iterP->phiSize);
1399  clusterSizeR[iterP->dID*4+iterP->quadID]->Fill(iterP->rSize);
1400  }
1401  // cout <<"get track10-31" <<endl;
1402  h_clusterChargeR[iterP->dID]->Fill(iterP->rCharge);
1403  // cout <<"get track10-4" <<endl;
1404  h_clusterChargePhi[iterP->dID]->Fill(iterP->phiCharge);
1405  h_clusterSizeR[iterP->dID]->Fill(iterP->rSize);
1406  // cout <<"get track10-5" <<endl;
1407  h_clusterSizePhi[iterP->dID]->Fill(iterP->phiSize);
1408  disksHit.insert(iterP->dID);
1409  // cout <<"get track10-6" <<endl;
1410  }
1411  // cout <<"get track11" <<endl;
1412  for(int i=0;i<6;i++)
1413  {
1414  Double_t xExp=mx*getLocDiscZ(i)+bx;
1415  Double_t yExp=my*getLocDiscZ(i)+by;
1416  // cout <<"expecting x: " << xExp << ", " << yExp<<endl;
1417  Int_t quad=-1;
1418  //x=r*cos(phi)
1419  //y=r*sin(phi)
1420  Double_t r=sqrt(xExp*xExp+yExp*yExp);
1421  //if both, x and y, are negative we get the wrong angle
1422  Double_t phi=atan(yExp/xExp);
1423  if(xExp <0 && yExp <0)
1424  phi-=TMath::Pi();
1425  // cout <<"phi: "<<phi << " from x: " << xExp <<" y: " << yExp<<endl;
1426  if(phi<-MY_PI)
1427  phi+=2*MY_PI;
1428  if(phi>MY_PI)
1429  phi-=2*MY_PI;
1430  // if(phi<-MY_PI)
1431  // phi+=2*MY_PI;
1432 
1433  // quad=getQuadFromCoo(xExp,yExp);
1434  quad=StFgtGeom::getQuad(phi);
1435  // cout <<"got quad : " << quad << " for phi: " << phi <<endl;
1436  //convert to phi in quad.., so we have to subtract that axis...
1437  phi-=StFgtGeom::phiQuadXaxis(quad);
1438 
1439  if(phi>TMath::Pi())
1440  phi-=(2*TMath::Pi());
1441  if(phi<((-1)*TMath::Pi()))
1442  phi+=(2*TMath::Pi());
1443 
1444  if(phi<0)
1445  phi+=MY_PI;
1446  (*outTxtFile) << " looking at Track with chi2/ndf *[cm}: " << it_lastTrack->chi2 << " z vertex: " << it_lastTrack->ipZ << endl;
1447 
1448  //this snipplet is only good for the disks in which we don't look for efficiencies, otherwise it overcounts the nonEff
1449  //so either use the findClosest point for all, or check if we look at efficient disk or not...
1450  /* if(disksHit.find(i)!=disksHit.end())
1451  {
1452  radioPlotsEff[i]->Fill(xExp,yExp);
1453  }
1454  else//not efficient
1455  {
1456  radioPlotsNonEff[i]->Fill(xExp,yExp);
1457 
1458  }*/
1460 
1461  // cout <<"fill strip histos " << endl;
1462  fillStripHistos(r,phi,i,quad);
1463  if(isSomewhatEff(r,phi,i,quad))
1464  {
1465  // cout <<" somewhat eff: " << xExp <<" y: " << yExp <<" disk: " << i <<endl;
1466  radioPlotsEffLoose[i]->Fill(xExp,yExp);
1467  }
1468  else
1469  {
1470  // cout <<"not somewhat eff: " << xExp <<" y: " << yExp <<endl;
1471  radioPlotsNonEffLoose[i]->Fill(xExp,yExp);
1472  }
1473 
1474  if(i==m_effDisk && quad==QUAD_EFF)
1475  {
1476  //do the r/phi thing
1477  if(findClosestStrip('R',r,i,quad)<MAX_DIST_STRIP_R)
1478  radioPlotsEffR[i]->Fill(xExp,yExp);
1479  else
1480  radioPlotsNonEffR[i]->Fill(xExp,yExp);
1481  //correct phi again, so it is absolute in space, like the quads
1482  if(findClosestStrip('P',phi+StFgtGeom::phiQuadXaxis(quad),i,quad)<MAX_DIST_STRIP_PHI)
1483  {
1484  // cout <<" found phi " <<endl;
1485  radioPlotsEffPhi[i]->Fill(xExp,yExp);
1486  }
1487  else
1488  {
1489  // cout << " not found, dist is:" <<findClosestStrip('P',phi,i,quad) <<endl;
1490  radioPlotsNonEffPhi[i]->Fill(xExp,yExp);
1491  }
1492  }
1493 
1494  Double_t closestPoint=findClosestPoint(mx,bx,my,by,xExp,yExp,i);
1495  // cout <<" closest point is " << closestPoint <<" away " <<endl;
1496  // (*outTxtFile) <<"closest point t " << xExp <<" , " << yExp << " is : " << closestPoint << " away " << endl;
1497 
1498  if(findClosestPoint(mx,bx,my,by,xExp,yExp,i)<MAX_DIST2_EFF)
1499  {
1500  // cout <<"found point on eff disk, x: " << xExp <<" y: " << yExp <<", dist: " <<closestPoint <<endl;
1501  radioPlotsEff[i]->Fill(xExp,yExp);
1502  if(i==m_effDisk)
1503  // if(i==1)
1504  {
1505  hResidua->Fill(sqrt(closestPoint));
1506  hResiduaX->Fill(xExp,sqrt(closestPoint));
1507  hResiduaY->Fill(yExp,sqrt(closestPoint));
1508  hResiduaR->Fill(r,sqrt(closestPoint));
1509  hResiduaP->Fill(phi,sqrt(closestPoint));
1510  }
1511 
1512  (*outTxtFile) <<"***** found hit in disk " <<i << " quad: " << quad << " at " << xExp<<", " << yExp<<" r: " << r <<" phi: " <<phi << endl;
1513 
1514 #ifdef DO_PRINT
1515  // if(i==m_effDisk)
1516  printArea(r,phi,i,quad);
1517 #endif
1518  // chargeCorrInEffDisk->Fill(rPhiRatio.first,rPhiRatio.second);
1519  }
1520  else
1521  {
1522  // cout <<"non eff disk, x: " << xExp <<" y: " << yExp <<endl;
1523  radioPlotsNonEff[i]->Fill(xExp,yExp);
1524  // if(i==m_effDisk)
1525  (*outTxtFile) <<"expected (but haven't found) point on disk " << i <<", x: " << xExp <<" y: " << yExp << " r: " << r <<" phi: " << phi << " quad:: " << quad << endl;
1527 #ifdef DO_PRINT
1528  // if(i==m_effDisk)
1529  printArea(r,phi,i,quad);
1530 #endif
1531  }
1532  pair<Double_t,Double_t> rPhiRatio=getChargeRatio(r,phi,i,quad);
1533 
1534  if(rPhiRatio.first>0 && rPhiRatio.second>0)
1535  {
1536  double asym=fabs((Double_t)1-rPhiRatio.first/rPhiRatio.second);
1537  double ratio=rPhiRatio.first/rPhiRatio.second;
1538 
1539  // cout <<" filling with : r charge: " << rPhiRatio.first <<" , " << rPhiRatio.second <<" ratio: " << ratio <<" asym: " << asym<<endl;
1540  if(asym<2 && ratio <2)
1541  {
1542  hChargeAsym->Fill(asym);
1543  hChargeRatio->Fill(ratio);
1544  chargeRatioInEffDisk->Fill(xExp,yExp,ratio);
1545  chargeAsymInEffDisk->Fill(xExp,yExp,asym);
1546  }
1547  }
1548 
1549  }
1550  }
1551  // cout <<" returning true " <<endl;
1552  return true;
1553  }
1554  // cout <<" false... " <<endl;
1555  return false;
1556 };
1557 
1558 Double_t StFgtGenAVEMaker::getRPhiRatio(vector<generalCluster>::iterator hitIterBegin, vector<generalCluster>::iterator hitIterEnd)
1559 {
1560  Short_t quad;
1561  Char_t layer;
1562  Int_t numR=0;
1563  Int_t numPhi=0;
1564  vector<generalCluster>::iterator hitIter=hitIterBegin;
1565  for(;hitIter!=hitIterEnd;hitIter++)
1566  {
1567  quad=hitIter->quad;
1568  layer=hitIter->layer;
1569 
1570  if(layer=='R')
1571  numR++;
1572  else
1573  numPhi++;
1574  }
1575 
1576  if(numR+numPhi>0)
1577  return (numR-numPhi)/((Double_t)(numR+numPhi));
1578  else
1579  return -1;
1580 };
1581 
1585 Int_t StFgtGenAVEMaker::Make()
1586 {
1587  // cout <<"ave make " <<endl;
1588  Int_t ierr = kStOk;
1589  (*outTxtFile) <<"----------------------------- Event Nr: " << evtNr<<" -----------------" <<endl;
1590  StFgtGeneralBase::Make();
1591  unsigned int oldNumTracks=m_tracks.size();
1592 
1593  // cout <<"eff disk : " << m_effDisk <<endl;
1594 
1595  // cout <<" mtracks size: " << m_tracks.size() << " oldnum: "<< oldNumTracks <<endl;
1596  for(int i=0;i<6;i++)
1597  {
1598  // cout <<"there are " << pClusters[i]->size() << " clusters in disk " << i <<endl;
1599  for(int j=0;j<pClusters[i]->size();j++)
1600  {
1601  /* if((*(pClusters[i]))[j].layer=='R')
1602  cout <<"R layer, ";
1603  else
1604  cout <<"Phi layer, ";*/
1605  Int_t seedType=(*(pClusters[i]))[j].seedType;
1606  Double_t posPhi=(*(pClusters[i]))[j].posPhi;
1607  Double_t posR=(*(pClusters[i]))[j].posR;
1608  Int_t clusSize=(*(pClusters[i]))[j].clusterSize;
1609  Double_t charge=(*(pClusters[i]))[j].clusterCharge;
1610  Int_t cntGeoId=(*(pClusters[i]))[j].centralStripGeoId;
1611  // cout <<"cluster pos phi: " << posPhi <<" posR: " << posR <<" clusSize: " << clusSize << " charge: "<< charge <<" geoId: "<< cntGeoId <<" seedT : " << seedType<<endl;
1612  }
1613  }
1614 
1615  for(int iD=0;iD<6;iD++)
1616  {
1617  (*outTxtFile) <<" In Disc " << iD << " we have clusters with geo id: ";
1618 
1619  // cout <<" In Disc " << iD << " we have clusters with geo id: ";
1620  vector<generalCluster>::iterator hitIter;
1621  vector<generalCluster> &hitVec=*(pClusters[iD]);
1622 
1623  for(hitIter=hitVec.begin();hitIter != hitVec.end();hitIter++)
1624  {
1625  (*outTxtFile) << hitIter->centralStripGeoId << ", ";
1626  clusterGeoId[iD]->Fill(hitIter->centralStripGeoId);
1627  //Adding to fill APV noise histograms for disk 1, quad A, RDO=1, ARM=0, APV 0-4, 17-21
1628  Int_t rdo, arm, apv, chan;
1629  mDb->getElecCoordFromGeoId(hitIter->centralStripGeoId, rdo,arm,apv,chan);
1630  //cout<<" WHAT IS APV="<<apv<<" rdo="<<rdo<<" iD="<<iD<<" arm="<<arm<<endl;
1631  if(iD==0 && rdo==1 && arm==0){
1632  disk1QuadA[apv]->Fill(chan);
1633  }
1634  if(hitIter->layer=='R')
1635  {
1636  (*outTxtFile) <<"R ordinate: " << hitIter->posR<<endl;
1637  //added this line under hitIter layer R
1638  clustersR[iD]->Fill(hitIter->posR);
1639  }
1640  else
1641  {
1642  Double_t phiQ=StFgtGeom::phiQuadXaxis(hitIter->quad);
1643  (*outTxtFile) <<"Phi ordinate: " << hitIter->posPhi-phiQ<<endl;
1644  //added this line under phi ordinate
1645  clustersP[iD]->Fill(hitIter->posPhi);
1646  }
1647 
1648  }
1649  (*outTxtFile)<<endl;
1650 
1651  for(int iq=0;iq<4;iq++)
1652  {
1653  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
1654  {
1655  generalStrip& pStrip=pStrips[iD*4+iq][i];
1656  if(pStrip.seedType==kFgtSeedType1|| pStrip.seedType==kFgtSeedType2 || pStrip.seedType==kFgtSeedType3)
1657  {
1658 
1659  for(hitIter=hitVec.begin();hitIter != hitVec.end();hitIter++)
1660  {
1661  if(abs(hitIter->centralStripGeoId-pStrip.geoId)<2)
1662  {
1663  (*outTxtFile) <<"found cluster with geo id: " << hitIter->centralStripGeoId;
1664  Double_t phiQ=StFgtGeom::phiQuadXaxis(iq);
1665  if(hitIter->layer=='R')
1666  (*outTxtFile) <<"R ordinate: " << hitIter->posR<<endl;
1667  else
1668  (*outTxtFile) <<"Phi ordinate: " << hitIter->posPhi-phiQ<<endl;
1669  }
1670  }
1671 #ifdef PRINT_1D
1672  printArea1D(iD,iq,pStrip.geoId);
1673 #endif
1674  }
1675  }
1676  }
1677  }
1678 
1679 
1680  Float_t x;
1681  Float_t y;
1682  // if(vtxRank<=1)
1683  // return kStOk;
1684 
1685  //look at the r phi ratio for each disk
1686  // for(int i=0;i<6;i++)
1687  // {
1688  // Double_t ratio=getRPhiRatio(pClusters[i]->begin(),pClusters[i]->end());
1689  // rPhiRatioPlots[i]->Fill(ratio);
1690  // cout << "ratio for disk: " << i << " is " << ratio <<" disk has: " << tmpClusterCol->getHitVec().size() << "hits" <<endl;
1691  // }
1692 
1693  //vector<generalCluster> &hitVecD1=*(pClusters[0]);
1694  //vector<generalCluster> &hitVecD6=*(pClusters[5]);
1695 
1696 
1697  for(int iD=0;iD<6;iD++)
1698  {
1699  // cout << " there are " << pClusters[iD]->size() <<" cluster in disk : " << iD+1 <<endl;
1700  int clusCounts[8];
1701  memset(clusCounts,0,sizeof(int)*8);
1702  vector<generalCluster> &tHitVec=*(pClusters[iD]);
1703  vector<generalCluster>::iterator tHit=tHitVec.begin();
1704  for(;tHit!=tHitVec.end();tHit++)
1705  {
1706 
1707  if(!useChargeMatch || tHit->hasMatch)
1708  {
1709  if(tHit->quad<3)
1710  {
1711  if(tHit->layer=='R')
1712  clusCounts[tHit->quad]++;
1713  else
1714  clusCounts[tHit->quad+4]++;
1715  }
1716  }
1717  }
1718  for(int iq=0;iq<4;iq++)
1719  {
1720  numClustersR[iD*4+iq]->Fill(clusCounts[iq]);
1721  numClustersPhi[iD*4+iq]->Fill(clusCounts[iq+4]);
1722  }
1723 
1724  int numClus=0;
1725  int numPulses=0;
1726 
1727  for(int iq=0;iq<4;iq++)
1728  {
1729  for(unsigned int i=0;i< pStrips[iD*4+iq].size();i++)
1730  {
1731  if(pStrips[iD*4+iq][i].charge>1000)
1732  {
1733  numClus++;
1734  }
1735  if(validPulse(pStrips[iD*4+iq][i]))
1736  numPulses++;
1737  }
1738  // cout <<"disk " << iD+1 << " quad " << iq <<" has " << numClus << " charges and " << numPulses <<" pulses " << endl;
1739  }
1740 
1741  }
1742  set<Int_t> usedPoints;//saves the points that have been used for tracks (and shouldn't be reused)
1743  Double_t D1Pos=getLocDiscZ(0);
1744  Double_t D6Pos=getLocDiscZ(5);
1745  Double_t zArm=D6Pos-D1Pos;
1746  vector<generalCluster>::iterator hitIterD1,hitIterD6, hitIterD1R, hitIterD6R, hitIter, hitIter2;
1747 
1748  int clusCountsTemp[24];
1749  int clusCounts[24];
1750  memset(clusCounts,0,sizeof(int)*24);
1751  memset(clusCountsTemp,0,sizeof(int)*24);
1752 
1753  for(int iSeed1=0;iSeed1<5;iSeed1++)
1754  {
1755  for(int iSeed2=iSeed1+1;iSeed2<6;iSeed2++)
1756  {
1757  // cout<< " using " << iSeed1 <<" and " << iSeed2 << " as seeds " << endl;
1758  // (*outTxtFile) << " using " << iSeed1 <<" and " << iSeed2 << " as seeds " << endl;
1759  if((iSeed2-iSeed1)<1)//to have large enough lever arm. Also, since three points are required shouldn't matter?
1760  continue;
1761  if(iSeed1==m_effDisk || iSeed2==m_effDisk)
1762  continue;
1763  if(pClusters[iSeed1]->size() > MAX_CLUSTERS || pClusters[iSeed2]->size() > MAX_CLUSTERS)
1764  {
1765  // cout <<"too many clusters in the disk!!!"<<endl<<endl;
1766  continue;
1767  }
1768  //track where we have hits in disks
1769 
1770  // cout <<"using " << iSeed1 << " and " << iSeed2 << " as seed " <<endl;
1771  vector<generalCluster> &hitVecSeed1=*(pClusters[iSeed1]);
1772  vector<generalCluster> &hitVecSeed2=*(pClusters[iSeed2]);
1773 
1774  D1Pos=getLocDiscZ(iSeed1);
1775  D6Pos=getLocDiscZ(iSeed2);
1776  zArm=D6Pos-D1Pos;
1777 
1778  for(hitIterD1=hitVecSeed1.begin();hitIterD1 != hitVecSeed1.end();hitIterD1++)
1779  {
1780  //this is from the loose clustering and the cluster doesn't have energy match
1781  if(useChargeMatch && !hitIterD1->hasMatch)
1782  continue;
1783  Short_t quadP=hitIterD1->quad;
1784  //Short_t disc=hitIterD1->disc;
1785  //Short_t strip=hitIterD1->strip;
1786  Char_t layer=hitIterD1->layer;
1787 
1788  Double_t seed1ChargePhi=hitIterD1->CHARGE_MEASURE;
1789  Double_t seed1SizePhi=hitIterD1->clusterSize;
1790  // if(quadP>2)
1791  // continue;
1792 
1793  //do 1D 'fit' with r strips and the (x,y) thing
1794  Int_t geoIdSeed1=hitIterD1->centralStripGeoId;
1795  if(usedPoints.find(geoIdSeed1)!=usedPoints.end())
1796  continue;
1797  if(layer!='P')
1798  continue;
1799  // cout <<"ave make1 " <<endl;
1800  Float_t phiD1=hitIterD1->posPhi;
1801  for(hitIterD6=hitVecSeed2.begin();hitIterD6 != hitVecSeed2.end();hitIterD6++)
1802  {
1803  if(useChargeMatch && !hitIterD6->hasMatch)
1804  continue;
1805 
1806  Int_t geoIdSeed2=hitIterD6->centralStripGeoId;
1807  Short_t quadP_2=hitIterD6->quad;
1808  //Short_t disc=hitIterD6->disc;
1809  //Short_t strip=hitIterD6->strip;
1810  Char_t layer=hitIterD6->layer;
1811  Double_t seed2ChargePhi=hitIterD6->CHARGE_MEASURE;
1812  Double_t seed2SizePhi=hitIterD6->clusterSize;
1813  if(layer!='P')
1814  continue;
1815  if(usedPoints.find(geoIdSeed2)!=usedPoints.end())
1816  continue;
1817  Float_t phiD6=hitIterD6->posPhi;
1818  if(fabs(phiD6-phiD1)>MAX_PHI_DIFF)
1819  continue;
1820 
1821  for(hitIterD1R=hitVecSeed1.begin();hitIterD1R != hitVecSeed1.end();hitIterD1R++)
1822  {
1823  if(useChargeMatch && !hitIterD1R->hasMatch)
1824  continue;
1825  Int_t geoIdSeed1R=hitIterD1R->centralStripGeoId;
1826  Short_t quadR=hitIterD1R->quad;
1827  //Short_t disc=hitIterD1R->disc;
1828  //Short_t strip=hitIterD1R->strip;
1829  Char_t layer=hitIterD1R->layer;
1830  Double_t seed1ChargeR=hitIterD1R->CHARGE_MEASURE;
1831  Double_t seed1SizeR=hitIterD1R->clusterSize;
1832  Int_t quadSeed1=-1;
1833  if(layer!='R')
1834  continue;
1835  if(quadR!=quadP)
1836  continue;
1837  quadSeed1=quadR;
1838  if(usedPoints.find(geoIdSeed1R)!=usedPoints.end())
1839  continue;
1840 
1841  Float_t rD1=hitIterD1R->posR;
1842  Float_t xD1=rD1*cos(phiD1);
1843  Float_t yD1=rD1*sin(phiD1);
1844  // cout <<"disk: " << iSeed1<<", phiD1: " << phiD1 <<" xD1: " << xD1 <<" yD1: " << yD1 <<" rD1: " << rD1 <<endl;
1845  // cout <<"ave make3 " <<endl;
1846 
1847  for(hitIterD6R=hitVecSeed2.begin();hitIterD6R != hitVecSeed2.end();hitIterD6R++)
1848  {
1849  if(useChargeMatch && !hitIterD6R->hasMatch)
1850  continue;
1851  Int_t geoIdSeed2R=hitIterD6R->centralStripGeoId;
1852  Short_t quadR_2=hitIterD6R->quad;
1853  //Short_t disc=hitIterD6R->disc;
1854  //Short_t strip=hitIterD6R->strip;
1855  Char_t layer=hitIterD6R->layer;
1856  Double_t seed2ChargeR=hitIterD6R->CHARGE_MEASURE;
1857 
1858  Double_t seed2SizeR=hitIterD6R->clusterSize;
1859  Int_t quadSeed2=-1;
1860 
1861  if(quadP_2!=quadR_2)
1862  continue;
1863  quadSeed2=quadP_2;
1864  if(layer!='R')
1865  continue;
1866 
1867  clusCountsTemp[iSeed1*4+quadR]++;
1868  clusCountsTemp[iSeed2*4+quadR_2]++;
1869  if(usedPoints.find(geoIdSeed2R)!=usedPoints.end())
1870  continue;
1871  Float_t rD6=hitIterD6R->posR;
1872  //track goes towards smaller radii
1873 #ifndef COSMIC
1874  if(rD1>rD6)
1875  continue;
1876 #endif
1877 
1878  vector<AVPoint> v_points;
1879  //add the seed points to the points
1880 
1881  Double_t xD6=rD6*cos(phiD6);
1882  Double_t yD6=rD6*sin(phiD6);
1883  // cout <<"Disk " << iSeed2 <<", phiD6: " << phiD6 <<" xD6: " << xD6 <<" yD6: " << yD6 <<" rD6: " << rD6 <<endl;
1884  v_points.push_back(AVPoint(xD1,yD1,D1Pos,rD1,phiD1,iSeed1,quadSeed1,seed1ChargeR, seed1ChargePhi, seed1SizeR, seed1SizePhi));
1885  v_points.push_back(AVPoint(xD6,yD6,D6Pos,rD6,phiD6,iSeed2,quadSeed2,seed2ChargeR, seed2ChargePhi, seed2SizeR, seed2SizePhi));
1887  vector< pair<Int_t,Double_t> > v_x;
1888  vector< pair<Int_t,Double_t> > v_y;
1889  vector< pair<Int_t,Double_t> > v_r;
1890 
1891  vector< pair<Int_t,Double_t> > v_xFail;
1892  vector< pair<Int_t,Double_t> > v_yFail;
1893  vector< pair<Int_t,Double_t> > v_rFail;
1894 
1895  vector< pair< Int_t, Double_t> > v_rCharge;
1896  vector< pair< Int_t, Double_t> > v_phiCharge;
1897 
1898  vector<Int_t> v_clusterSizeR;
1899  vector<Int_t> v_clusterSizePhi;
1900 
1901  vector<Int_t> v_geoIDsR;
1902  vector<Int_t> v_geoIDsPhi;
1903 
1904  int iFound=0;
1905  int iFoundR=0;
1906  // cout <<"ave make4 " <<endl;
1907  //zarm is d6 position - D1
1908  //Double_t xIpExp=xD1+(xD6-xD1)*(-D1Pos)/zArm;
1909  //Double_t yIpExp=yD1+(yD6-yD1)*(-D1Pos)/zArm;
1911  //Double_t zIpExpX0=(D6Pos-(xD6/xD1)*D1Pos)*1/(1-xD6/xD1);
1913  //Double_t zIpExpY0=(D6Pos-yD6/yD1*D1Pos)*1/(1-yD6/yD1);
1914  // cout <<"
1915  //now find other points
1916  vector<generalCluster>::iterator iterClosestPhi;
1917  vector<generalCluster>::iterator iterClosestR;
1918 
1919  Double_t closestDist=999999;
1920  Int_t quadTestR=-1;
1921  for(int iD=0;iD<6;iD++)
1922  {
1923  // cout <<"testting disk: " << iD <<endl;
1924  if(iD==iSeed1 || iD==iSeed2 || iD==m_effDisk)
1925  continue;
1926  // cout <<"not seed " << endl;
1927  Bool_t found=false;
1928  Bool_t foundR=false;
1929  //check for hit
1930  Double_t diskZ=getLocDiscZ(iD);
1931  //expected
1932 
1933  Double_t xPosExp=xD1+(xD6-xD1)*(diskZ-D1Pos)/zArm;
1934  Double_t yPosExp=yD1+(yD6-yD1)*(diskZ-D1Pos)/zArm;
1935  // cout <<"using disk: " << iD <<" diskZ: " << diskZ <<endl;
1936  // cout <<"hope to see something x: " << xPosExp <<" y: " << yPosExp <<endl;
1937  // cout <<"phi1: " << phiD1 <<" phi6: " << phiD6 <<endl;
1938  Double_t rPosExp=rD1+(rD6-rD1)*(diskZ-D1Pos)/zArm;
1939  vector<generalCluster> &hitVec=*(pClusters[iD]);
1940  for(hitIter=hitVec.begin();hitIter!=hitVec.end();hitIter++)
1941  {
1942  if(useChargeMatch && !hitIter->hasMatch)
1943  continue;
1944  //do 1D 'fit' with r strips and the (x,y) thing
1945  Int_t geoIdPhi=hitIter->centralStripGeoId;
1946  Short_t quad=hitIter->quad;
1947  Int_t quadTestPhi=quad;
1948  Short_t disc=hitIter->disc;
1949  Short_t strip=hitIter->strip;
1950  Char_t layer=hitIter->layer;
1951 
1952  if(usedPoints.find(geoIdPhi)!=usedPoints.end())
1953  continue;
1954  if(layer!='P')
1955  continue;
1956  Float_t phi=hitIter->posPhi;
1957 
1958  if(fabs(phi-phiD1)>MAX_PHI_DIFF)
1959  continue;
1960  if(fabs(phi-phiD6)>MAX_PHI_DIFF)
1961  continue;
1962 
1963  //Double_t phiCharge=hitIter->CHARGE_MEASURE;
1964 
1965  // Double_t phiCharge=hitIter->maxAdc;
1966  // Int_t clusterSizePhi=hitIter->clusterSize;
1967  // if(clusterSizePhi<=1)
1968  // continue;
1969  // cout <<"ave make5 " <<endl;
1970  for(hitIter2=hitVec.begin();hitIter2!=hitVec.end();hitIter2++)
1971  {
1972  if(useChargeMatch && !hitIter2->hasMatch)
1973  continue;
1974  Int_t geoIdR=hitIter2->centralStripGeoId;
1975  StFgtGeom::decodeGeoId(geoIdR,disc, quad, layer, strip);//ok
1976  // cout <<" r? " <<endl;
1977  if(usedPoints.find(geoIdR)!=usedPoints.end())
1978  continue;
1979  // cout <<"not used yet " <<endl;
1980  if(layer!='R')
1981  continue;
1982  quadTestR=quad;
1983  if(quadTestR!=quadTestPhi)
1984  continue;
1985  Float_t r=hitIter2->posR;
1986  //make sure that the radius makes sense for a track that goes from inner to outer radious
1987 #ifndef COSMIC
1988  if(r>rD6 || r<rD1)
1989  continue;
1990 #endif
1991  x=r*cos(phi);
1992  y=r*sin(phi);
1993  // cout <<"checking with x: " << x << " y: " << y << " phi: " << phi <<" r: " << r <<endl;
1994  // cout <<" x, y: " << x <<", " << y << " exp: " << xPosExp << " , " << yPosExp <<endl;
1995  Double_t dist2=(x-xPosExp)*(x-xPosExp)+(y-yPosExp)*(y-yPosExp);
1996  // cout <<" dist2: " << dist2 <<endl;
1997 
1998 #ifdef ADD_MULTIPLE
1999  if(dist2<MAX_DIST2)
2000  {
2001  Double_t rCharge=hitIter2->CHARGE_MEASURE;
2002  Double_t phiCharge=hitIter->CHARGE_MEASURE;
2003  Int_t clusterSizeR=hitIter2->clusterSize;
2004  Int_t clusterSizePhi=hitIter->clusterSize;
2005  v_points.push_back(AVPoint(x,y,diskZ,r,phi,iD,quadTestR, rCharge,phiCharge, clusterSizeR,clusterSizePhi));
2006  }
2007 
2008 #endif
2009  if(dist2<closestDist)
2010  {
2011  closestDist=dist2;
2012  iterClosestPhi=hitIter;
2013  iterClosestR=hitIter2;
2014  }
2015  }
2016  }
2017  // cout <<"closest dist for disk: "<< iD <<" is : " << closestDist <<endl;
2018  // cout <<"ave make6 " <<endl;
2019  if(closestDist<1000 && closestDist<MAX_DIST2)
2020  {
2021  found=true;
2022  clusCountsTemp[iD*4+quadTestR]++;
2023 
2024  // cout <<"accepted " <<endl;
2025  double r=iterClosestR->posR;
2026  double phi=iterClosestPhi->posPhi;
2027  // cout <<"found closest phi in disk " << iD <<" : " << phi << ", r: " << r << " x: " << x <<" y: " << y << ", this is good enough, distance: " << closestDist <<endl;
2028  Int_t geoIdR=iterClosestR->centralStripGeoId;
2029  Int_t geoIdPhi=iterClosestPhi->centralStripGeoId;
2030 
2031  double x=r*cos(phi);
2032  double y=r*sin(phi);
2033 
2034  Double_t rCharge=iterClosestR->CHARGE_MEASURE;
2035  Double_t phiCharge=iterClosestPhi->CHARGE_MEASURE;
2036  Int_t clusterSizeR=iterClosestR->clusterSize;
2037  Int_t clusterSizePhi=iterClosestPhi->clusterSize;
2038  // cout <<"charge R of middle disk: "<< iD<<": "<< rCharge <<" phicharge: " << phiCharge<<endl;
2039 
2040  v_x.push_back(pair<Int_t,Double_t>(iD,x));
2041  v_y.push_back(pair<Int_t,Double_t>(iD,y));
2042  v_rCharge.push_back(pair<Int_t, Double_t>(iD,rCharge));
2043  v_phiCharge.push_back(pair<Int_t, Double_t>(iD,phiCharge));
2044  //to avoid double counting, t
2045  v_geoIDsPhi.push_back(geoIdPhi);
2046  v_geoIDsR.push_back(geoIdR);
2047  v_clusterSizeR.push_back(clusterSizeR);
2048  v_clusterSizePhi.push_back(clusterSizePhi);
2049  v_points.push_back(AVPoint(x,y,diskZ,r,phi,iD,quadTestR, rCharge,phiCharge, clusterSizeR,clusterSizePhi));
2050  // cout<<" adding point with r: "<< r <<" phi: " << phi <<" x: " << x <<" y: " << y <<endl;
2051  }
2052  // cout <<"ave make8 " <<endl;
2053  //only one per disk
2054  if(found)
2055  iFound++;
2056  else
2057  {
2058  v_xFail.push_back(pair<Int_t,Double_t>(iD,xPosExp));
2059  v_yFail.push_back(pair<Int_t,Double_t>(iD,yPosExp));
2060  }
2061  if(foundR)
2062  iFoundR++;
2063  else
2064  {
2065  // cout <<"failed to find r " << rPosExp<<endl;
2066  v_rFail.push_back(pair<Int_t, Double_t>(iD,rPosExp));
2067  }
2068  closestDist=999999;
2069 
2070  }
2071 
2072  // cout <<"ave make9 " <<endl;
2073  // 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
2074  //with hard cuts on the track we can get the whole vertex distribution
2075  // cout << " we found " << iFound <<" points " << endl;
2076  if(iFound>=(MIN_NUM_POINTS-2)) //at least one hit plus seed and pointing roughly to the interaction region
2077  {
2078  // cout <<"found " <<endl;
2079  Bool_t validTrack=false;
2080  // if(v_x.size()>iFound)
2081  {
2082  Double_t ipZ;
2083  // cout <<"about to get track" <<endl;
2084  // cout <<"check for valid track " << endl;
2085  validTrack=getTrack(v_points, ipZ);
2086 
2087  // cout <<"found 2" <<endl;
2088  if(validTrack)
2089  {
2090  for(Int_t iD=0;iD<kFgtNumDiscs;iD++){
2091  for(int iq=0;iq<4;iq++){
2092  clusCounts[iD*4+iq]+=clusCountsTemp[iD*4+iq];
2093  }}
2094 
2095 
2096  // cout <<"track was valid, phi1: " << phiD1 <<" phiD6: " << phiD6;
2097  // cout <<", rD1: " << rD1 <<" rD6: " << rD6<<endl;
2098  // cout <<"was valid " << endl;
2099  }
2100  }
2101 
2102  // cout <<"found4 " <<endl;
2103  // else
2104  {
2105  // cout <<"ave make8]10 " <<endl;
2106  // cout <<"filled hip with " << xIpExp << " / " << yIpExp <<endl;
2107  for(unsigned int i=0;i<v_x.size();i++)
2108  {
2109  if(validTrack)
2110  {
2111  //at least don't duplicate seeds..., the other ones might belong to multiple tracks
2112  usedPoints.insert(geoIdSeed1);
2113  usedPoints.insert(geoIdSeed1R);
2114  usedPoints.insert(geoIdSeed2);
2115  usedPoints.insert(geoIdSeed2R);
2116  // usedPoints.insert(v_geoIDsPhi[i]);
2117  // usedPoints.insert(v_geoIDsR[i]);
2118  // Int_t disk=v_x[i].first;
2119  // Double_t x=v_x[i].second;
2120  // Double_t y=v_y[i].second;
2121  // Double_t rCharge=v_rCharge[i].second;
2122  // Double_t phiCharge=v_phiCharge[i].second;
2123  // cout <<"filling disk: " << disk <<" with: " << x <<" / " <<y <<endl;
2124  // radioPlotsEff[disk]->Fill(x,y);
2125  // chargeCorr[disk]->Fill(rCharge,phiCharge);
2126  // h_clusterSizeR[disk]->Fill(v_clusterSizeR[i]);
2127  // h_clusterSizePhi[disk]->Fill(v_clusterSizePhi[i]);
2128  // h_clusterChargeR[disk]->Fill(rCharge);
2129  // h_clusterChargePhi[disk]->Fill(phiCharge);
2130  }
2131  }
2132  hitCounter++;
2133  }
2134  }
2135  // cout <<"ave make8]11 " <<endl;
2136 
2137  //start over
2138  iFound=0;
2139  iFoundR=0;
2140  v_x.clear();
2141  v_points.clear();
2142  v_y.clear();
2143  v_r.clear();
2144  v_xFail.clear();
2145  v_yFail.clear();
2146  v_rFail.clear();
2147  memset(clusCountsTemp,0,sizeof(int)*24);
2148  }
2149 
2150  }
2151 
2152  }
2153 
2154  }
2155 
2156  }
2157 
2158  }
2159  for(Int_t iD=0;iD<kFgtNumDiscs;iD++){
2160  for(int iq=0;iq<4;iq++){
2161  // cout <<"filling clus2: " << iD<<" iq: " << iq << " overall: " << iD*4+iq <<endl;
2162  // cout <<"fill disk " << iD <<" quad: " << iq << " with " << clusCounts[iD*4+iq] <<" clusters from tracks"<<endl;
2163  numTrackHits[iD*4+iq]->Fill(clusCounts[iD*4+iq]);
2164 
2165  }}
2166  // cout <<"we have " << m_tracks.size()-oldNumTracks<< " tracks in this event " <<endl;
2167  // cout <<"oldNumtracsk: " << oldNumTracks << " new size: " << m_tracks.size() <<endl;
2168  numTracks->Fill(m_tracks.size()-oldNumTracks);
2169 
2170  return ierr;
2171 
2172 };
2173 
2174 StFgtGenAVEMaker::StFgtGenAVEMaker( const Char_t* name): StFgtGeneralBase( name ),useChargeMatch(false),runningEvtNr(0),hitCounter(0),hitCounterR(0),printCounter(0),fitCounter(0)
2175 {
2176 
2177  // cout <<"AVE constructor!!" <<endl;
2178  int numTb=7;
2179  mPulseShapePtr=new TF1("pulseShape","[0]*(x>[4])*(x-[4])**[1]*exp(-[2]*(x-[4]))+[3]",0,numTb);
2180  mPulseShapePtr->SetParName( 0, "C" );
2181  mPulseShapePtr->SetParName( 1, "a" );
2182  mPulseShapePtr->SetParName( 2, "b" );
2183  mPulseShapePtr->SetParName( 3, "ped" );
2184  mPulseShapePtr->SetParName( 4, "t0" );
2185 
2186  mHistPtr=new TH1F("tempFitHist","tempFitHist",numTb,0,numTb);
2187  mHistPtr2=new TH1F("tempFitHist2","tempFitHist2",numTb,0,numTb);
2188  mHistPtr2->SetLineColor(kBlue);
2189  mHistPtr2->SetFillColor(kBlue);
2190  mHistPtr2->SetMarkerColor(kBlue);
2191  mCanvas=new TCanvas("tmpCnvs","tmpCnvs");
2192 };
2193 
2194 StFgtGenAVEMaker::~StFgtGenAVEMaker()
2195 {
2196 
2197  //delete histogram arrays
2198 };
2199 
2200 Int_t StFgtGenAVEMaker::Finish(){
2201  StFgtGeneralBase::Finish();
2202  cout <<" closing txt file " << endl;
2203  outTxtFile->close();
2204  gStyle->SetPalette(1);
2205  cout <<"AVE finish function " <<endl;
2206  Int_t ierr = kStOk;
2207 
2209  vector<AVTrack>::iterator it=m_tracks.begin();
2210  cout <<"we found " << m_tracks.size() <<" tracks" <<endl;
2211  int counter=0;
2212  int goodTracks=0;
2213  for(;it!=m_tracks.end();it++)
2214  {
2215  cout <<" looking at track " << counter <<endl;
2216  counter++;
2217  // cout <<"This track has parameters: ";
2218  cout <<" mx: " << it->mx <<" my: " << it->my <<" bx: " << it->ax << " by: " << it->ay << " chi2: " << it->chi2 <<endl;
2219  Double_t vertZ = ( -( it->mx*it->ax + it->my*it->ay )/(it->mx*it->mx+it->my*it->my));
2220 
2221  pair<double,double> dca=getDca(it);
2222  // cout <<"dca: " << dca.second <<" z vertex: " << vertZ <<" or " << dca.first <<endl;
2223  if(it->chi2<MAX_DIST_CHI && fabs(vertZ)< VERTEX_CUT )
2224  {
2225  goodTracks++;
2226  cout <<"looking at track with mx: " << it->mx <<" ax: " << it->ax <<" my: " << it->my <<" ay: " << it->ay <<endl;
2227  hIpZ->Fill(dca.first);
2228  hIp->Fill(dca.first,dca.second);
2229  hTrkZ->Fill(vertZ);
2230  hMx->Fill(it->mx);
2231  hMy->Fill(it->my);
2232  hBx->Fill(it->ax);
2233  hBy->Fill(it->ay);
2234  hChi2->Fill(it->chi2);
2235  tpcFgtZVertexCorr->Fill(dca.first,it->ipZEv);
2236  tpcFgtZVertexCorr2->Fill(vertZ,it->ipZEv);
2237  tpcFgtZVertexCorr3->Fill(vertZ,dca.first);
2238  hIpDca->Fill(dca.second);
2239  }
2240  }
2241  cout <<" we found " << goodTracks << " good Tracks " <<endl;
2244  TCanvas* cRadio=new TCanvas("radioPlots","radioPlot",1000,1500);
2245  TCanvas* cRadioLoose=new TCanvas("radioPlotsLoose","radioPlotLoose",1000,1500);
2246  TCanvas* cRadioR=new TCanvas("radioPlotsR","radioPlotR",1000,1500);
2247  TCanvas* cRadioPhi=new TCanvas("radioPlotsPhi","radioPlotPhi",1000,1500);
2248 
2249  TCanvas* cRadioHits=new TCanvas("radioPlotsHits","radioPlotHits",1000,1500);
2250  TCanvas* cRadioNonHits=new TCanvas("radioPlotsNonHits","radioPlotNonHits",1000,1500);
2251  // cout <<"divide "<<endl;
2252  cRadio->Divide(2,3); //6 discs
2253  cRadioR->Divide(2,3); //6 discs
2254  cRadioPhi->Divide(2,3); //6 discs
2255  cRadioLoose->Divide(2,3); //6 discs
2256  cRadioHits->Divide(2,3); //6 discs
2257  cRadioNonHits->Divide(2,3); //6 discs
2258  TCanvas* cRPRatio=new TCanvas("rPhiRatio","rPhiRatios",1000,1500);
2259  cRPRatio->Divide(2,3); //6 discs
2260 
2261  TCanvas* cREff=new TCanvas("crEff","crEff",1000,1500);
2262 
2263  cREff->Divide(2,3); //6 discs
2264  // cout <<"drawing hits " <<endl;
2265  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
2266  {
2267  // cRadio->cd(iD+1)->SetLogz();
2268  cRadioHits->cd(iD+1);
2269  radioPlotsEff[iD]->Draw("colz");
2270  cRadioNonHits->cd(iD+1);
2271  radioPlotsNonEff[iD]->Draw("colz");
2272  }
2273  char buffer[100];
2274 
2275 
2276  sprintf(buffer,"%s/clusterPics.root",fileBase);
2277  cout <<"setting cluster pic file to : " << buffer <<endl;
2278  TFile *fClu = new TFile(buffer,"recreate");
2279  fClu->cd();
2280 
2281  for(unsigned int i=0;i<v_hClusP.size();i++)
2282  {
2283  (v_hClusP[i])->Write();
2284  }
2285  for(unsigned int i=0;i<v_hClusR.size();i++)
2286  {
2287  (v_hClusR[i])->Write();
2288  }
2289  fClu->Write();
2290  fClu->Close();
2291  sprintf(buffer,"%s/signalShapes.root",fileBase);
2292  cout <<"setting signal shapes file to : " << buffer <<endl;
2293  TFile *f1 = new TFile(buffer,"recreate");
2294  f1->cd();
2295 
2296  //normalize
2297  for(int iB=1;iB<8;iB++)
2298  {
2299  exPulseMaxAdcNormP->SetBinContent(iB,exPulseMaxAdcNormP->GetBinContent(iB)/(float)pulseCounterP);
2300  exPulseSigP->SetBinContent(iB,exPulseSigP->GetBinContent(iB)/(float)pulseCounterP);
2301 
2302  exPulseMaxAdcNormR->SetBinContent(iB,exPulseMaxAdcNormR->GetBinContent(iB)/(float)pulseCounterR);
2303  exPulseSigR->SetBinContent(iB,exPulseSigR->GetBinContent(iB)/(float)pulseCounterR);
2304 
2305  exPulseMaxAdcNormTrackP->SetBinContent(iB,exPulseMaxAdcNormTrackP->GetBinContent(iB)/(float)pulseCounterTP);
2306  exPulseSigTrackP->SetBinContent(iB,exPulseSigTrackP->GetBinContent(iB)/(float)pulseCounterTP);
2307 
2308  exPulseMaxAdcNormTrackR->SetBinContent(iB,exPulseMaxAdcNormTrackR->GetBinContent(iB)/(float)pulseCounterTR);
2309  exPulseSigTrackR->SetBinContent(iB,exPulseSigTrackR->GetBinContent(iB)/(float)pulseCounterTR);
2310  }
2311 
2312 
2313  chargeCorrSum3->Write();
2314  chargeCorrMaxStrip->Write();
2315  chargeCorrMaxAdc->Write();
2316 
2317  exPulseMaxAdcNormP->Write();
2318  exPulseSigP->Write();
2319 
2320  exPulseMaxAdcNormR->Write();
2321  exPulseSigR->Write();
2322 
2323  exPulseMaxAdcNormTrackP->Write();
2324  exPulseSigTrackP->Write();
2325 
2326  exPulseMaxAdcNormTrackR->Write();
2327  exPulseSigTrackR->Write();
2328  numTracks->Write();
2329 
2330 
2331  for(int xx=0; xx<22; xx++){
2332 
2333  disk1QuadA[xx]->Write();
2334  }
2335 
2336 
2337 
2338  for(int iD=0;iD<kFgtNumDiscs;iD++)
2339  {
2340 
2341  //Added for writing geoId, R, phi cluster histrograms
2342 
2343  clusterGeoId[iD]->Write();
2344  clustersR[iD]->Write();
2345  clustersP[iD]->Write();
2346 
2347 
2348  for(int iq=0;iq<4;iq++)
2349  {
2350  numClustersR[iD*4+iq]->Write();
2351  numClustersPhi[iD*4+iq]->Write();
2352  numTrackHits[iD*4+iq]->Write();
2353 
2354 
2355  maxTbCloseClusterP[iD*4+iq]->Write();
2356  maxTbCloseClusterR[iD*4+iq]->Write();
2357  maxAdcCloseClusterP[iD*4+iq]->Write();
2358  maxSigTrackClusterP[iD*4+iq]->Write();
2359  numFSigCloseClusterP[iD*4+iq]->Write();
2360  numFirstHighCloseClusterP[iD*4+iq]->Write();
2361  maxAdcCloseClusterP[iD*4+iq]->Write();
2362  maxSigCloseClusterP[iD*4+iq]->Write();
2363  secondToLastRatioCloseClusterP[iD*4+iq]->Write();
2364  maxAdcCloseClusterR[iD*4+iq]->Write();
2365  maxSigCloseClusterR[iD*4+iq]->Write();
2366  numFSigCloseClusterR[iD*4+iq]->Write();
2367  numFirstHighCloseClusterR[iD*4+iq]->Write();
2368  maxAdcCloseClusterR[iD*4+iq]->Write();
2369  maxSigCloseClusterR[iD*4+iq]->Write();
2370  secondToLastRatioCloseClusterR[iD*4+iq]->Write();
2371  firstTbSigCloseClusterR[iD*4+iq]->Write();
2372  firstTbSigCloseClusterP[iD*4+iq]->Write();
2373 
2374  maxTbTrackClusterP[iD*4+iq]->Write();
2375  maxTbTrackClusterR[iD*4+iq]->Write();
2376  maxAdcTrackClusterP[iD*4+iq]->Write();
2377  maxSigTrackClusterP[iD*4+iq]->Write();
2378  numFSigTrackClusterP[iD*4+iq]->Write();
2379  numFirstHighTrackClusterP[iD*4+iq]->Write();
2380  maxAdcTrackClusterP[iD*4+iq]->Write();
2381  maxSigTrackClusterP[iD*4+iq]->Write();
2382  secondToLastRatioTrackClusterP[iD*4+iq]->Write();
2383  maxAdcTrackClusterR[iD*4+iq]->Write();
2384  maxSigTrackClusterR[iD*4+iq]->Write();
2385  numFSigTrackClusterR[iD*4+iq]->Write();
2386  numFirstHighTrackClusterR[iD*4+iq]->Write();
2387  maxAdcTrackClusterR[iD*4+iq]->Write();
2388  maxSigTrackClusterR[iD*4+iq]->Write();
2389  secondToLastRatioTrackClusterR[iD*4+iq]->Write();
2390  firstTbSigTrackClusterR[iD*4+iq]->Write();
2391  firstTbSigTrackClusterP[iD*4+iq]->Write();
2392 
2393  }
2394  }
2395 
2396  for(int iD=0;iD<6;iD++)
2397  {
2398  for(int binAPVi=0;binAPVi<40;binAPVi++)
2399  {
2400  APVfitChi2P[iD*40+binAPVi]->Write();
2401  APVfitChi2R[iD*40+binAPVi]->Write();
2402  APVfitAmpP[iD*40+binAPVi]->Write();
2403  APVfitAmpR[iD*40+binAPVi]->Write();
2404  APVfitT0P[iD*40+binAPVi]->Write();
2405  APVfitT0R[iD*40+binAPVi]->Write();
2406 
2407  APVfirstTbSigCloseClusterP[iD*40+binAPVi]->Write();
2408  APVmaxAdcCloseClusterP[iD*40+binAPVi]->Write();
2409  APVmaxTbCloseClusterP[iD*40+binAPVi]->Write();
2410  APVnumFSigCloseClusterP[iD*40+binAPVi]->Write();
2411  APVnumFirstHighCloseClusterP[iD*40+binAPVi]->Write();
2412  APVmaxSigCloseClusterP[iD*40+binAPVi]->Write();
2413  APVsecondToLastRatioCloseClusterP[iD*40+binAPVi]->Write();
2414  APVfirstTbSigCloseClusterR[iD*40+binAPVi]->Write();
2415  APVmaxAdcCloseClusterR[iD*40+binAPVi]->Write();
2416  APVmaxTbCloseClusterR[iD*40+binAPVi]->Write();
2417  APVnumFSigCloseClusterR[iD*40+binAPVi]->Write();
2418  APVnumFirstHighCloseClusterR[iD*40+binAPVi]->Write();
2419  APVmaxSigCloseClusterR[iD*40+binAPVi]->Write();
2420  APVsecondToLastRatioCloseClusterR[iD*40+binAPVi]->Write();
2421  }
2422  }
2423 
2424 
2425  cout <<"writen and closed " << endl;
2426 
2427  cout <<"drawing hits2 " <<endl;
2430 
2431  TCanvas* cClusterSizeR=new TCanvas("clusterSizeR","clusterSizeR",1000,1500);
2432  cClusterSizeR->Divide(2,3);
2433  TCanvas* cClusterSizePhi=new TCanvas("clusterSizePhi","clusterSizePhi",1000,1500);
2434  cClusterSizePhi->Divide(2,3);
2435  TCanvas* cChargeCorr=new TCanvas("chargeCorr","chargeCorr",1000,1500);
2436  cChargeCorr->Divide(3,4);
2437 
2438  TCanvas* cClusterChargePhi=new TCanvas("clusterChargePhi","clusterChargePhi",1000,1500);
2439  cClusterChargePhi->Divide(2,3);
2440  TCanvas* cClusterChargeR=new TCanvas("clusterChargeR","clusterChargeR",1000,1500);
2441  cClusterChargeR->Divide(2,3);
2442 
2443  TCanvas cIPProj;
2444  hIp->Draw("colz");
2445  hIp->Write();
2447 
2448  hBx->Draw();
2449  hBx->Write();
2451  hBy->Draw();
2452  hBy->Write();
2454  hMx->Draw();
2455  hMx->Write();
2457  hMy->Draw();
2458  hMy->Write();
2460 
2461 
2462  hIpZ->Draw();
2463  hIpZ->Write();
2465 
2466 
2467  hIpDca->Draw();
2468  hIpDca->Write();
2470 
2471  hTrkZ->Draw();
2472  hTrkZ->Write();
2474 
2475  hResidua->Draw();
2476  hResidua->Write();
2477  hResiduaX->Draw();
2478  hResiduaX->Write();
2479 
2480  hResiduaY->Draw();
2481  hResiduaY->Write();
2482  hResiduaR->Draw();
2483  hResiduaR->Write();
2484  hResiduaP->Draw();
2485  hResiduaP->Write();
2487 
2488  hChi2->Draw();
2489  hChi2->Write();
2491 
2492  tpcFgtZVertexCorr->Draw("colz");
2493  tpcFgtZVertexCorr->Write();
2495 
2496  tpcFgtZVertexCorr2->Draw("colz");
2497  tpcFgtZVertexCorr2->Write();
2499  tpcFgtZVertexCorr3->Draw("colz");
2500  tpcFgtZVertexCorr3->Write();
2502 
2503 
2504  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
2505  {
2506  cClusterSizeR->cd(iD+1);
2507  h_clusterSizeR[iD]->Draw();
2508  cClusterSizePhi->cd(iD+1);
2509  h_clusterSizePhi[iD]->Draw();
2510  for(int iq=0;iq<4;iq++)
2511  {
2512  cChargeCorr->cd(iD*4+iq+1);
2513  chargeCorr[iD*4+iq]->Draw("colz");
2514  chargeCorr[iD*4+iq]->Write();
2515  clusterSizeR[iD*4+iq]->Write();
2516  clusterSizeP[iD*4+iq]->Write();
2517  }
2518  cClusterChargeR->cd(iD+1);
2519  h_clusterChargeR[iD]->Draw();
2520  cClusterChargePhi->cd(iD+1);
2521  h_clusterChargePhi[iD]->Draw();
2522  h_clusterChargeR[iD]->Write();
2523  h_clusterChargePhi[iD]->Write();
2524 
2525  }
2526  f1->Write();
2527 
2533 
2534  cout <<"saving .." <<endl;
2535  doNormalize(radioPlotsEffR, radioPlotsNonEffR);
2536  doNormalize(radioPlotsEffPhi, radioPlotsNonEffPhi);
2537 
2538 
2539  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
2540  {
2541  sprintf(buffer,"allCountsLooseDisk_%d",iD+1);
2542  TH2D* tmpAllCountsLoose=(TH2D*)radioPlotsEffLoose[iD]->Clone(buffer);
2543  tmpAllCountsLoose->Write();
2544  sprintf(buffer,"allClusterCountsDisk_%d",iD+1);
2545  TH2D* tmpAllClusterCounts=(TH2D*)radioPlotsEff[iD]->Clone(buffer);
2546  tmpAllClusterCounts->Write();
2547  }
2548 
2549  doNormalize(radioPlotsEffLoose, radioPlotsNonEffLoose);
2550 
2551 
2552  for(Int_t iD=0;iD<kFgtNumDiscs;iD++)
2553  {
2554  cRadio->cd(iD+1);
2555  sprintf(buffer,"allCountsDisk_%d",iD+1);
2556  TH2D* tmpAllCounts=(TH2D*)radioPlotsEff[iD]->Clone(buffer);
2557  tmpAllCounts->Write();
2558  radioPlotsEff[iD]->Add(radioPlotsNonEff[iD]);//all counts
2559  // radioPlotsEff[iD]->Add(radioPlotsNonEff[iD],-1); //subtract non eff
2560 
2562  for(int nx=1;nx<radioPlotsEff[iD]->GetNbinsX()+1;nx++)
2563  {
2564  for(int ny=1;ny<radioPlotsEff[iD]->GetNbinsY()+1;ny++)
2565  {
2566  Double_t denom=radioPlotsEff[iD]->GetBinContent(nx,ny);
2567  if(denom>0 && (tmpAllCounts->GetBinContent(nx,ny)/denom)<=1.0)
2568  {
2569  // cout <<" efficiency bin nx: " << nx << " ny: " << ny << ", num counts " << tmpAllCounts->GetBinContent(nx,ny) << " / " << denom << " : " << tmpAllCounts->GetBinContent(nx,ny)/denom <<endl;
2570  radioPlotsEff[iD]->SetBinContent(nx,ny,tmpAllCounts->GetBinContent(nx,ny)/denom);
2571  if(iD==m_effDisk)
2572  {
2573  // cout <<"chargeRatio is: " << chargeRatioInEffDisk->GetBinContent(nx,ny)<<endl;
2574  chargeRatioInEffDisk->SetBinContent(nx,ny,chargeRatioInEffDisk->GetBinContent(nx,ny)/denom);
2575  chargeAsymInEffDisk->SetBinContent(nx,ny,chargeAsymInEffDisk->GetBinContent(nx,ny)/denom);
2576  }
2577  }
2578  else
2579  {
2580  // cout <<" efficiency bin nx: " << nx << " ny: " << ny << ", num counts " << tmpAllCounts->GetBinContent(nx,ny) << " / " << denom << " : 0.0" <<endl;
2581  radioPlotsEff[iD]->SetBinContent(nx,ny,0.0);
2582  }
2583  }
2584  }
2585  // radioPlotsEff[iD]->Divide(tmpAllCounts);
2586  radioPlotsEff[iD]->SetMaximum(1.0);
2587  radioPlotsEff[iD]->Draw("colz");
2588 
2589  radioPlotsEffR[iD]->SetMaximum(1.0);
2590  radioPlotsEffPhi[iD]->SetMaximum(1.0);
2591  radioPlotsEffLoose[iD]->SetMaximum(1.0);
2592  cRadioR->cd(iD+1);
2593  radioPlotsEffR[iD]->Draw("colz");
2594  radioPlotsEffR[iD]->Write();
2595  cRadioPhi->cd(iD+1);
2596  radioPlotsEffPhi[iD]->Draw("colz");
2597  radioPlotsEffPhi[iD]->Write();
2598  cRadioLoose->cd(iD+1);
2599  radioPlotsEffLoose[iD]->Draw("colz");
2600  radioPlotsEffLoose[iD]->Write();
2601  radioPlotsNonEffLoose[iD]->Write();
2602  radioPlotsNonEff[iD]->Write();
2603  cRPRatio->cd(iD+1);
2604  rPhiRatioPlots[iD]->Draw();
2605  cREff->cd(iD+1);
2606  f1->Write();
2607 
2608  TH1D* tmpR=(TH1D*)rEff[iD]->Clone("tmpR");
2609  rEff[iD]->Add(rNonEff[iD]);
2610  for(int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
2611  {
2612  Double_t denom=rEff[iD]->GetBinContent(nx);
2613  if(denom>0)
2614  rEff[iD]->SetBinContent(nx,tmpR->GetBinContent(nx)/denom);
2615  else
2616  rEff[iD]->SetBinContent(nx,0.0);
2617  }
2618  rEff[iD]->Draw();
2619  }
2620 
2623 
2627 
2628  cRadio->cd(0);
2629  chargeRatioInEffDisk->Draw("colz");
2631  chargeAsymInEffDisk->Draw("colz");
2633  chargeCorrInEffDisk->Draw("colz");
2634  chargeCorrInEffDisk->Write();
2636  hChargeAsym->Draw();
2638  hChargeRatio->Draw();
2640 
2641 
2644 
2647  f1->Write();
2649  f1->Close();
2650  myRootFile->Write();
2651  myRootFile->Close();
2652 
2653  pulsePictureFile->Write();
2654  cout <<"returning after finish" <<endl;
2655  return ierr;
2656 };
2657 
2658 
2659 
2660 // construct histograms
2661 
2662 
2663 Int_t StFgtGenAVEMaker::Init(){
2664  outTxtFile=new ofstream;
2665  outTxtFile->open("clusExpectations.txt");
2666  cluNotFoundTxt=new ofstream;
2667  cluNotFoundTxt->open("clusNotFound.txt");
2668  cout <<"AVE!!" <<endl;
2669  myRootFile=new TFile("clusterEff.root","RECREATE");
2670  pulsePictureFile=new TFile("pulsePics.root","RECREATE");
2671 
2672  char buffer[100];
2673  for(int i=0;i<6;i++){
2674  for(int iq=0;iq<4;iq++){
2675  sprintf(buffer,"d%d_quad%d",i,iq);
2676  pulsePictureFile->mkdir(buffer);
2677  pulsePictureFile->cd(buffer);
2678  for(int iA=0;iA<24;iA++)
2679  {
2680  sprintf(buffer,"apv%d_R",iA);
2681  gDirectory->mkdir(buffer);
2682  sprintf(buffer,"apv%d_P",iA);
2683  gDirectory->mkdir(buffer);
2684  }
2685  pulsePictureFile->cd();
2686 
2687  }}
2688  myRootFile->cd();
2689  // outTxtFile=new ofstream;
2690  // outTxtFile->open("clusters.txt");
2691 
2692 
2693  Int_t ierr = kStOk;
2694 
2695  chargeRatioInEffDisk=new TH2D("chargeRatioInEffDisk","chargeRatioInEffDisk",NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2696  chargeRatioInEffDisk->SetMaximum(2.0);
2697  chargeAsymInEffDisk=new TH2D("chargeAsymInEffDisk","chargeAsymInEffDisk",NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2698  chargeAsymInEffDisk->SetMaximum(1.0);
2699  chargeCorrInEffDisk=new TH2D("chargeCorrInEffDisk","chargeCorrInEffDisk",500,0,50000,500,0,50000);
2700  hChargeAsym=new TH1D("chargeAsym","chargeAsym",100,0,50);
2701  hChargeRatio=new TH1D("chargeRatio","chargeRatio",100,0,50);
2702 
2703 
2704  chargeCorrSum3=new TH2D("chargeCorrSum3","chargeCorrSum3",500,0,50000,500,0,50000);
2705  chargeCorrMaxStrip=new TH2D("chargeCorrMaxStrip","chargeCorrMaxStrip",500,0,20000,500,0,20000);
2706  chargeCorrMaxAdc=new TH2D("chargeCorrMaxAdc","chargeCorrMaxAdc",500,0,5000,500,0,5000);
2707 
2708 
2709  radioPlotsEff=new TH2D*[kFgtNumDiscs];
2710  radioPlotsNonEff=new TH2D*[kFgtNumDiscs];
2711  radioPlotsEffR=new TH2D*[kFgtNumDiscs];
2712  radioPlotsNonEffR=new TH2D*[kFgtNumDiscs];
2713  radioPlotsEffPhi=new TH2D*[kFgtNumDiscs];
2714  radioPlotsEffLoose=new TH2D*[kFgtNumDiscs];
2715  radioPlotsNonEffPhi=new TH2D*[kFgtNumDiscs];
2716  radioPlotsNonEffLoose=new TH2D*[kFgtNumDiscs];
2717  rPhiRatioPlots=new TH1D*[kFgtNumDiscs];
2718 
2719 
2720  exPulseMaxAdcNormP=new TH1F("pulseMaxAdcNormP","pulseMaxAdcNormP",7,0,6);
2721  exPulseSigP=new TH1F("pulseSigNormP","pulseSigNormP",7,0,6);
2722 
2723  exPulseMaxAdcNormR=new TH1F("pulseMaxAdcNormR","pulseMaxAdcNormR",7,0,6);
2724  exPulseSigR=new TH1F("pulseSigNormR","pulseSigNormR",7,0,6);
2725 
2726  exPulseMaxAdcNormTrackP=new TH1F("pulseMaxAdcNormTrackP","pulseMaxAdcNormTrackP",7,0,6);
2727  exPulseSigTrackP=new TH1F("pulseSigNormTrackP","pulseSigNormTrackP",7,0,6);
2728 
2729  exPulseMaxAdcNormTrackR=new TH1F("pulseMaxAdcNormTrackR","pulseMaxAdcNormTrackR",7,0,6);
2730  exPulseSigTrackR=new TH1F("pulseSigNormTrackR","pulseSigNormTrackR",7,0,6);
2731 
2732 
2733  pulseCounterP=0;
2734  pulseCounterR=0;
2735 
2736  pulseCounterTP=0;
2737  pulseCounterTR=0;
2738 
2739  createPlots(&numClustersR,kFgtNumDiscs*4,"numClustersR",101,0,100);
2740  createPlots(&numClustersPhi,kFgtNumDiscs*4,"numClustersPhi",101,0,100);
2741  createPlots(&numTrackHits,kFgtNumDiscs*4,"numTrackHits",101,0,100);
2742  numTracks=new TH1I("numTracksPerEvent","numTracksPerEvent",101,0,100);
2743 
2744  createPlots(&firstTbSigCloseClusterR,kFgtNumDiscs*4,"firstTbSigCloseClusterR",100,0,20);
2745  createPlots(&firstTbSigCloseClusterP,kFgtNumDiscs*4,"firstTbSigCloseClusterP",100,0,20);
2746  createPlots(&firstTbSigTrackClusterR,kFgtNumDiscs*4,"firstTbSigTrackClusterR",100,0,20);
2747  createPlots(&firstTbSigTrackClusterP,kFgtNumDiscs*4,"firstTbSigTrackClusterP",100,0,20);
2748 
2749 
2750  createPlots(&maxAdcTrackClusterR,kFgtNumDiscs*4,"maxAdcTrackClusterR",100,0,5000);
2751  createPlots(&maxAdcCloseClusterR,kFgtNumDiscs*4,"maxAdcCloseClusterR",100,0,5000);
2752  createPlots(&maxSigTrackClusterR,kFgtNumDiscs*4,"maxSigTrackClusterR",100,1,200);
2753  createPlots(&maxSigCloseClusterR,kFgtNumDiscs*4,"maxSigCloseClusterR",100,1,200);
2754  createPlots(&numFSigTrackClusterR,kFgtNumDiscs*4,"numFSigTrackClusterR",9,1,8);
2755  createPlots(&maxTbCloseClusterR,kFgtNumDiscs*4,"maxTbCloseClusterR",8,0,7);
2756 
2757  createPlots(&maxTbCloseClusterP,kFgtNumDiscs*4,"maxTbCloseClusterP",8,0,7);
2758 
2759  createPlots(&maxTbTrackClusterR,kFgtNumDiscs*4,"maxTbTrackClusterR",8,0,7);
2760 
2761  createPlots(&maxTbTrackClusterP,kFgtNumDiscs*4,"maxTbTrackClusterP",9,0,8);
2762  createPlots(&numFSigCloseClusterR,kFgtNumDiscs*4,"numFSigCloseClusterR",9,0,8);
2763 
2764  createPlots(&numFirstHighTrackClusterR,kFgtNumDiscs*4,"numFirstHighTrackClusterR",8,0,7);
2765  createPlots(&numFirstHighCloseClusterR,kFgtNumDiscs*4,"numFirstHighCloseClusterR",8,0,7);
2766 
2767  createPlots(&maxAdcTrackClusterP,kFgtNumDiscs*4,"maxAdcTrackClusterP",100,0,5000);
2768  createPlots(&maxAdcCloseClusterP,kFgtNumDiscs*4,"maxAdcCloseClusterP",100,0,5000);
2769  createPlots(&maxSigTrackClusterP,kFgtNumDiscs*4,"maxSigTrackClusterP",100,1,200);
2770  createPlots(&maxSigCloseClusterP,kFgtNumDiscs*4,"maxSigCloseClusterP",100,1,200);
2771  createPlots(&numFSigTrackClusterP,kFgtNumDiscs*4,"numFSigTrackClusterP",9,0,8);
2772  createPlots(&numFSigCloseClusterP,kFgtNumDiscs*4,"numFSigCloseClusterP",9,0,8);
2773  createPlots(&numFirstHighTrackClusterP,kFgtNumDiscs*4,"numFirstHighTrackClusterP",8,0,7);
2774  createPlots(&numFirstHighCloseClusterP,kFgtNumDiscs*4,"numFirstHighCloseClusterP",8,0,7);
2775  createPlots(&secondToLastRatioCloseClusterP,kFgtNumDiscs*4,"secondToLastRatioClosClusterP",100,0,5);
2776 
2777  createPlots(&secondToLastRatioCloseClusterR,kFgtNumDiscs*4,"secondToLastRatioClosClusterR",100,0,5);
2778 
2779  createPlots(&secondToLastRatioTrackClusterP,kFgtNumDiscs*4,"secondToLastRatioTrackClusterP",100,0,5);
2780 
2781  createPlots(&secondToLastRatioTrackClusterR,kFgtNumDiscs*4,"secondToLastRatioTrackClusterR",100,0,5);
2782 
2783 
2784 
2785  createPlots(&APVfitChi2P,kFgtNumDiscs*40,"APVfitChi2P",100,0,30);
2786  createPlots(&APVfitChi2R,kFgtNumDiscs*40,"APVfitChi2R",100,0,30);
2787  createPlots(&APVfitAmpP,kFgtNumDiscs*40,"APVfitAmpP",100,0,30);
2788  createPlots(&APVfitAmpR,kFgtNumDiscs*40,"APVfitAmpR",100,0,30);
2789  createPlots(&APVfitT0P,kFgtNumDiscs*40,"APVfitT0P",100,0,30);
2790  createPlots(&APVfitT0R,kFgtNumDiscs*40,"APVfitT0R",100,0,30);
2791 
2792  createPlots(&APVfirstTbSigCloseClusterP,kFgtNumDiscs*40,"APVfirstTbSigCloseClusterP",100,0,20);
2793  createPlots(&APVfirstTbSigCloseClusterR,kFgtNumDiscs*40,"APVfirstTbSigCloseClusterR",100,0,20);
2794  createPlots(&APVmaxAdcCloseClusterP,kFgtNumDiscs*40,"APVmaxAdcCloseClusterP",100,0,5000);
2795  createPlots(&APVmaxAdcCloseClusterR,kFgtNumDiscs*40,"APVmaxAdcCloseClusterR",100,0,5000);
2796  createPlots(&APVmaxTbCloseClusterP,kFgtNumDiscs*40,"APVmaxTbCloseClusterP",8,0,7);
2797  createPlots(&APVmaxTbCloseClusterR,kFgtNumDiscs*40,"APVmaxTbCloseClusterR",8,0,7);
2798  createPlots(&APVnumFSigCloseClusterP,kFgtNumDiscs*40,"APVnumFSigCloseClusterP",9,0,8);
2799  createPlots(&APVnumFSigCloseClusterR,kFgtNumDiscs*40,"APVnumFSigCloseClusterR",9,0,8);
2800  createPlots(&APVnumFirstHighCloseClusterP,kFgtNumDiscs*40,"APVnumFirstHighCloseClusterP",8,0,7);
2801  createPlots(&APVnumFirstHighCloseClusterR,kFgtNumDiscs*40,"APVnumFirstHighCloseClusterR",8,0,7);
2802  createPlots(&APVmaxSigCloseClusterP,kFgtNumDiscs*40,"APVmaxSigCloseClusterP",100,1,200);
2803  createPlots(&APVmaxSigCloseClusterR,kFgtNumDiscs*40,"APVmaxSigCloseClusterR",100,1,200);
2804  createPlots(&APVsecondToLastRatioCloseClusterP,kFgtNumDiscs*40,"APVsecondToLastRatioCloseClusterP",100,0,5);
2805  createPlots(&APVsecondToLastRatioCloseClusterR,kFgtNumDiscs*40,"APVsecondToLastRatioCloseClusterR",100,0,5);
2806 
2807 
2808  //Joes plots
2809  createPlots(&clusterGeoId, kFgtNumDiscs,"clusterGeoId",32000,0,32000);
2810  createPlots(&clustersR, kFgtNumDiscs,"clustersR", 500, 0, 50);
2811  createPlots(&clustersP, kFgtNumDiscs,"clustersP", 100, -3.14159, 3.14159);
2812  //createPlots(&disk1QuadA, 22,"disk1QuadA",128,0,128);
2813  for (Int_t iii=0;iii< 22;iii++)
2814  {
2815  char buffer[100];
2816  sprintf(buffer, "%s_APV%d", "disk1QuadA",iii);
2817  disk1QuadA[iii]=new TH1I(buffer,buffer,128,0,128);
2818  }
2819 
2820  rEff=new TH1D*[kFgtNumDiscs];
2821  rNonEff=new TH1D*[kFgtNumDiscs];
2822  // cout <<"ave2" << endl;
2823 
2824  clusterSizeP=new TH1D*[kFgtNumDiscs*4];
2825  clusterSizeR=new TH1D*[kFgtNumDiscs*4];
2826  chargeCorr=new TH2D*[kFgtNumDiscs*4];
2827  h_clusterSizeR=new TH1D*[kFgtNumDiscs];
2828  h_clusterSizePhi=new TH1D*[kFgtNumDiscs];
2829 
2830  h_clusterChargeR=new TH1D*[kFgtNumDiscs];
2831  h_clusterChargePhi=new TH1D*[kFgtNumDiscs];
2832 
2833 
2834  hIp=new TH2D("Proj_to_IP","Proj_to_Ip",50,-100,100,50,-100,100);
2835  hBx=new TH1D("hBx","hBx",50,-100,100);
2836  hBy=new TH1D("hBy","hBy",50,-100,100);
2837  hMx=new TH1D("hMx","hMx",50,-100,100);
2838  hMy=new TH1D("hMy","My",50,-0.1,0.1);
2839  hIpZ=new TH1D("IP_Z","IP_Z",50,-100,100);
2840 
2841  hIpDca=new TH1D("ipDCA","ipDCA",50,-100,100);
2842  hTrkZ=new TH1D("z_Vtx_From_trk_fit","z_Vtx_From_trk_fit",50,-100,100);
2843  hResidua=new TH1D("residua","residua",100,0,2);
2844  hResiduaX=new TH2D("residuaX","residuaX",100,-10,10,200,0,1.0);
2845  hResiduaY=new TH2D("residuaY","residuaY",100,-40,-20,200,0,1.0);
2846  hResiduaR=new TH2D("residuaR","residuaR",100,12,32,200,0,1.0);
2847  hResiduaP=new TH2D("residuaP","residuaP",100,0,0.8,200,0,1.0);
2848 
2849 
2850  hChi2=new TH1D("chi2","chi2",50,0,2);
2851  tpcFgtZVertexCorr=new TH2D("tpc_fgt_corr","tpc_fgt_corr",100,-120,120,100,-120,120);
2852  tpcFgtZVertexCorr2=new TH2D("tpc_fgt_corr2","tpc_fgt_corr2",100,-120,120,100,-120,120);
2853  tpcFgtZVertexCorr3=new TH2D("fgt_fgt_corr","fgt_fgt_corr",50,-50,50,50,-50,50);
2854 
2855  for(int iD=0;iD<kFgtNumDiscs;iD++)
2856  {
2857 
2858  sprintf(buffer,"radioDiskEff_%d",iD);
2859  radioPlotsEff[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2860 
2861  sprintf(buffer,"radioDiskEffR_%d",iD);
2862  radioPlotsEffR[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2863 
2864  sprintf(buffer,"radioDiskEffPhi_%d",iD);
2865  radioPlotsEffPhi[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2866 
2867  sprintf(buffer,"radioDiskEffLoose_%d",iD);
2868  radioPlotsEffLoose[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2869  // cout <<"1" <<endl;
2870  sprintf(buffer,"rEff_%d",iD);
2871  rEff[iD]=new TH1D(buffer,buffer,100,0,DISK_DIM);
2872  sprintf(buffer,"rNonEff_%d",iD);
2873  rNonEff[iD]=new TH1D(buffer,buffer,100,0,DISK_DIM);
2874  sprintf(buffer,"clusterSizeR_Disk_%d",iD);
2875  h_clusterSizeR[iD]=new TH1D(buffer,buffer,20,0,20);
2876  h_clusterSizeR[iD]->SetFillColor(kYellow);
2877  sprintf(buffer,"clusterSizePhi_Disk_%d",iD);
2878  h_clusterSizePhi[iD]=new TH1D(buffer,buffer,20,0,20);
2879  h_clusterSizePhi[iD]->SetFillColor(kYellow);
2880  sprintf(buffer,"clusterChargeR_Disk_%d",iD);
2881  h_clusterChargeR[iD]=new TH1D(buffer,buffer,100,0,5000);
2882  h_clusterChargeR[iD]->SetFillColor(kYellow);
2883 
2884  sprintf(buffer,"clusterChargePhi_Disk_%d",iD);
2885  h_clusterChargePhi[iD]=new TH1D(buffer,buffer,100,0,5000);
2886  h_clusterChargePhi[iD]->SetFillColor(kYellow);
2887  // cout <<"2" <<endl;
2888  for(int nx=0;nx<radioPlotsEff[iD]->GetNbinsX();nx++)
2889  {
2890  for(int ny=0;ny<radioPlotsEff[iD]->GetNbinsY();ny++)
2891  {
2892  // radioPlotsEff[iD]->SetBinContent(nx,ny,0.1);//so that there is no divide by zero
2893  }
2894  }
2895  // cout <<"3" <<endl;
2896  for(int iq=0;iq<4;iq++)
2897  {
2898  sprintf(buffer,"r_phi_ChargeCorrelationInDisk_%d_quad_%d",iD+1,iq);
2899  chargeCorr[iD*4+iq]=new TH2D(buffer,buffer,200,0,70000,200,0,70000);
2900  sprintf(buffer,"clusterSizeRInDisk_%d_quad_%d",iD+1,iq);
2901  clusterSizeR[iD*4+iq]=new TH1D(buffer,buffer,100,0,100);
2902  sprintf(buffer,"clusterSizePInDisk_%d_quad_%d",iD+1,iq);
2903  clusterSizeP[iD*4+iq]=new TH1D(buffer,buffer,100,0,100);
2904  }
2905  sprintf(buffer,"radioDiskNonEff_%d",iD);
2906  radioPlotsNonEff[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2907  // cout <<" created non eff histo " << endl;
2908  sprintf(buffer,"radioDiskNonEffR_%d",iD);
2909  radioPlotsNonEffR[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2910  sprintf(buffer,"radioDiskNonEffPhi_%d",iD);
2911  radioPlotsNonEffPhi[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2912  sprintf(buffer,"radioDiskNonEffLoose_%d",iD);
2913  radioPlotsNonEffLoose[iD]=new TH2D(buffer,buffer,NUM_EFF_BIN,-DISK_DIM,DISK_DIM,NUM_EFF_BIN,-DISK_DIM,DISK_DIM);
2914  for(int nx=0;nx<rEff[iD]->GetNbinsX();nx++)
2915  {
2916  // rEff[iD]->SetBinContent(nx,0.1);
2917  }
2918  sprintf(buffer,"rPhiRatio_%d",iD);
2919  rPhiRatioPlots[iD]=new TH1D(buffer,buffer,100,-2,10);
2920  }
2921 
2922  return ierr;
2923 };
2924 ClassImp(StFgtGenAVEMaker);
StFgtGenAVEMaker::getQuadFromCoo
Short_t getQuadFromCoo(Double_t x, Double_t y)
this is too naive..., assumes non-rotated quads
Definition: StFgtGenAVEMaker.cxx:921
StFgtGenAVEMaker::findClosestPoint
Double_t findClosestPoint(float mx, float bx, float my, float by, double xE, double yE, Int_t iD)
Definition: StFgtGenAVEMaker.cxx:860
TDataSet::Clone
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone
Definition: TDataSet.cxx:308
TDataSet::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Definition: TDataSet.cxx:893
StFgtGenAVEMaker::fillStripHistos
void fillStripHistos(Float_t r, Float_t phi, Int_t iD, Int_t iq)
Definition: StFgtGenAVEMaker.cxx:315
StFgtGenAVEMaker::getTrack
Bool_t getTrack(vector< AVPoint > &points, Double_t ipZ)
Definition: StFgtGenAVEMaker.cxx:1174
StFgtGenAVEMaker::getDca
pair< double, double > getDca(vector< AVTrack >::iterator it)
Definition: StFgtGenAVEMaker.cxx:258
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5