SiCalo/INTTCalo/EventDisplay.cpp

0001 #include <iostream>
0002 #include "INTT_Calo_trkReco.h"
0003 #include "EventDisplay.h"
0004
0005 const float PI = 3.141592;
0006
0007 void EventSelection( int& ievent) {
0008   // Menu
0009   int imenu;
0010   cout <<" <0: quit"<<endl;
0011   cout <<"  0: display current data"<<endl;
0012   cout <<"  1: display next data"<<endl;
0013   cout <<"  2: display specified event"<<endl;
0014   cin >> imenu;
0015   if(imenu<0) {
0016     cout <<" quit "<<endl;
0017     ievent = -1;
0018     return;
0019   }
0020   if(imenu==0) {
0021     cout <<"Event "<<ievent<<endl;
0022   }
0023   if(imenu==1) {
0024     ievent ++;
0025     cout <<"Event "<<ievent<<endl;
0026   }
0027   if(imenu==2) {
0028     cout <<"Enter event number to display (negative to quit) "<<endl;
0029     cin>>ievent;
0030     if(ievent<0 || ievent >40000) {
0031       ievent = -1;
0032       return;
0033     }
0034   }
0035   if(imenu>2) {
0036     cout <<" out of range. Keep current event"<<endl;
0037   }
0038
0039   cout<<" enter phiSlice"<<endl;
0040   cin>>phiSlice;
0041   cout << "phiSlice="<<phiSlice<<endl;
0042 }
0043
0044
0045 void EventDisplay(int iframe){
0046
0047   int nEmcClust=energy->size();
0048   int ntrack=pt0->size();
0049   int nSiClust=si_x->size();
0050   int nSiClustA=a_si_x->size();
0051
0052   if(nEmcClust>nMax) nEmcClust=nMax;
0053
0054   // Draw hits (EMC and Silicon) and seedtrk in c1
0055
0056   if(c1==nullptr) {
0057     c1 = new TCanvas("c1","XY View",500.,500.);
0058     frame1 = c1->DrawFrame(-240.,-240.,240.,240.);
0059   }
0060   if(c2==nullptr) {
0061     c2 = new TCanvas("c2","RZ View",500.,250.);
0062     frame2 = c2->DrawFrame(-240.,0.,240.,240.);
0063   }
0064   if(c3==nullptr) {
0065     c3 = new TCanvas("c3","XY View",500.,500.);
0066     c3->SetGridx();
0067     c3->SetGridy();
0068     frame3 = c3->DrawFrame(-12,-12,12,12);
0069   }
0070   if(c4==nullptr) {
0071     c4 = new TCanvas("c4","RZ View",750.,250.);
0072     frame4 = c4->DrawFrame(-18,0,18,12);
0073   }
0074
0075   if(nhits>0) {
0076     for(int ihit=0;ihit<nhits;ihit++) {
0077       hits[ihit]->Delete();
0078     }
0079     nhits=0;
0080   }
0081
0082   if(nhitsRZ>0) {
0083     for(int ihitRZ=0;ihitRZ<nhitsRZ;ihitRZ++) {
0084       hitsRZ[ihitRZ]->Delete();
0085     }
0086     nhitsRZ=0;
0087   }
0088
0089   if(ntrks>0) {
0090     for(int itrk=0;itrk<ntrks;itrk++) {
0091       seedtrk[itrk]->Delete();
0092     }
0093     ntrks=0;
0094   }
0095   if(ntrksRZ>0){
0096     for(int itrkRZ=0;itrkRZ<ntrksRZ;itrkRZ++) {
0097       seedtrkRZ[itrkRZ]->Delete();
0098     }
0099     ntrksRZ=0;
0100   }
0101
0102   float s0;    //size of EMC and top hit
0103   float s1;    //size of Silicon hit
0104   float s2;    //size of SiliconAll hit
0105   const float sLA=0.8;  //half length of strip. type-A
0106   const float sLB=1.0;  //half lenght of strip. type-B
0107
0108   s0=20.0;   // unit size of calo cluster
0109   if(iframe==0) {
0110     s1=0.5;
0111     s2=0.2;
0112   } else {
0113     s1=0.2;
0114     s2=0.1;
0115   }
0116
0117   if(abs(zgvtx)>50) {
0118     zgvtx = 0;
0119   }
0120
0121   // Create ZVTX marker
0122   if(ZvtxMarker !=nullptr) ZvtxMarker->Delete();
0123   ZvtxMarker= new TArrow(zgvtx,2.0,zgvtx,0.0,0.01);
0124   ZvtxMarker->SetLineWidth(2);
0125   ZvtxMarker->SetLineColor(kRed);
0126
0127   //  Create EM hits
0128   if(nEmcClust>0) {
0129     for(int j=0;j<nEmcClust;j++){
0130       x[j]= (emc_r->at(j))*cos(emc_phi->at(j));
0131       y[j]= (emc_r->at(j))*sin(emc_phi->at(j));
0132
0133       float sE = sqrt(s0*energy->at(j));
0134
0135       hits[nhits]=new TBox(x[j]-sE,y[j]-sE,x[j]+sE,y[j]+sE);
0136       hits[nhits]->SetFillColor(kBlue);
0137       nhits++;
0138
0139       if(abs(emc_phi->at(j)-phiSlice)<DphiSlice) {
0140     float Remc = emc_r->at(j);
0141     float Zemc = emc_z->at(j);
0142     hitsRZ[nhitsRZ]=new TBox(Zemc-sE,Remc-sE,Zemc+sE,Remc+sE);
0143     hitsRZ[nhitsRZ]->SetFillColor(kBlue);
0144     nhitsRZ++;
0145       }//for
0146     }//if(nEmcClust)
0147   }
0148
0149   // Create topoCluster hits
0150   int nTop=top_e->size();
0151   if(nTop>0) {
0152     for(int itop=0;itop< nTop ;itop++){
0153       x[itop] = top_x->at(itop);
0154       y[itop] = top_y->at(itop);
0155       float sE = sqrt(s0*top_e->at(itop));
0156
0157       hits[nhits]=new TBox(x[itop]-sE,y[itop]-sE,x[itop]+sE,y[itop]+sE);
0158       hits[nhits]->SetFillColor(kRed);
0159       nhits++;
0160
0161       if(abs(top_phi->at(itop)-phiSlice)<DphiSlice) {
0162     float Rtop = top_r->at(itop);
0163     float Ztop = top_z->at(itop);
0164     hitsRZ[nhitsRZ]=new TBox(Ztop-sE,Rtop-sE,Ztop+sE,Rtop+sE);
0165     hitsRZ[nhitsRZ]->SetFillColor(kRed);
0166     nhitsRZ++;
0167       }//if(phiSlice)
0168     }//for itop;
0169   }//if(nTop)
0170
0171   // Creates Silicon hits
0172   // Note that the track reconstruciton is done after shifting
0173   // the beam center (xBC, yBC) to (0.0). To be consitent with
0174   // tracks, silicon hit positions are shifted by the same amount
0175   // in the event display.
0176   if(nSiClust>0) {
0177     for(int jj=0;jj<nSiClust;jj++) {
0178       x[jj]=si_x->at(jj) - xBC;
0179       y[jj]=si_y->at(jj) - yBC;
0180       hits[nhits]=new TBox(x[jj]-s1,y[jj]-s1,x[jj]+s1,y[jj]+s1);
0181       if(layer->at(jj) < 3) hits[nhits]->SetFillColor(kGreen);
0182       else          hits[nhits]->SetFillColor(kRed);
0183       nhits++;
0184
0185       float phiHit=atan2(y[jj],x[jj]);
0186       if(abs(phiHit-phiSlice)<DphiSlice) {
0187     float Rhit =  sqrt(x[jj]*x[jj]+y[jj]*y[jj]);
0188     float Zhit =  si_z->at(jj);
0189
0190     if(layer->at(jj) < 3) { //MVTX hit
0191       hitsRZ[nhitsRZ]=new TBox(Zhit-s1,Rhit-s1,Zhit+s1,Rhit+s1);
0192       hitsRZ[nhitsRZ]->SetFillColor(kGreen);
0193     } else { //INTT hit
0194       if(abs(Zhit)<12.8) { // type-A strip
0195         hitsRZ[nhitsRZ]=new TBox(Zhit-sLA,Rhit-s1,Zhit+sLA,Rhit+s1);
0196       } else {  // type-B strip
0197         hitsRZ[nhitsRZ]=new TBox(Zhit-sLB,Rhit-s1,Zhit+sLB,Rhit+s1);
0198       }
0199       hitsRZ[nhitsRZ]->SetFillColor(kRed);
0200     }
0201     nhitsRZ++;
0202       }//for(jj)
0203     }//if(nSiClust)
0204   }
0205
0206   // Creates SiliconAll hits
0207   if(nSiClustA>0) {
0208     for(int jj=0;jj<nSiClustA;jj++) {
0209       x[jj]=a_si_x->at(jj) - xBC;
0210       y[jj]=a_si_y->at(jj) - yBC;
0211       hits[nhits]=new TBox(x[jj]-s2,y[jj]-s2,x[jj]+s2,y[jj]+s2);
0212       if(a_layer->at(jj) < 3) hits[nhits]->SetFillColor(kBlack);
0213       else          hits[nhits]->SetFillColor(kBlack);
0214       nhits++;
0215
0216       float phiHit=atan2(y[jj],x[jj]);
0217       if(abs(phiHit-phiSlice)<DphiSlice) {
0218     float Rhit =  sqrt(x[jj]*x[jj]+y[jj]*y[jj]);
0219     float Zhit =  a_si_z->at(jj);
0220
0221     if(a_layer->at(jj) < 3) { // MVTX hit
0222       hitsRZ[nhitsRZ]=new TBox(Zhit-s2,Rhit-s2,Zhit+s2,Rhit+s2);
0223       hitsRZ[nhitsRZ]->SetFillColor(kBlack);
0224     } else {// INTT hit
0225       if(abs(Zhit)<12.8) { // type-A strip
0226         hitsRZ[nhitsRZ]=new TBox(Zhit-sLA,Rhit-s2,Zhit+sLA,Rhit+s2);
0227       } else { // type-B strip
0228         hitsRZ[nhitsRZ]=new TBox(Zhit-sLB,Rhit-s2,Zhit+sLB,Rhit+s2);
0229       }
0230       hitsRZ[nhitsRZ]->SetFillColor(kBlack);
0231     }
0232     nhitsRZ++;
0233       }
0234     }//if(nSiClustA)
0235   }
0236
0237   // Creates SiSeed tracklets (arrow)
0238   float Rtrk;
0239   const float Asize = 0.02;
0240   if(iframe==1) Rtrk=2.0; else Rtrk=80;
0241   if(iframe==1) Rmax=12.0; else Rmax=80;
0242   for(int k=0;k<ntrack;k++) {
0243     if(pt0->at(k)>pt0min) {
0244       float dx=Rtrk*cos(phi0->at(k));
0245       float dy=Rtrk*sin(phi0->at(k));
0246       seedtrk[ntrks]=new TArrow(vx0->at(k),vy0->at(k),vx0->at(k)+dx,vy0->at(k)+dy,Asize);
0247       ntrks++;
0248
0249       if(abs(phi0->at(k)-phiSlice)<DphiSlice) {
0250     float r0=sqrt((vx0->at(k))*(vx0->at(k))+(vy0->at(k))*(vy0->at(k)));
0251     float eta = eta0->at(k);
0252     float theta = 2*atan(exp(-eta));
0253     float Zend = (Rmax-r0)/tan(theta)+vz0->at(k);
0254
0255     seedtrkRZ[ntrksRZ]=new TArrow(vz0->at(k),r0,Zend,Rmax,Asize);
0256     ntrksRZ++;
0257       }
0258     }
0259   }
0260
0261   // clear track orbits of the previous event;
0262   if(vEmcINTTorbit.size()>0) {
0263     for(int io=0;io<(int)vEmcINTTorbit.size();io++) {
0264       vEmcINTTorbit.at(io)->Delete();
0265     }//for(io)
0266   }//if
0267   vEmcINTTorbit.clear();
0268
0269   if(vEmcINTT_RZ.size()>0) {
0270     for(int io=0;io<(int)vEmcINTT_RZ.size();io++) {
0271       vEmcINTT_RZ.at(io)->Delete();
0272     }
0273   }
0274   vEmcINTT_RZ.clear();
0275   // create Circle orbit of CaloInttMvtx tracks
0276   int ntrk = vCaloInttMvtx.size();
0277
0278   if(ntrk>0) {
0279     for(int itrk=0;itrk<ntrk;itrk++) {
0280       CaloInttMvtx CIMtrack = vCaloInttMvtx.at(itrk);
0281       // create emc_intt tracks
0282       float r_e  = CIMtrack.r_emc_intt;
0283       if(r_e>0) { //this track is a EMC-INTT
0284     float xc_e = CIMtrack.xc_emc_intt-xBC;
0285     float yc_e = CIMtrack.yc_emc_intt-yBC;
0286     float xemc = CIMtrack.xemc;
0287     float yemc = CIMtrack.yemc;
0288         float zemc = CIMtrack.zemc;
0289     float phi0e = CIMtrack.phi0e;
0290     float remc = sqrt(xemc*xemc+yemc*yemc);
0291     float zvtx_e = CIMtrack.zvtx_emc_intt;
0292     float zsign_e = xc_e*yemc-yc_e*xemc;
0293     float phi0_e = atan2(-yc_e,-xc_e)*180./PI;
0294     TArc *orbit_e;
0295     if(zsign_e<0) {
0296       orbit_e = new TArc(xc_e,yc_e,r_e,phi0_e,phi0_e+10000/r_e);
0297     } else {
0298       orbit_e = new TArc(xc_e,yc_e,r_e,phi0_e,phi0_e-10000/r_e);
0299     }
0300     orbit_e->SetFillStyle(0);
0301     orbit_e->SetLineColor(kGreen);
0302     vEmcINTTorbit.push_back(orbit_e);
0303
0304     //create RZ view of the track
0305     if(abs(phi0e - phiSlice)<DphiSlice) {
0306       TArrow *RZ_trk = new TArrow(zvtx_e,0,zemc,remc,0.03);
0307       RZ_trk->SetLineColor(kGreen);
0308       vEmcINTT_RZ.push_back(RZ_trk);
0309     }
0310       }//if(r_e>0)
0311
0312       // create intt-mvtx tracks
0313       float r_m  = CIMtrack.r_intt_mvtx;
0314       if(r_m>0) {  //this track is a INTT-Mvtx track
0315     float xc_m = CIMtrack.xc_intt_mvtx;
0316     float yc_m = CIMtrack.yc_intt_mvtx;
0317     float x_m = CIMtrack.x0m;
0318     float y_m = CIMtrack.y0m;
0319     float z_m = CIMtrack.z0m;
0320     float phi0m = CIMtrack.phi0m;
0321     float r0m = sqrt(x_m*x_m+y_m*y_m);
0322     float zvtx_m = CIMtrack.zvtx_intt_mvtx;
0323     float dzdr = (z_m-zvtx_m)/r0m;
0324     float zsign_m = xc_m*y_m-yc_m*x_m;
0325     float phi0_m = atan2(-yc_m,-xc_m)*180./PI;
0326     const float R_emc = 100.0;
0327     TArc *orbit_m;
0328     if(zsign_m<0) {
0329       orbit_m = new TArc(xc_m,yc_m,r_m,phi0_m,phi0_m+10000/r_m);
0330     } else {
0331       orbit_m = new TArc(xc_m,yc_m,r_m,phi0_m,phi0_m-10000/r_m);
0332     }
0333     orbit_m->SetFillStyle(0);
0334     orbit_m->SetLineColor(kBlue);
0335     vEmcINTTorbit.push_back(orbit_m);
0336
0337     //create RZ view of the track
0338     if(abs(phi0m - phiSlice)<DphiSlice) {
0339       TArrow *RZ_trk = new TArrow(zvtx_m,0,zvtx_m+dzdr*R_emc,R_emc,0.03);
0340       RZ_trk->SetLineColor(kBlue);
0341       vEmcINTT_RZ.push_back(RZ_trk);
0342     }
0343       }//if(r_m)
0344     }//for(itrk)
0345   }
0346   // Draw X-Y view
0347   if(iframe==0) {
0348     c1->cd();
0349   } else {
0350     c3->cd();
0351   }
0352
0353   for(int kk=0;kk<nhits;kk++){
0354     hits[kk]->Draw();
0355   }
0356   for(int itk=0;itk<ntrks;itk++) {
0357     seedtrk[itk]->Draw();
0358   }
0359   // Draw phi slice region
0360   if(dphiArc != nullptr) dphiArc->Delete();
0361   float phimin=(phiSlice-DphiSlice)*180./PI;
0362   float phimax=(phiSlice+DphiSlice)*180./PI;
0363   dphiArc = new TArc(0.,0.,Rmax,phimin,phimax);
0364   dphiArc->SetLineColorAlpha(0,0.1);
0365   dphiArc->SetFillColorAlpha(kBlue,0.1);
0366   dphiArc->Draw();
0367
0368   // Draw Circle for reference size of INTT and MVTX
0369   if(Circle2==nullptr) {
0370     Circle2=new TArc(0.,0.,2.3,0,360);
0371     Circle2->SetFillStyle(0);
0372   }
0373   if(Circle4==nullptr) {
0374     Circle4=new TArc(0.,0.,4.0,0,360);
0375     Circle4->SetFillStyle(0);
0376   }
0377   if(Circle8==nullptr) {
0378     Circle8=new TArc(0.,0.,7.6,0.,360);
0379     Circle8->SetFillStyle(0);
0380   }
0381   if(Circle10==nullptr) {
0382     Circle10=new TArc(0.,0.,10.2,0.,360);
0383     Circle10->SetFillStyle(0);
0384   }
0385   Circle2->Draw();
0386   Circle4->Draw();
0387   Circle8->Draw();
0388   Circle10->Draw();
0389
0390   // draw EmcINTT tracks
0391   int norbit = vEmcINTTorbit.size();
0392   if(norbit>0) {
0393     for(int i=0;i<norbit;i++) {
0394       vEmcINTTorbit.at(i)->Draw("only");
0395     }//for
0396   }//if
0397   if(iframe==0) {
0398     //    c1->Draw();
0399     c1->Modified();
0400     c1->Update();
0401   } else {
0402     //    c3->Draw();
0403     c3->Modified();
0404     c3->Update();
0405   }
0406   //Draw RZ slice
0407   if(iframe==0) {
0408     c2->cd();
0409   } else {
0410     c4->cd();
0411   }
0412
0413   ZvtxMarker->Draw();
0414   for(int kk=0;kk<nhitsRZ;kk++){
0415     hitsRZ[kk]->Draw();
0416   }
0417   for(int itk=0;itk<ntrksRZ;itk++) {
0418     seedtrkRZ[itk]->Draw();
0419   }
0420   int norbitRZ = vEmcINTT_RZ.size();
0421   if(norbitRZ>0) {
0422     for(int i=0;i<norbitRZ;i++) {
0423       vEmcINTT_RZ.at(i)->Draw();
0424     }//for(i)
0425   }//if(ntrk)
0426   if(iframe==0) {
0427     c2->Modified();
0428     c2->Update();
0429   } else {
0430     c4->Modified();
0431     c4->Update();
0432   }
0433 }