StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFmsOfflineQaMaker.cxx
1 // \class StFmsOfflineQaMaker
2 // \author Akio Ogawa
3 //
4 // $Id: StFmsOfflineQaMaker.cxx,v 1.2 2016/06/08 19:55:11 akio Exp $
5 // $Log: StFmsOfflineQaMaker.cxx,v $
6 // Revision 1.2 2016/06/08 19:55:11 akio
7 // applying coverity report
8 //
9 // Revision 1.1 2016/01/26 19:54:33 akio
10 // Separated from StFmsFpsMaker... This is for FMS offline QA and also FMS-FPS alignments
11 //
12 //
13 
14 #include "StFmsOfflineQaMaker.h"
15 
16 #include "StMessMgr.h"
17 #include "Stypes.h"
18 
19 #include "StFmsDbMaker/StFmsDbMaker.h"
20 #include "StEnumerations.h"
21 #include "StEventTypes.h"
22 #include "StEvent/StEvent.h"
23 #include "StEvent/StFmsCollection.h"
24 #include "StEvent/StFmsHit.h"
25 #include "StEvent/StFmsPoint.h"
26 #include "StEvent/StFmsPointPair.h"
27 #include "StMuDSTMaker/COMMON/StMuTypes.hh"
28 
29 #include "TFile.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 
33 static const float EThr=5.0;
34 static const float ETotThr1=20.0;
35 
36 inline float project(float x, float z, float zp, float vz){
37  return x/(z-vz)*(zp-vz);
38 }
39 
40 static const int NTRG=12; //123=FMS-sm-bs123,456=FMS-lg-bs123,7=FMS-DiBS,8910=FMS-JP012,11=FMS-DiJP,13=LED
41 static const char *CTRG[NTRG] = {"SmBS1","SmBS2","SmBS3",
42  "LgBS1","LgBS2","LgBS3","DiBS",
43  "JP2","JP1","JP0","DiJP","LED"};
44 
45 int getFmsTrigId(const StTriggerId& trgid, int print=1){
46  const int TIDBASE=480800;
47  const int NBEAM=4;
48  const int MAXVERSION=3;
49  int trig=0;
50  // LOG_INFO << trgid <<endm;
51  LOG_INFO << "TRGID = ";
52  for(int k=0; k<NBEAM; k++){
53  for(int j=0; j<MAXVERSION; j++){
54  for(int i=1; i<=NTRG; i++){
55  int l=i;
56  if(i==12) l=13;
57  int id=TIDBASE + 10000*k + 20*j + l;
58  if(trgid.isTrigger(id)){
59  trig |= (1<<(i-1));
60  LOG_INFO << CTRG[i-1] << " ";
61  }
62  }
63  }
64  }
65  //include none-production id for now.... could be dangerous
66  static int nonProdId[NTRG]={42,43,44,45,46,47,48,49,50,51,52,53};
67  for(int i=0; i<NTRG; i++){
68  if(trgid.isTrigger(nonProdId[i])){
69  trig |= (1<<i);
70  LOG_INFO << CTRG[i] << " ";
71  }
72  }
73  LOG_INFO << endm;
74  return trig;
75 }
76 
77 ClassImp(StFmsOfflineQaMaker);
78 
79 StFmsOfflineQaMaker::StFmsOfflineQaMaker(const Char_t* name):
80  StMaker(name),mFilename((char *)"fmsqa.root")
81 {}
82 
83 StFmsOfflineQaMaker::~StFmsOfflineQaMaker(){}
84 
85 Int_t StFmsOfflineQaMaker::Init(){
86  mFmsDbMaker=static_cast<StFmsDbMaker*>(GetMaker("fmsDb"));
87  if(!mFmsDbMaker){
88  LOG_ERROR << "StFmsOfflineQaMaker::InitRun Failed to get StFmsDbMaker" << endm;
89  return kStFatal;
90  }
91 
92  char c[100];
93  mFile=new TFile(mFilename,"RECREATE");
94  static const int pm[kFpsNQuad][kFpsNLayer]={{1,1,1},{1,-1,-1},{-1,1,1},{-1,-1,-1}}; //+or- measured coordnates
95 
96  mERatio[0]= new TH1F("ERatioCluster-Hit","ERatioCH",50.0,0.0,1.2);
97  mERatio[1]= new TH1F("ERatioPoint-Cluster","ERatioPC",50.0,0.0,1.2);
98  mBC = new TH1F("BC","BC",120,0.0,120.0);
99  mTrig[0] = new TH1F("Trig","Trig", 13,0.0,13.0);
100  mTrig[1] = new TH1F("TrigAg","TrigAg",13,0.0,13.0);
101 
102  mFmsAdc = new TH1F("FmsAdc","FmsAdc",4096,0.0,4096.0);
103  mFmsHitLarge[0] = new TH2F("FmsHitLarge","FmsHitLarge",34,-17.0,17.0,34,0.0,34.0);
104  mFmsHitSmall[0] = new TH2F("FmsHitSmall","FmsHitSmall",24,-12.0,12.0,24,0.0,24.0);
105  mFmsHitLarge[1] = new TH2F("FmsHitLargeAg","FmsHitLargeAg",34,-17.0,17.0,34,0.0,34.0);
106  mFmsHitSmall[1] = new TH2F("FmsHitSmallAg","FmsHitSmallAg",24,-12.0,12.0,24,0.0,24.0);
107  mFpsMip[0] = new TH1F("FpsMipL1","FpsMipL1",100,0.0,5.0);
108  mFpsMip[1] = new TH1F("FpsMipL2","FpsMipL2",100,0.0,5.0);
109  mFpsMip[2] = new TH1F("FpsMipL3","FpsMipL3",180,0.0,90.0);
110 
111  mNTow[0] = new TH1F("NTowL","NTowL",50,0.0,50.0);
112  mNTow[1] = new TH1F("NTowS","NTowS",50,0.0,50.0);
113  mNTowE[0] = new TH2F("NTowEL","NTowEL",50,0.0,50.0,100.0,0.0,100.0);
114  mNTowE[1] = new TH2F("NTowES","NTowES",50,0.0,50.0,100.0,0.0,100.0);
115  mSigmax[0]= new TH1F("SigMaxL","SigMaxL",100.0,0.0,4.0);
116  mSigmax[1]= new TH1F("SigMaxS","SigMaxS",100.0,0.0,4.0);
117  mSigmin[0]= new TH1F("SigMinL","SigMinL",100.0,0.0,4.0);
118  mSigmin[1]= new TH1F("SigMinS","SigMinS",100.0,0.0,4.0);
119  mSigmaxE[0]= new TH2F("SigMaxEL","SigMaxEL",100.0,0.0,100.0,100,0.0,2.0);
120  mSigmaxE[1]= new TH2F("SigMaxES","SigMaxES",100.0,0.0,100.0,100,0.0,2.0);
121  mChi2[0] = new TH1F("Chi2L","Chi2L",200.0,0.0,200.0);
122  mChi2[1] = new TH1F("Chi2S","Chi2S",200.0,0.0,200.0);
123  mCluXY[0][0]= new TH2F("CluXYL","CluXYL",100.0,-17.0,17.0,100,0.0,34.0);
124  mCluXY[1][0]= new TH2F("CluXYS","CluXYS",100.0,-12.0,12.0,100,0.0,24.0);
125  mCluXY[0][1]= new TH2F("CluXYLAg","CluXYLAg",100.0,-17.0,17.0,100,0.0,34.0);
126  mCluXY[1][1]= new TH2F("CluXYSAg","CluXYSAg",100.0,-12.0,12.0,100,0.0,24.0);
127  mCluEta[0]=new TH1F("CluEta0","CluEta0",100,2.5,4.5);
128  mCluEta[1]=new TH1F("CluEta1","CluEta1",100,2.5,4.5);
129  mCluEta[2]=new TH1F("CluEta2","CluEta2",100,2.5,4.5);
130  mCluEta[3]=new TH1F("CluEta3","CluEta3",100,2.5,4.5);
131  mCluEta[4]=new TH1F("CluEta4","CluEta4",100,2.5,4.5);
132  mCluEta[5]=new TH1F("CluEta5","CluEta5",100,2.5,4.5);
133 
134  for(int cut=0; cut<NCUT1; cut++){
135  for(int q=0; q<kFpsNQuad; q++){
136  for(int l=0; l<kFpsNLayer; l++){
137  float min,max;
138  if(pm[q][l]==1){min=-10.0; max=110;}
139  else {min=-110.0; max=10.0;}
140  sprintf(c,"FMSFPS_Q%1dL%1d_c%1d",q+1,l+1,cut);
141  mH2[q][l][cut]=new TH2F(c,c,120,min,max,21,0.5,21.5);
142  sprintf(c,"FMSFPSd_Q%1dL%1d_c%1d",q+1,l+1,cut);
143  mHd2[q][l][cut]=new TH2F(c,c,100,min,max,100,-50.0,50.0);
144  sprintf(c,"FMS-FPS_Q%1dL%1d_c%1d",q+1,l+1,cut);
145  mHd[q][l][cut]=new TH1F(c,c,100,-50.0,50.0);
146  }//loop over layer
147  }//loop over quad
148  sprintf(c,"NP_c%d",cut); mHn [cut]=new TH1F(c,c,100,0.0,100.0);
149  sprintf(c,"e_c%d",cut); mHene[cut]=new TH1F(c,c,100,0.0,100.0);
150  sprintf(c,"elo_c%d",cut); mHelo[cut]=new TH1F(c,c,100,0.0,4.0);
151  sprintf(c,"pt_c%d",cut); mHpt [cut]=new TH1F(c,c,100,0.0,10.0);
152  sprintf(c,"ept_c%d",cut); mHept[cut]=new TH2F(c,c,100,0.0,100.0,100,0.0,10.0);
153  sprintf(c,"eta_c%d",cut); mHeta[cut]=new TH1F(c,c,100,2.5,4.5);
154  sprintf(c,"phi_c%d",cut); mHphi[cut]=new TH1F(c,c,100,-3.2,3.2);
155  sprintf(c,"x_c%d",cut); mHx [cut]=new TH1F(c,c,100,-100.0,100.0);
156  sprintf(c,"y_c%d",cut); mHy [cut]=new TH1F(c,c,100,-100.0,100.0);
157  sprintf(c,"dxl_c%d",cut); mHdxL[cut]=new TH1F(c,c,120,-0.1,1.1);
158  sprintf(c,"dxs_c%d",cut); mHdxS[cut]=new TH1F(c,c,120,-0.1,1.1);
159  sprintf(c,"xy_c%d",cut); mHxy [cut]=new TH2F(c,c,200,-100.0,100.0,200,-100.0,100.0);
160  sprintf(c,"pid_c%d",cut); mHpid[cut]=new TH1F(c,c,41,-0.5,40.5);
161  sprintf(c,"pid2_c%d",cut); mHpid2[cut]=new TH1F(c,c,5,-0.5,4.5);
162  }//loop over cuts
163 
164  for(int cut=0; cut<NCUT2; cut++){
165  sprintf(c,"p_NP_c%d",cut); mPn [cut]=new TH1F(c,c, 16,0.0,16.0);
166  sprintf(c,"p_e_c%d",cut); mPene[cut]=new TH1F(c,c,100,0.0,100.0);
167  sprintf(c,"p_pt_c%d",cut); mPpt [cut]=new TH1F(c,c,100,0.0,10.0);
168  sprintf(c,"p_ept_c%d",cut); mPept[cut]=new TH2F(c,c,100,0.0,100.0,100,0.0,10.0);
169  sprintf(c,"p_eta_c%d",cut); mPeta[cut]=new TH1F(c,c,100,2.5,5.0);
170  sprintf(c,"p_phi_c%d",cut); mPphi[cut]=new TH1F(c,c,100,-3.2,3.2);
171  sprintf(c,"p_pid_c%d",cut); mPpid[cut]=new TH2F(c,c,5,-0.5,4.5,5,-0.5,4.5);
172  sprintf(c,"p_m1_c%d",cut); mPm1 [cut]=new TH1F(c,c,100,0.0,1.0);
173  sprintf(c,"p_m2_c%d",cut); mPm2 [cut]=new TH1F(c,c,100,0.0,5.0);
174  sprintf(c,"p_zgg_c%d",cut); mPzgg[cut]=new TH1F(c,c,100,0.0,1.0);
175  sprintf(c,"p_dgg_c%d",cut); mPdgg[cut]=new TH1F(c,c,200,0.0,20.0);
176  sprintf(c,"p_r30_c%d",cut); mPr30[cut]=new TH1F(c,c,100,0.0,3.0);
177  sprintf(c,"p_r100_c%d",cut); mPr100[cut]=new TH1F(c,c,100,0.0,2.0);
178  sprintf(c,"p_xy_c%d",cut); mPxy[cut]=new TH2F(c,c,200,-100.0,100.0,200,-100.0,100.0);
179  }
180  return kStOK;
181 }
182 
183 Int_t StFmsOfflineQaMaker::Finish(){
184  LOG_INFO << Form("Writing and closing %s",mFilename) << endm;
185  mFile->Write();
186  mFile->Close();
187  return kStOK;
188 }
189 
190 Int_t StFmsOfflineQaMaker::Make(){
191  mBunch=-1;
192  mTrigger=0;
193  StEvent* event = (StEvent*)GetInputDS("StEvent");
194  if(!event) {LOG_ERROR << "StFmsOfflineQaMaker::Make did not find StEvent"<<endm; return kStErr;}
195  if(event->triggerData()) mBunch = event->triggerData()->bunchId7Bit();
196  if(event->triggerIdCollection() && event->triggerIdCollection()->nominal()){
197  mTrigger = getFmsTrigId(*(event->triggerIdCollection()->nominal()));
198  }
199  StMuDst* mudst = (StMuDst*)GetInputDS("MuDst");
200  if(mudst){
201  LOG_INFO<<"StFmsFpsMaker::Make found Mudst, getting buunchId and triggerId from MuDst"<<endm;
202  mBunch = mudst->event()->triggerData()->bunchId7Bit();
203  mTrigger = getFmsTrigId(mudst->event()->triggerIdCollection().nominal());
204  }
205 
206  mFmsColl = event->fmsCollection();
207  if(!mFmsColl) {LOG_ERROR << "StFmsOfflineQaMaker::Make did not find StEvent->FmsCollection"<<endm; return kStErr;}
208 
209  if(mPrint) print();
210 
211  StSPtrVecFmsHit& hits = mFmsColl->hits();
212  StSPtrVecFmsCluster& clusters = mFmsColl->clusters();
213  StSPtrVecFmsPoint& points = mFmsColl->points();
214  StSPtrVecFpsSlat& slats = mFmsColl->fpsSlats();
215  vector<StFmsPointPair*>& pairs = mFmsColl->pointPairs();
216  int nh=mFmsColl->numberOfHits();
217  int nc=mFmsColl->numberOfClusters();
218  int np=mFmsColl->numberOfPoints();
219  int ns=slats.size();
220  int npair=mFmsColl->numberOfPointPairs();
221  float etot[3]={0.0,0.0,0.0};
222 
223  //Bunch counter and AG
224  int ag=0;
225  if (mBunch>=30 && mBunch< 40) {ag=1;}
226  else if(mBunch>=110 && mBunch<120) {ag=2;}
227  mBC->Fill(float(mBunch));
228  if(ag==0) mTrig[0]->Fill(0);
229  if(ag==1) mTrig[1]->Fill(0);
230  for(int i=0; i<NTRG; i++){
231  if(mTrigger & (1<<i)){
232  if(ag==0) mTrig[0]->Fill(float(i+1));
233  if(ag==1) mTrig[1]->Fill(float(i+1));
234  }
235  }
236 
237  //First fill some hists for hits for just QA
238  for(int i=0; i<nh; i++){
239  StFmsHit* hit=hits[i];
240  int det=hit->detectorId();
241  int ch=hit->channel();
242  if(det>=kFmsNorthLargeDetId && det<=kFmsSouthSmallDetId){
243  mFmsAdc->Fill(float(hit->adc()));
244  float row=mFmsDbMaker->getRowNumber(det,ch)-0.5;
245  float col=mFmsDbMaker->getColumnNumber(det,ch)-0.5;
246  if(mFmsDbMaker->northSouth(det)==0) col*=-1.0;
247  if(mFmsDbMaker->largeSmall(det)==0){
248  if(ag==0) mFmsHitLarge[0]->Fill(col,row,hit->energy());
249  if(ag==1) mFmsHitLarge[1]->Fill(col,row,hit->energy());
250  }else{
251  if(ag==0) mFmsHitSmall[0]->Fill(col,row,hit->energy());
252  if(ag==1) mFmsHitSmall[1]->Fill(col,row,hit->energy());
253  }
254  etot[0]+=hit->energy();
255  }else if(det==kFpsDetId){
256  int q,l,s;
257  mFmsDbMaker->fpsQLSfromSlatId(ch,&q,&l,&s);
258  if(l>=1 && l<=3) mFpsMip[l-1]->Fill(hit->energy());
259  }
260  }
261 
262  //Hists for clusters
263  for(int i=0; i<nc; i++){
264  StFmsCluster* clu=clusters[i];
265  etot[1]+=clu->energy();
266  if(clu->energy()<5.0) continue;
267  int det=clu->detectorId();
268  int ls=mFmsDbMaker->largeSmall(det);
269  if(ls<0) continue;
270  float nt=float(clu->nTowers());
271  mNTow[ls]->Fill(nt);
272  mNTowE[ls]->Fill(nt,clu->energy());
273  mSigmax[ls]->Fill(clu->sigmaMax());
274  mSigmin[ls]->Fill(clu->sigmaMin());
275  mSigmaxE[ls]->Fill(clu->energy(),clu->sigmaMax());
276  if(det%2==0){
277  if(ag==0) mCluXY[ls][0]->Fill(-clu->x(),clu->y());
278  if(ag==1) mCluXY[ls][1]->Fill(-clu->x(),clu->y());
279  }else{
280  if(ag==0) mCluXY[ls][0]->Fill(clu->x(),clu->y());
281  if(ag==1) mCluXY[ls][1]->Fill(clu->x(),clu->y());
282  }
283  int nphoton=clu->nPhotons();
284  float c2=0.0;
285  if(nphoton==0) {c2=0.0;}
286  else if(nphoton==1) {c2=clu->chi2Ndf1Photon();}
287  else if(nphoton==2) {c2=clu->chi2Ndf2Photon();}
288  else {LOG_INFO << Form("NPHOTON=%d ???",nphoton) << endm;}
289  mChi2[ls]->Fill(c2);
290  int ctgy=clu->category();
291  int type=(nphoton-1)*3 + ctgy;
292  //cout << Form("NP=%1d Ctgry=%1d Type=%d\n",nphoton,ctgy,type);
293  mCluEta[type]->Fill(clu->fourMomentum().pseudoRapidity());
294  }
295 
296  //Counting points and slats per quad
297  int npoint[kFpsNQuad], nfps[kFpsNQuad][kFpsNLayer];
298  memset(npoint,0,sizeof(npoint));
299  memset(nfps,0,sizeof(nfps));
300  //Get # of FMS point per quad
301  for(int i=0; i<np; i++) {
302  float x=points[i]->XYZ().x();
303  float y=points[i]->XYZ().y();
304  int q=0;
305  if (x>=0.0 && y>=0.0) {q=1;}
306  else if(x>=0.0 && y< 0.0) {q=2;}
307  else if(x< 0.0 && y>=0.0) {q=3;}
308  else if(x< 0.0 && y< 0.0) {q=4;}
309  npoint[q-1]++;
310  etot[2]+=points[i]->energy();
311  }
312  //Get # of FPS slats hitted per quad/layer
313  for(int i=0; i<ns; i++) {
314  if(slats[i]->mip()>0.5){
315  int slatid=slats[i]->slatId();
316  int q,l,s;
317  mFmsDbMaker->fpsQLSfromSlatId(slatid,&q,&l,&s);
318  nfps[q-1][l-1]++;
319  }
320  }
321  if(etot[0]>0.0) mERatio[0]->Fill(etot[1]/etot[0]);
322  if(etot[1]>0.0) mERatio[1]->Fill(etot[2]/etot[1]);
323 
324  //Loop over FMS points to set up cuts and counts
325  int cut1[NCUT1], n1[NCUT1];
326  memset(n1,0,sizeof(n1));
327  for(int i=0; i<np; i++) {
328  memset(cut1,0,sizeof(cut1));
329  float e=points[i]->energy();
330  float x=points[i]->XYZ().x();
331  float y=points[i]->XYZ().y();
332  float z=points[i]->XYZ().z();
333  int pid=points[i]->fpsPid();
334  int pid2=pid/10;
335  StLorentzVectorF v1=points[i]->fourMomentum();
336  int q=0;
337  if (x>=0.0 && y>=0.0) {q=1;}
338  else if(x>=0.0 && y< 0.0) {q=2;}
339  else if(x< 0.0 && y>=0.0) {q=3;}
340  else if(x< 0.0 && y< 0.0) {q=4;}
341 
342  cut1[0]=1; //all
343  int edgeType;
344  float distance=mFmsDbMaker->distanceFromEdge(points[i],edgeType);
345  /*
346  if(points[i]->detectorId()<=9) cut1[1]=1;
347  if(points[i]->detectorId()>=10) cut1[2]=1;
348  if(points[i]->detectorId()<=9 && edgeType==4) cut1[3]=1;
349  if(points[i]->detectorId()>=10 && edgeType==4) cut1[4]=1;
350  */
351  if(e>=EThr){
352  cut1[1]=1;
353  if(distance<-0.51 || edgeType==4){
354  cut1[2]=1;
355  if(distance<-0.51 || (mFmsColl->isMergeSmallToLarge() && edgeType==4)){
356  if(ag==0) {
357  cut1[3]=1;
358  if(pid2==1) {cut1[6]=1;}
359  else if(pid2==2) {cut1[7]=1;}
360  else if(pid2==3) {cut1[8]=1;}
361  else {cut1[9]=1;}
362  }
363  if(ag==1) cut1[5]=1;
364  }
365  }
366  if(edgeType==4 && distance>=-0.51) cut1[4]=1;
367  }
368  for(int c=0; c<NCUT1; c++){
369  if(cut1[c]==0) continue;
370  n1[c]++;
371  mHene[c]->Fill(e);
372  mHelo[c]->Fill(e);
373  mHpt[c]->Fill(v1.perp());
374  mHept[c]->Fill(e,v1.perp());
375  mHeta[c]->Fill(v1.pseudoRapidity());
376  mHphi[c]->Fill(v1.phi());
377  mHx[c]->Fill(x);
378  mHy[c]->Fill(y);
379  mHxy[c]->Fill(x,y);
380  int det=points[i]->detectorId();
381  float localx=points[i]->x()/mFmsDbMaker->getXWidth(det);
382  float localy=points[i]->y()/mFmsDbMaker->getYWidth(det);
383  float dlx=localx-int(localx);
384  float dly=localy-int(localy);
385  if(points[i]->detectorId()<=9){
386  mHdxL[c]->Fill(dlx);
387  mHdxL[c]->Fill(dly);
388  }else{
389  mHdxS[c]->Fill(dlx);
390  mHdxS[c]->Fill(dly);
391  }
392  mHpid[c]->Fill(float(pid));
393  mHpid2[c]->Fill(float(pid2));
394  }//loop over cuts
395 
396  //Now FPS
397  for(int l=1; l<=kFpsNLayer; l++) {
398  for(int s=1; s<=kFpsNSlat; s++) {
399  int slatid = mFmsDbMaker->fpsSlatId(q,l,s);
400  if(slatid<0) continue;
401  float xyz[3],width[3];
402  mFmsDbMaker->fpsPosition(q,l,s,xyz,width);
403  //printf("Q%1dL%1dS%02d hit=%6.3f xyz=%6.2f %6.2f %6.2f\n",q,l,s,mHit[q-1][l-1][s-1],xyz[0],xyz[1],xyz[2]);
404  if(slats[slatid]->mip()>0.5){
405  for(int c=0; c<NCUT1; c++){
406  if(cut1[c]==0) continue;
407  static const float zvertex=0.0;
408  if(l==1){
409  float xx=project(x,z,xyz[2],zvertex);
410  mH2[q-1][l-1][c]->Fill(xx,float(s));
411  mHd[q-1][l-1][c]->Fill(xx-xyz[0]);
412  mHd2[q-1][l-1][c]->Fill(xx,xx-xyz[0]);
413  }else{
414  float yy=project(y,z,xyz[2],zvertex);
415  mH2[q-1][l-1][c]->Fill(yy,float(s));
416  mHd[q-1][l-1][c]->Fill(yy-xyz[1]);
417  mHd2[q-1][l-1][c]->Fill(yy,y-xyz[1]);
418  }
419  }//loop over cuts
420  }//if FPS has mip
421  }//loop slat
422  }//loop layer
423  }//loop over points
424  for(int c=0; c<NCUT1; c++) mHn[c]->Fill(float(n1[c]));
425 
426  //Now pairs
427  int cut2[NCUT2], n2[NCUT2];
428  memset(n2,0,sizeof(n2));
429  for(int i=0; i<npair; i++) {
430  memset(cut2,0,sizeof(cut2));
431  StFmsPointPair* pair=pairs[i];
432  StFmsPoint* p0=pair->point(0);
433  StFmsPoint* p1=pair->point(1);
434  float e0=p0->energy();
435  float e1=p1->energy();
436  int pid=pair->fpsPid();
437  cut2[0]=1; //all
438  int edgeType0, edgeType1;
439  float d0=mFmsDbMaker->distanceFromEdge(p0,edgeType0);
440  float d1=mFmsDbMaker->distanceFromEdge(p1,edgeType1);
441  if(e0>=EThr && e1>=EThr && pair->energy()>ETotThr1){
442  cut2[1]=1;
443  if( (d0<-0.51 || edgeType0==4) &&
444  (d1<-0.51 || edgeType1==4) ){
445  cut2[2]=1;
446  if( (d0<-0.51 || (mFmsColl->isMergeSmallToLarge() && edgeType0==4)) &&
447  (d1<-0.51 || (mFmsColl->isMergeSmallToLarge() && edgeType1==4)) ){
448  if(ag==0){
449  cut2[3]=1;
450  if (pid==11) {cut2[6]=1;}
451  else if(pid==22) {cut2[7]=1;}
452  else if(pid==33) {cut2[8]=1;}
453  else {cut2[9]=1;}
454  }
455  if(ag==1) cut2[5]=1;
456  if(n1[2]==2) cut2[10]=1;
457  }
458  }
459  if( (edgeType0==4 && d0>=-0.51) || (edgeType1==4 && d1>=-0.51) ){
460  cut2[4]=1;
461  }
462  }
463  for(int c=0; c<NCUT2; c++){
464  if(cut2[c]==0) continue;
465  n2[c]++;
466  mPene[c]->Fill(pair->energy());
467  mPpt[c]->Fill(pair->pT());
468  mPept[c]->Fill(pair->energy(),pair->pT());
469  mPeta[c]->Fill(pair->eta());
470  mPphi[c]->Fill(pair->phi());
471  mPpid[c]->Fill(float(pid/10),float(pid%10));
472  mPm1[c]->Fill(pair->mass());
473  mPm2[c]->Fill(pair->mass());
474  mPzgg[c]->Fill(pair->zgg());
475  mPdgg[c]->Fill(pair->dgg());
476  mPr30[c]->Fill(pair->coneEnergyFraction(2));
477  mPr100[c]->Fill(pair->coneEnergyFraction(0));
478  mPxy[c]->Fill(pair->x(),pair->y());
479  }//loop over cuts
480  }//loop over pairs
481  for(int c=0; c<NCUT2; c++) mPn[c]->Fill(float(n2[c]));
482 
483  return kStOk;
484 }
485 
486 void StFmsOfflineQaMaker::print(){
487  mFmsColl->print(4);
488  StSPtrVecFpsSlat& slats = mFmsColl->fpsSlats();
489  for(int l=1; l<=kFpsNLayer; l++){
490  printf("NMIP FPS layer%1d ",l);
491  for(int q=1; q<=kFpsNQuad; q++){
492  printf(" Q%1d ",q);
493  for(int s=1; s<=kFpsNSlat; s++){
494  int slatid = mFmsDbMaker->fpsSlatId(q,l,s);
495  if(slatid<0) {;
496  printf("x");
497  }else{
498  float mip=slats[slatid]->mip();
499  int n=int(mip+0.5);
500  if(n>9) n=9;
501  printf("%1d",n);
502  }
503  }
504  }
505  printf("\n");
506  printf("NFmsPoint FPS layer%1d ",l);
507  for(int q=1; q<=kFpsNQuad; q++){
508  printf(" Q%1d ",q);
509  for(int s=1; s<=kFpsNSlat; s++){
510  int slatid = mFmsDbMaker->fpsSlatId(q,l,s);
511  if(slatid<0) {
512  printf("x");
513  }else{
514  int n=slats[slatid]->nPoint(0);
515  if(n>9) n=9;
516  printf("%1d",n);
517  }
518  }
519  }
520  printf("\n");
521  }
522 
523  /*
524  //testing StFmsDbMaker::distanceFromEdge()
525  const int nn=100;
526  const float minx=-10.0, maxx=110.0;
527  for(int i=0; i<nn; i++){
528  float x=minx + i*(maxx-minx)/float(nn);
529  float y=20.0;
530  int edge=0;
531  float d=mFmsDbMaker->distanceFromEdge(8,x,y,edge);
532  printf("EDGE i=%4d det=%2d x=%8.2f y=%8.2f d=%8.2f edge=%1d\n",i,8,x,y,d,edge);
533  }
534  for(int i=0; i<nn; i++){
535  float x=minx + i*(maxx-minx)/float(nn);
536  float y=80.0;
537  int edge=0;
538  float d=mFmsDbMaker->distanceFromEdge(8,x,y,edge);
539  printf("EDGE i=%4d det=%2d x=%8.2f y=%8.2f d=%8.2f edge=%1d\n",i,8,x,y,d,edge);
540  }
541  for(int i=0; i<nn; i++){
542  float x=minx + i*(maxx-minx)/float(nn);
543  float y=10.0;
544  int edge=0;
545  float d=mFmsDbMaker->distanceFromEdge(10,x,y,edge);
546  printf("EDGE i=%4d det=%2d x=%8.2f y=%8.2f d=%8.2f edge=%1d\n",i,10,x,y,d,edge);
547  }
548  for(int i=0; i<nn; i++){
549  float x=minx + i*(maxx-minx)/float(nn);
550  float y=40.0;
551  int edge=0;
552  float d=mFmsDbMaker->distanceFromEdge(10,x,y,edge);
553  printf("EDGE i=%4d det=%2d x=%8.2f y=%8.2f d=%8.2f edge=%1d\n",i,10,x,y,d,edge);
554  }
555  */
556 }
StFmsDbMaker::getXWidth
Float_t getXWidth(Int_t detectorId)
get the offset of the detector
Definition: StFmsDbMaker.cxx:1370
StFmsDbMaker::largeSmall
Int_t largeSmall(Int_t detectorId)
north or south side
Definition: StFmsDbMaker.cxx:1297
StFmsDbMaker::getRowNumber
Int_t getRowNumber(Int_t detectorId, Int_t ch)
maximum number of channels
Definition: StFmsDbMaker.cxx:1349
StFmsDbMaker::getColumnNumber
Int_t getColumnNumber(Int_t detectorId, Int_t ch)
get the row number for the channel
Definition: StFmsDbMaker.cxx:1354
StFmsDbMaker::distanceFromEdge
Float_t distanceFromEdge(Int_t det, Float_t x, Float_t y, int &edge)
get 4 vector assuing m=0 and taking beamline from DB
Definition: StFmsDbMaker.cxx:2142
StFmsDbMaker::northSouth
Int_t northSouth(Int_t detectorId)
east or west to the STAR IP
Definition: StFmsDbMaker.cxx:1289
TDataSet::ls
virtual void ls(Option_t *option="") const
Definition: TDataSet.cxx:495
StFmsDbMaker::getYWidth
Float_t getYWidth(Int_t detectorId)
get the X width of the cell
Definition: StFmsDbMaker.cxx:1375
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
kStOK
Definition: Stypes.h:40
kStErr
Definition: Stypes.h:44
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5