SiCalo/INTTCalo/display_data.cpp

0001
0002 #include <vector>
0003 #include <iostream>
0004 #include <TFile.h>
0005 #include <TTree.h>
0006 #include <TCanvas.h>
0007 #include <TH1F.h>
0008 #include <TGraph.h>
0009 #include <TArc.h>
0010 #include <TArrow.h>
0011 #include <TBox.h>
0012
0013 //
0014 // Variables for Event Display
0015 //
0016 const float PI = 3.141592;
0017
0018 string idatafile;   //filename of the input data file
0019 int NmaxEvent;      // # of events in the data file
0020 int idata_Type;     // data type. 1=Hachiya 2=Jingyu
0021
0022
0023 //Global variables
0024 //
0025 // imode==1 : run AnalyzeEvent()
0026 // imode==2 : Event Display
0027 int imode = 2;
0028
0029 // display frame;
0030 static int iframe=-1; //0: All view 1: Si View
0031 static TH1F *frame1 = nullptr;
0032 static TH1F *frame2 = nullptr;
0033 static TCanvas *c1=nullptr;
0034 static TCanvas *c2=nullptr;
0035
0036 static TArc *Circle2 = nullptr;
0037 static TArc *Circle4 = nullptr;
0038 static TArc *Circle8 = nullptr;
0039 static TArc *Circle10 = nullptr;
0040
0041 //pointer to the TTree for the data.0
0042 TTree *evtTree;
0043 TTree *caloTree;  //EMCal clusters
0044 TTree *topoTree;  //topoClusters (EMC+iHCal+oHCal)
0045 TTree *SiClusTree;
0046 TTree *SiClusAllTree;
0047 TTree *trackTree;
0048
0049 // data in evtTree
0050 Int_t evt;
0051 float xvtx;  //SiSeed X
0052 float yvtx;  //SiSeed Y
0053 float zvtx;  //SiSeed Z
0054 float gxvtx;  //Global vertex X
0055 float gyvtx;  //Global vertex Y
0056 float gzvtx;  //Global vertex Z
0057
0058 // data in caloTree
0059 vector<float> *energy=nullptr;
0060 vector<float> *emc_phi=nullptr;
0061 vector<float> *emc_r=nullptr;
0062 vector<float> *emc_z  =nullptr;
0063
0064 // data in topoTree
0065 vector<float> *top_e  =nullptr;
0066 vector<float> *top_x  =nullptr;
0067 vector<float> *top_y  =nullptr;
0068 vector<float> *top_z  =nullptr;
0069 vector<float> *top_r  =nullptr;
0070 vector<float> *top_phi=nullptr;
0071 vector<float> *emc_e  =nullptr;
0072 vector<float> *ihc_e  =nullptr;
0073 vector<float> *ohc_e  =nullptr;
0074
0075 // data in SiClustTree
0076 vector<int> *layer=nullptr;
0077 vector<int> *trkid=nullptr;
0078 vector<float> *si_x=nullptr;
0079 vector<float> *si_y=nullptr;
0080 vector<float> *si_z=nullptr;
0081
0082 // data in SiClustAllTree
0083 vector<int> *a_layer=nullptr;
0084 vector<int> *a_trkid=nullptr;
0085 vector<float> *a_si_x=nullptr;
0086 vector<float> *a_si_y=nullptr;
0087 vector<float> *a_si_z=nullptr;
0088
0089 // data in trackTree
0090 vector<float> *eta0=nullptr;
0091 vector<float> *phi0=nullptr;
0092 vector<float> *pt0=nullptr;
0093 vector<float> *vx0=nullptr;
0094 vector<float> *vy0=nullptr;
0095 vector<float> *vz0=nullptr;
0096
0097 // buffer for hits to display
0098 const int nMax=4000;
0099 float x[nMax], y[nMax];
0100 TBox *hits[nMax];
0101 TBox *hitsRZ[nMax];
0102 static int nhits = 0;
0103 static int nhitsRZ = 0;
0104
0105 // ZVTX maker
0106 TArrow *ZvtxMarker = nullptr;
0107
0108 // buffer for tracks to display
0109 const int nMaxTrk=300;
0110 TArrow *seedtrk[nMaxTrk];
0111 TArrow *seedtrkRZ[nMaxTrk];
0112 static int ntrks = 0;
0113 static int ntrksRZ = 0;
0114
0115 // buffer for EMC-INTT tracklet orbit
0116 vector<TArc*> vEmcINTTorbitL;   //Local Orbit
0117 vector<TArc*> vEmcINTTorbitG;   //Global Orbit
0118
0119 // constants
0120 const float pt0min=0.5; // GeV/c. minimum pT0 for tracks
0121 const float EMScale=0.14;  //EMC hit scaling factor
0122
0123 // RZ slice angle
0124 static float phiSlice = 3.141592/2;
0125 const float DphiSlice = 0.3;
0126 static float Rmax = 12.0;
0127 static TArc *dphiArc = nullptr;
0128
0129 //
0130 // variables for EMC-INTT tracking
0131 //
0132 struct EmcINTTtriplet {
0133   int iEmc;       //index in vEMr, vEMphi, vEMz
0134   int iINTT0;     //index in vINTT0r, vINTT0phi, vINTT0z
0135   int iINTT1;     //index in vINTT1r, vINTT1phi, vINTT1z
0136 };
0137
0138 // Emc-INTT1-INTT0 triplet
0139 vector<EmcINTTtriplet> vEmcINTT;  //triplet
0140 vector<int> vEmcINTT_Mvtx1;  //index of Mvtx1 hit assciated with triplet
0141 vector<int> vEmcINTT_Mvtx2;  //index of Mvtx2 hit assciated with triplet
0142 vector<int> vEmcINTT_Mvtx3;  //index of Mvtx3 hit assciated with triplet
0143 vector<float> vEmcINTT_R;    //local radius of EmcINTT tracklet orbit (circle)
0144 vector<float> vEmcINTT_xO;   // x of the center of the tracklet orbit
0145 vector<float> vEmcINTT_yO;   // y of the center of the tracklet orbit
0146 vector<float> vEmcINTT_R_EI; //local radius of EmcINTT tracklet orbit (circle)
0147 vector<float> vEmcINTT_xO_EI;// x of the center of the tracklet orbit
0148 vector<float> vEmcINTT_yO_EI;// y of the center of the tracklet orbit
0149
0150 vector<float> vEMr;
0151 vector<float> vEMphi;
0152 vector<float> vEMz;
0153 vector<float> vEMe;
0154
0155 vector<float> vINTT0r;
0156 vector<float> vINTT0phi;
0157 vector<float> vINTT0z;
0158
0159 vector<float> vINTT1r;
0160 vector<float> vINTT1phi;
0161 vector<float> vINTT1z;
0162
0163 vector<float> vMVTXr;
0164 vector<float> vMVTXphi;
0165 vector<float> vMVTXx;
0166 vector<float> vMVTXy;
0167 vector<float> vMVTXz;
0168
0169
0170 void LoadData(string filename);
0171 void GetEvent(int ievent);
0172 void EventSelection(int&);
0173 void EventDisplay(void);
0174
0175 void AnalyzeInit(void);
0176 void AnalyzeEvent(void);
0177 void AnalyzeEnd(void);
0178 void LoadTrackerData(void);
0179 void FindEmcINTTtriplet(void);
0180 void INTT_zvtx_Finder(float&);
0181 void FindEmcINTT_Mvtx(void);
0182
0183 void display_data(void){
0184   int ievent = 0;
0185   const int MaxLoop = 100;
0186
0187
0188   // open data file and load the data
0189
0190   //input data file
0191   //  idatafile = "run3pp/79516_old/ana_79516_00000.root";
0192   idatafile = "run2pp/53876/ana_53876_00001.root";
0193   // idatafile ="ana_53879_1kevt_4evtdisplay.root";
0194
0195   // type of the data
0196   idata_Type = 1;  // Hachiya nDST
0197   //  idata_Type = 2;  // Jingy's new nDST
0198
0199   LoadData(idatafile.c_str());  //field ON and new
0200   NmaxEvent = evtTree->GetEntries();
0201
0202   cout <<"Data file:"<<idatafile.c_str()<<endl;
0203   cout <<"# of events in the file:"<<NmaxEvent<<endl;
0204
0205   cout << "Select Mode: " <<endl;
0206   cout << " 1: run AnalyzeEvent" << endl;
0207   cout << " 2: event display" <<endl;
0208   cin >> imode;
0209   if(imode == 1) {
0210     int istart = 0;
0211     int nevent=0;
0212
0213     cout << "enter istart (first event)"<<endl;
0214     cin >> istart;
0215     if(istart<0) istart=0;
0216
0217     cout << "how many event:"<<endl;
0218     cin >> nevent;
0219     if(nevent<0) nevent=0;
0220     if(istart+nevent > NmaxEvent) {
0221       cout <<" Exceeds the # of events in the file" <<endl;
0222       nevent = NmaxEvent - istart;
0223     }
0224     cout << "analyze "<< nevent << " events" <<endl;
0225
0226     AnalyzeInit();
0227     for(int i=istart;i<istart+nevent;i++) {
0228       //      cout << "event = "<<i << endl;
0229       GetEvent(i);
0230       AnalyzeEvent();
0231       if(vEmcINTT.size()>0) {
0232     cout <<"ievent="<<i<<" # of triplet="<<vEmcINTT.size()<<endl<<endl;
0233       }
0234       if(i%1000 == 0) cout << i << endl;
0235     }//for
0236     AnalyzeEnd();
0237   } else if(imode == 2 ) { // Event Display
0238     gStyle->SetCanvasPreferGL();
0239
0240     AnalyzeInit();
0241     for(int i=0;i<MaxLoop;i++) {
0242       EventSelection( ievent);
0243       if(ievent<0) break;
0244       GetEvent(ievent);
0245       AnalyzeEvent();  // do additional analysis
0246       EventDisplay();
0247     }//for
0248   }//if
0249 }
0250
0251 void EventSelection( int& ievent) {
0252   // Menu
0253   int imenu;
0254   cout <<" <0: quit"<<endl;
0255   cout <<"  0: display current data"<<endl;
0256   cout <<"  1: display next data"<<endl;
0257   cout <<"  2: display specified event"<<endl;
0258   cin >> imenu;
0259   if(imenu<0) {
0260     cout <<" quit "<<endl;
0261     ievent = -1;
0262     return;
0263   }
0264   if(imenu==0) {
0265     cout <<"Event "<<ievent<<endl;
0266   }
0267   if(imenu==1) {
0268     ievent ++;
0269     cout <<"Event "<<ievent<<endl;
0270   }
0271   if(imenu==2) {
0272     cout <<"Enter event number to display (negative to quit) "<<endl;
0273     cin>>ievent;
0274     if(ievent<0 || ievent >40000) {
0275       ievent = -1;
0276       return;
0277     }
0278   }
0279   if(imenu>2) {
0280     cout <<" out of range. Keep current event"<<endl;
0281   }
0282
0283   cout<<" enter phiSlice"<<endl;
0284   cin>>phiSlice;
0285   cout << "phiSlice="<<phiSlice<<endl;
0286 }
0287
0288 bool Is_Hot_EMC(int iemc) {
0289   float phiEMC = emc_phi->at(iemc);
0290   float zEMC = emc_z->at(iemc);
0291   if( (2.4 < phiEMC && phiEMC < 2.7) && (zEMC < 0)) return true;
0292   else return false;
0293 }
0294
0295 void LoadData(string filename) {
0296   TFile *file = TFile::Open(filename.c_str(),"READ");
0297   if(!file || file->IsZombie()) return;
0298
0299   // connect to the TTrees in the nDST file
0300   evtTree       = dynamic_cast<TTree*>(file->Get("evtTree"));
0301   caloTree      = dynamic_cast<TTree*>(file->Get("caloTree"));
0302   SiClusTree    = dynamic_cast<TTree*>(file->Get("SiClusTree"));
0303   SiClusAllTree = dynamic_cast<TTree*>(file->Get("SiClusAllTree"));
0304   trackTree     = dynamic_cast<TTree*>(file->Get("trackTree"));
0305   topoTree      = dynamic_cast<TTree*>(file->Get("TopoClusTree"));
0306
0307   //data in evtTree
0308   evtTree->SetBranchAddress("evt",&evt);
0309   evtTree->SetBranchAddress("xgvtx",&xvtx);
0310   evtTree->SetBranchAddress("ygvtx",&yvtx);
0311   if(idata_Type == 1) {
0312     evtTree->SetBranchAddress("zgvtx",&zvtx);
0313   } else if(idata_Type == 2) {
0314     evtTree->SetBranchAddress("zvtx",&zvtx);
0315   }
0316
0317   //data in caloTree
0318   if(idata_Type == 1) {// data by Hachiya
0319     caloTree->SetBranchAddress("energy",&energy);
0320     caloTree->SetBranchAddress("phi",&emc_phi);
0321     caloTree->SetBranchAddress("r",&emc_r);
0322     caloTree->SetBranchAddress("z",&emc_z);
0323   } else if(idata_Type ==2) {// data by Jingyu. This includes topoTree
0324     caloTree->SetBranchAddress("Calo_Clus_energy",&energy);
0325     caloTree->SetBranchAddress("Calo_Clus_phi",&emc_phi);
0326     caloTree->SetBranchAddress("Calo_Clus_r",&emc_r);
0327     caloTree->SetBranchAddress("Calo_Clus_z",&emc_z);
0328   }
0329   // now Hachiya data includes topoTree
0330     topoTree->SetBranchAddress("clus_e",&top_e);
0331     topoTree->SetBranchAddress("clus_x",&top_x);
0332     topoTree->SetBranchAddress("clus_y",&top_y);
0333     topoTree->SetBranchAddress("clus_z",&top_z);
0334     topoTree->SetBranchAddress("clus_r",&top_r);
0335     topoTree->SetBranchAddress("clus_z",&top_z);
0336     topoTree->SetBranchAddress("clus_phi",&top_phi);
0337     topoTree->SetBranchAddress("clus_emcal_e",&emc_e);
0338     topoTree->SetBranchAddress("clus_ihcal_e",&ihc_e);
0339     topoTree->SetBranchAddress("clus_ohcal_e",&ohc_e);
0340
0341   //data in SiClusTree
0342   SiClusTree->SetBranchAddress("Siclus_layer",&layer);
0343   SiClusTree->SetBranchAddress("Siclus_trackid",&trkid);
0344   SiClusTree->SetBranchAddress("Siclus_x",&si_x);
0345   SiClusTree->SetBranchAddress("Siclus_y",&si_y);
0346   SiClusTree->SetBranchAddress("Siclus_z",&si_z);
0347
0348   //data in SiClusAllTree
0349   SiClusAllTree->SetBranchAddress("Siclus_layer",&a_layer);
0350   SiClusAllTree->SetBranchAddress("Siclus_trackid",&a_trkid);
0351   SiClusAllTree->SetBranchAddress("Siclus_x",&a_si_x);
0352   SiClusAllTree->SetBranchAddress("Siclus_y",&a_si_y);
0353   SiClusAllTree->SetBranchAddress("Siclus_z",&a_si_z);
0354
0355   //data in trackTree
0356   trackTree->SetBranchAddress("eta0",&eta0);
0357   trackTree->SetBranchAddress("phi0",&phi0);
0358   trackTree->SetBranchAddress("pt0",&pt0);
0359   trackTree->SetBranchAddress("x0",&vx0);
0360   trackTree->SetBranchAddress("y0",&vy0);
0361   trackTree->SetBranchAddress("z0",&vz0);
0362 }
0363
0364 void GetEvent(int i) {
0365   evtTree->GetEntry(i);
0366   caloTree->GetEntry(i);
0367   trackTree->GetEntry(i);
0368   SiClusTree->GetEntry(i);
0369   SiClusAllTree->GetEntry(i);
0370 }
0371
0372 void EventDisplay(void){
0373
0374   int nEmcClust=energy->size();
0375   int ntrack=pt0->size();
0376   int nSiClust=si_x->size();
0377   int nSiClustA=a_si_x->size();
0378
0379   if(nEmcClust>nMax) nEmcClust=nMax;
0380
0381   // Draw hits (EMC and Silicon) and seedtrk in c1
0382
0383   { // select frame.
0384     int new_frame;
0385     cout << "select frame (0 or 1)" << endl;
0386     cin >> new_frame;
0387     if(new_frame != 0) new_frame=1;
0388     if(new_frame != iframe) {
0389       cout << "frame is changed to "<<new_frame<<endl;
0390       iframe = new_frame;
0391       if(c1==nullptr){
0392     c1 = new TCanvas("c1","XY View",800.,800.);
0393       } else {
0394     c1->cd();
0395       }
0396       if(iframe==0) {
0397     frame1 = c1->DrawFrame(-160.,-160.,160.,160.);
0398       } else if(iframe==1) {
0399     frame1 = c1->DrawFrame(-16.,-16.,16.,16.);
0400       }
0401
0402       if(c2==nullptr){
0403     c2= new TCanvas("c2","RZ Slice",800,270);
0404       } else {
0405     c2->cd();
0406       }
0407       if(iframe==0) {
0408     frame2 = c2->DrawFrame(-240.,0.,240.,160.);
0409       } else if(iframe==1) {
0410     frame2 = c2->DrawFrame(-24.,0.,24.,16.);
0411       }
0412     }
0413   }
0414
0415   if(nhits>0) {
0416     for(int ihit=0;ihit<nhits;ihit++) {
0417       hits[ihit]->Delete();
0418     }
0419     nhits=0;
0420   }
0421
0422   if(nhitsRZ>0) {
0423     for(int ihitRZ=0;ihitRZ<nhitsRZ;ihitRZ++) {
0424       hitsRZ[ihitRZ]->Delete();
0425     }
0426     nhitsRZ=0;
0427   }
0428
0429   //  cout << "ntrks="<<ntrks<<endl;
0430   if(ntrks>0) {
0431     for(int itrk=0;itrk<ntrks;itrk++) {
0432       seedtrk[itrk]->Delete();
0433     }
0434     ntrks=0;
0435   }
0436   if(ntrksRZ>0){
0437     for(int itrkRZ=0;itrkRZ<ntrksRZ;itrkRZ++) {
0438       seedtrkRZ[itrkRZ]->Delete();
0439     }
0440     ntrksRZ=0;
0441   }
0442
0443   float s0;    //size of EMC hit
0444   float s1;    //size of Silicon hit
0445   float s2;    //size of SiliconAll hit
0446   const float sLA=0.8;  //half length of strip. type-A
0447   const float sLB=1.0;  //half lenght of strip. type-B
0448   if(iframe==0) {
0449     s0=1.0;
0450     s1=0.5;
0451     s2=0.2;
0452   } else {
0453     s0=0.3;
0454     s1=0.2;
0455     s2=0.1;
0456   }
0457
0458   /*
0459   cout <<evt<< ":   nEmcClust="<<nEmcClust;
0460   cout <<" nSiClust="<<nSiClust;
0461   cout <<" nSiClustA="<<nSiClustA;
0462   cout <<" ntrack="<<ntrack;
0463   cout <<" zvtx="<<zvtx;
0464   if(abs(zvtx)>50) {
0465     zvtx = 0;
0466     cout << "--> "<<zvtx;
0467   }
0468   cout <<endl;
0469   */
0470
0471   // Create ZVTX marker
0472   if(ZvtxMarker !=nullptr) ZvtxMarker->Delete();
0473   ZvtxMarker= new TArrow(zvtx,2.0,zvtx,0.0,0.01);
0474   ZvtxMarker->SetLineWidth(2);
0475   ZvtxMarker->SetLineColor(kRed);
0476
0477   //  Create EM hits
0478   for(int j=0;j<nEmcClust;j++){
0479     if(!Is_Hot_EMC(j)) {
0480       x[j]= (emc_r->at(j))*cos(emc_phi->at(j));
0481       y[j]= (emc_r->at(j))*sin(emc_phi->at(j));
0482
0483
0484       float sE = sqrt(s0*energy->at(j));
0485       // rescale EMC hits if iframe==1
0486       if(iframe == 1) {
0487     x[j]=EMScale*x[j];
0488     y[j]=EMScale*y[j];
0489       }
0490       hits[nhits]=new TBox(x[j]-sE,y[j]-sE,x[j]+sE,y[j]+sE);
0491       hits[nhits]->SetFillColor(kBlue);
0492       nhits++;
0493
0494       if(abs(emc_phi->at(j)-phiSlice)<DphiSlice) {
0495     float Remc = emc_r->at(j);
0496     float Zemc = emc_z->at(j);
0497     if(iframe == 1) {
0498       Remc = EMScale*Remc;
0499       Zemc = EMScale*Zemc;
0500     }
0501     hitsRZ[nhitsRZ]=new TBox(Zemc-sE,Remc-sE,Zemc+sE,Remc+sE);
0502     hitsRZ[nhitsRZ]->SetFillColor(kBlue);
0503     nhitsRZ++;
0504       }//if(!Is_Hot_EMC)
0505     }//for
0506   }
0507
0508   // Creates Silicon hits
0509   for(int jj=0;jj<nSiClust;jj++) {
0510     x[jj]=si_x->at(jj);
0511     y[jj]=si_y->at(jj);
0512     hits[nhits]=new TBox(x[jj]-s1,y[jj]-s1,x[jj]+s1,y[jj]+s1);
0513     if(layer->at(jj) < 3) hits[nhits]->SetFillColor(kGreen);
0514     else          hits[nhits]->SetFillColor(kRed);
0515     nhits++;
0516
0517     float phiHit=atan2(y[jj],x[jj]);
0518     if(abs(phiHit-phiSlice)<DphiSlice) {
0519       float Rhit =  sqrt(x[jj]*x[jj]+y[jj]*y[jj]);
0520       float Zhit =  si_z->at(jj);
0521
0522       if(layer->at(jj) < 3) { //MVTX hit
0523     hitsRZ[nhitsRZ]=new TBox(Zhit-s1,Rhit-s1,Zhit+s1,Rhit+s1);
0524     hitsRZ[nhitsRZ]->SetFillColor(kGreen);
0525       } else { //INTT hit
0526     if(abs(Zhit)<12.8) { // type-A strip
0527       hitsRZ[nhitsRZ]=new TBox(Zhit-sLA,Rhit-s1,Zhit+sLA,Rhit+s1);
0528     } else {  // type-B strip
0529       hitsRZ[nhitsRZ]=new TBox(Zhit-sLB,Rhit-s1,Zhit+sLB,Rhit+s1);
0530     }
0531     hitsRZ[nhitsRZ]->SetFillColor(kRed);
0532       }
0533       nhitsRZ++;
0534     }
0535   }
0536
0537   // Creates SiliconAll hits
0538   for(int jj=0;jj<nSiClustA;jj++) {
0539     x[jj]=a_si_x->at(jj);
0540     y[jj]=a_si_y->at(jj);
0541     hits[nhits]=new TBox(x[jj]-s2,y[jj]-s2,x[jj]+s2,y[jj]+s2);
0542     if(a_layer->at(jj) < 3) hits[nhits]->SetFillColor(kBlack);
0543     else          hits[nhits]->SetFillColor(kBlack);
0544     nhits++;
0545
0546     float phiHit=atan2(y[jj],x[jj]);
0547     if(abs(phiHit-phiSlice)<DphiSlice) {
0548       float Rhit =  sqrt(x[jj]*x[jj]+y[jj]*y[jj]);
0549       float Zhit =  a_si_z->at(jj);
0550
0551       if(a_layer->at(jj) < 3) { // MVTX hit
0552     hitsRZ[nhitsRZ]=new TBox(Zhit-s2,Rhit-s2,Zhit+s2,Rhit+s2);
0553     hitsRZ[nhitsRZ]->SetFillColor(kBlack);
0554       } else {// INTT hit
0555     if(abs(Zhit)<12.8) { // type-A strip
0556       hitsRZ[nhitsRZ]=new TBox(Zhit-sLA,Rhit-s2,Zhit+sLA,Rhit+s2);
0557     } else { // type-B strip
0558       hitsRZ[nhitsRZ]=new TBox(Zhit-sLB,Rhit-s2,Zhit+sLB,Rhit+s2);
0559     }
0560     hitsRZ[nhitsRZ]->SetFillColor(kBlack);
0561       }
0562       nhitsRZ++;
0563     }
0564   }
0565
0566   // Creates SiSeed tracklets (arrow)
0567   float Rtrk;
0568   const float Asize = 0.02;
0569   if(iframe==1) Rtrk=2.0; else Rtrk=80;
0570   if(iframe==1) Rmax=12.0; else Rmax=80;
0571   for(int k=0;k<ntrack;k++) {
0572     if(pt0->at(k)>pt0min) {
0573       //  cout<<k<<":"<<eta0->at(k)<<"  "<<phi0->at(k)<<"  "<<pt0->at(k)<<endl;
0574       float dx=Rtrk*cos(phi0->at(k));
0575       float dy=Rtrk*sin(phi0->at(k));
0576       seedtrk[ntrks]=new TArrow(vx0->at(k),vy0->at(k),vx0->at(k)+dx,vy0->at(k)+dy,Asize);
0577       ntrks++;
0578
0579       if(abs(phi0->at(k)-phiSlice)<DphiSlice) {
0580     float r0=sqrt((vx0->at(k))*(vx0->at(k))+(vy0->at(k))*(vy0->at(k)));
0581     float eta = eta0->at(k);
0582     float theta = 2*atan(exp(-eta));
0583     float Zend = (Rmax-r0)/tan(theta)+vz0->at(k);
0584
0585     seedtrkRZ[ntrksRZ]=new TArrow(vz0->at(k),r0,Zend,Rmax,Asize);
0586     ntrksRZ++;
0587     // if iframe==1, also plot "scaled" RZ track
0588     if(iframe==1){
0589       float Rmax2=80;
0590       float Zend2 = (Rmax2-r0)/tan(theta)+vz0->at(k);
0591       seedtrkRZ[ntrksRZ]=new TArrow(EMScale*(vz0->at(k)),
0592                     EMScale*r0,
0593                     EMScale*Zend2,
0594                     EMScale*Rmax2,Asize);
0595       seedtrkRZ[ntrksRZ]->SetLineStyle(2);
0596       ntrksRZ++;
0597     }
0598       }
0599     }
0600   }
0601
0602   // Create Orbit of EmcINTT tracklet if EmcINTT triplet exsits
0603   if(vEmcINTT.size()>0) {
0604     cout <<"EmcINTT tracklet found: # of tracklet = "<<vEmcINTT.size()<<endl;
0605     // clear previously existed track orbit;
0606     for(int io=0;io<vEmcINTTorbitL.size();io++) {
0607       vEmcINTTorbitL.at(io)->Delete();
0608     }
0609     vEmcINTTorbitL.clear();
0610
0611     for(int io=0;io<vEmcINTTorbitG.size();io++) {
0612       vEmcINTTorbitG.at(io)->Delete();
0613     }
0614     vEmcINTTorbitG.clear();
0615
0616     for(int i=0;i<vEmcINTT.size();i++) {
0617       TArc *orbitL = new TArc(vEmcINTT_xO.at(i),vEmcINTT_yO.at(i),
0618                  vEmcINTT_R.at(i),0.,360.);
0619       orbitL->SetFillStyle(0);
0620       orbitL->SetLineColor(kGreen);
0621       vEmcINTTorbitL.push_back(orbitL);
0622
0623       TArc *orbitG = new TArc(vEmcINTT_xO_EI.at(i),vEmcINTT_yO_EI.at(i),
0624                  vEmcINTT_R_EI.at(i),0.,360.);
0625       orbitG->SetFillStyle(0);
0626       orbitG->SetLineColor(kBlue);
0627       vEmcINTTorbitG.push_back(orbitG);
0628     }
0629   }
0630
0631   // Draw X-Y view
0632   c1->cd();
0633   for(int kk=0;kk<nhits;kk++){
0634     hits[kk]->Draw();
0635   }
0636   for(int itk=0;itk<ntrks;itk++) {
0637     seedtrk[itk]->Draw();
0638   }
0639   // Draw phi slice region
0640   if(dphiArc != nullptr) dphiArc->Delete();
0641   float phimin=(phiSlice-DphiSlice)*180./3.141592;
0642   float phimax=(phiSlice+DphiSlice)*180./3.141592;
0643   dphiArc = new TArc(0.,0.,Rmax,phimin,phimax);
0644   dphiArc->SetLineColorAlpha(0,0.1);
0645   dphiArc->SetFillColorAlpha(kBlue,0.1);
0646   dphiArc->Draw();
0647
0648   // Draw Circle for reference size of INTT and MVTX
0649   if(Circle2==nullptr) {
0650     Circle2=new TArc(0.,0.,2.3,0,360);
0651     Circle2->SetFillStyle(0);
0652   }
0653   if(Circle4==nullptr) {
0654     Circle4=new TArc(0.,0.,4.0,0,360);
0655     Circle4->SetFillStyle(0);
0656   }
0657   if(Circle8==nullptr) {
0658     Circle8=new TArc(0.,0.,7.6,0.,360);
0659     Circle8->SetFillStyle(0);
0660   }
0661   if(Circle10==nullptr) {
0662     Circle10=new TArc(0.,0.,10.2,0.,360);
0663     Circle10->SetFillStyle(0);
0664   }
0665   Circle2->Draw();
0666   Circle4->Draw();
0667   Circle8->Draw();
0668   Circle10->Draw();
0669
0670   for(int i=0;i<vEmcINTT.size();i++) {
0671     vEmcINTTorbitL.at(i)->Draw();
0672     vEmcINTTorbitG.at(i)->Draw();
0673   }
0674
0675   c1->Draw();
0676   c1->Modified();
0677   c1->Update();
0678   //  cout<<endl;
0679
0680   //Draw RZ slice
0681   c2->cd();
0682
0683   ZvtxMarker->Draw();
0684   for(int kk=0;kk<nhitsRZ;kk++){
0685     hitsRZ[kk]->Draw();
0686   }
0687   for(int itk=0;itk<ntrksRZ;itk++) {
0688     seedtrkRZ[itk]->Draw();
0689   }
0690   c2->Modified();
0691   c2->Update();
0692 }
0693
0694 //
0695 TH1F* hdphi1;
0696 TH1F* hdphi0;
0697 TH1F* hdeta1;
0698 TH1F* hdeta0;
0699 TH1F* hdphi1_0;
0700 TH1F* hdphi1t;
0701 TH1F* hdeta1t;
0702
0703 TH1F* hEMe;
0704 TH1F* hEMeM;
0705 TH1F* hEMeMt;
0706 TH2F* hEMC;
0707 TH2F* hEMCt;
0708 TH2F* hdphiM;
0709 TH2F* hdphiM2;
0710 TH2F* hdphiM3;
0711
0712 TH2F* hdphiMvtx_r;
0713 TH1F* hdphiMvtx1;
0714 TH1F* hdphiMvtx2;
0715 TH1F* hdphiMvtx3;
0716
0717 TH1F* hRorbit;
0718 TH1F* hpT;
0719
0720 TH1F* hzvtx;
0721 TH1F* hzvtxINTT;
0722
0723 const float xBC= -0.018; //  -0.018;
0724 const float yBC=0.126; //0.126;
0725
0726 void AnalyzeInit(void) {
0727   cout << "AnalyzeInit()"<<endl;
0728   hdphi1=new TH1F("hdphi1","dphi1",100,-0.5,0.5);
0729   hdphi0=new TH1F("hdphi0","dphi0",100,-0.5,0.5);
0730   hdeta1=new TH1F("hdeta1","eta1-EMCeta",150,-0.03,0.03);
0731   hdeta0=new TH1F("hdeta0","eta1-EMCeta",150,-0.03,0.03);
0732   hdphi1_0=new TH1F("hdphi1_0","phi1-phi0",100,-0.1,0.1);
0733   hdphi1t=new TH1F("hdphi1t","dphi1 (tight)",100,-0.5,0.5);
0734   hdeta1t=new TH1F("hdeta1t","eta1-EMCeta (tight)",150,-0.03,0.03);
0735
0736   hEMe = new TH1F("hEMe","EMC energy (inclusive)",250,0.,5.0);
0737   hEMeM = new TH1F("hEMeM","MEC energy (matched)",250,0.,5.0);
0738   hEMeMt = new TH1F("hEMeMt","MEC energy (tight matched)",250,0.,5.0);
0739   hEMC= new TH2F("hEMC","EMz vs EMphi matched w INTT",100,-150.,150.,100,-3.5,3.5);
0740   hEMCt= new TH2F("hEMCt","EMz vs EMphi tight matched w INTT",100,-150.,150.,100,-3.5,3.5);
0741   hdphiM = new TH2F("hdphiM","EMphi-INTT1phi vs INTT1phi-INTT0phi",100,-0.2,0.2,100,-0.025,0.025);
0742   hdphiM2 = new TH2F("hdphiM2","EMphi-INTT1phi vs INTT1phi-INTT0phi Rotate",100,-0.2,0.2,100,-0.025,0.025);
0743   hdphiM3 = new TH2F("hdphiM3","EMphi-INTT1phi vs INTT1phi-INTT0phi Rot, Tight",100,-0.2,0.2,100,-0.025,0.025);
0744
0745   hdphiMvtx_r = new TH2F("hdphiMvtx_r","dphi(MVTX) vs rMVTX",100,1.5,5.5,100,-0.05,0.05);
0746
0747   hdphiMvtx1 = new TH1F("hdphiMvtx1","dphi(MVTX) 2.5<r<3.5",100,-0.05,0.05);
0748   hdphiMvtx2 = new TH1F("hdphiMvtx2","dphi(MVTX) 3.5<r<4.3",100,-0.05,0.05);
0749   hdphiMvtx3 = new TH1F("hdphiMvtx3","dphi(MVTX) 4.3<r<5.5",100,-0.05,0.05);
0750
0751   hRorbit   = new TH1F("hRorbit","track orbit radius",200,0.,400);
0752   hpT       = new TH1F("hpT","track pT",100,0.,5.0);
0753   hzvtx     = new TH1F("hzvtx","zvtx in evtTree",41,-20.5,20.5);
0754   hzvtxINTT = new TH1F("hzvtxINTT","zvtx by INTT",41,-20.5,20.5);
0755
0756 }
0757
0758 void AnalyzeEvent(void) {
0759   float zvtxINTT;
0760   INTT_zvtx_Finder(zvtxINTT);
0761   LoadTrackerData();
0762   FindEmcINTTtriplet();
0763   FindEmcINTT_Mvtx();
0764
0765   hzvtxINTT->Fill(zvtxINTT);
0766 }
0767
0768 void AnalyzeEnd(void) {
0769   gStyle->SetOptFit(1111);
0770   TF1* f1= new TF1("f1","gaus(0)+pol0(3)",-0.03,0.03);
0771
0772   TCanvas *c3=new TCanvas("c3");
0773   hdphi1->SetMinimum(0);
0774   hdphi1->Draw();
0775
0776   TCanvas *c4=new TCanvas("c4");
0777   hdphi0->SetMinimum(0);
0778   hdphi0->Draw();
0779
0780   TCanvas *c5=new TCanvas("c5");
0781   hdeta1->SetMinimum(0);
0782   hdeta1->Draw();
0783
0784   TCanvas *c6=new TCanvas("c6");
0785   hdeta0->SetMinimum(0);
0786   hdeta0->Draw();
0787
0788   TCanvas *c7=new TCanvas("c7");
0789   f1->SetParameters(300,0.,0.01,0.);
0790   hdphi1_0->Draw();
0791   hdphi1_0->Fit("f1");
0792
0793   TCanvas *c8=new TCanvas("c8");
0794   hEMC->Draw();
0795
0796   TCanvas *c9=new TCanvas("c9");
0797   hdphiM->Draw();
0798
0799   TCanvas *c10=new TCanvas("c10");
0800   hdphiM2->Draw();
0801
0802   TCanvas *c11=new TCanvas("c11");
0803   hdphiM3->Draw();
0804
0805   TCanvas *c12=new TCanvas("c12");
0806   hEMCt->Draw();
0807
0808   TCanvas *c13=new TCanvas("c13");
0809   hdphi1t->SetMinimum(0);
0810   hdphi1t->Draw();
0811
0812   TCanvas *c14=new TCanvas("c14");
0813   hdeta1t->SetMinimum(0);
0814   hdeta1t->Draw();
0815
0816   TCanvas *c15=new TCanvas("c15");
0817   hdphiMvtx_r->Draw();
0818
0819   TCanvas *c16=new TCanvas("c16");
0820   f1->SetParameters(600,0.,0.007,50);
0821   hdphiMvtx1->Draw();
0822   hdphiMvtx1->Fit("f1","","",-0.03,0.03);
0823
0824   TCanvas *c17=new TCanvas("c17");
0825   f1->SetParameters(700,0.,0.005,50);
0826   hdphiMvtx2->Draw();
0827   hdphiMvtx2->Fit("f1","","",-0.03,0.03);
0828
0829   TCanvas *c18=new TCanvas("c18");
0830   f1->SetParameters(500,0.,0.003,50);
0831   hdphiMvtx3->Draw();
0832   hdphiMvtx3->Fit("f1","","",-0.03,0.03);
0833
0834   TFile *file = new TFile("AnaEvent.root","RECREATE");
0835   hdphi1->Write();
0836   hdphi0->Write();
0837   hdeta1->Write();
0838   hdeta0->Write();
0839   hdphi1_0->Write();
0840   hdphi1t->Write();
0841   hdeta1t->Write();
0842   hEMC->Write();
0843   hEMCt->Write();
0844   hEMe->Write();
0845   hEMeM->Write();
0846   hEMeMt->Write();
0847   hdphiM->Write();
0848   hdphiM2->Write();
0849   hdphiM3->Write();
0850   hdphiMvtx_r->Write();
0851   hdphiMvtx1->Write();
0852   hdphiMvtx2->Write();
0853   hdphiMvtx3->Write();
0854   hRorbit->Write();
0855   hpT->Write();
0856
0857   hzvtx->Write();
0858   hzvtxINTT->Write();
0859
0860   file->Close();
0861 }
0862
0863 void LoadTrackerData(void) {
0864
0865   int nEmc = energy->size();
0866   int nSiA  = a_layer->size();
0867
0868   //  cout << "nEmc = "<<nEmc <<endl;
0869   //  cout << "nSiA = "<<nSiA <<endl;
0870
0871   vEMr.clear();
0872   vEMphi.clear();
0873   vEMz.clear();
0874   vEMe.clear();
0875
0876   vMVTXr.clear();
0877   vMVTXphi.clear();
0878   vMVTXx.clear();
0879   vMVTXy.clear();
0880   vMVTXz.clear();
0881
0882   vINTT0r.clear();
0883   vINTT0phi.clear();
0884   vINTT0z.clear();
0885
0886   vINTT1r.clear();
0887   vINTT1phi.clear();
0888   vINTT1z.clear();
0889
0890   //Fill vector of EMC data
0891   for(int iemc=0;iemc<nEmc;iemc++) {
0892     if(!Is_Hot_EMC(iemc)){
0893       vEMr.push_back(emc_r->at(iemc));
0894       vEMphi.push_back(emc_phi->at(iemc));
0895       vEMz.push_back(emc_z->at(iemc));
0896       vEMe.push_back(energy->at(iemc));
0897     }
0898   }
0899
0900   //Fill vector of MVTX and INTT data
0901   for(int i=0;i<nSiA;i++) {
0902     float xhit = a_si_x->at(i)-xBC;
0903     float yhit = a_si_y->at(i)-yBC;
0904     float zhit = a_si_z->at(i);
0905     float phi_hit = atan2(yhit,xhit);
0906     float rhit = sqrt(xhit*xhit+yhit*yhit);
0907     if(a_layer->at(i)<3) {       // MVTX
0908       vMVTXr.push_back(rhit);
0909       vMVTXphi.push_back(phi_hit);
0910       vMVTXx.push_back(xhit);
0911       vMVTXy.push_back(yhit);
0912       vMVTXz.push_back(zhit);
0913     }
0914     else if(a_layer->at(i)< 5) { //INTT0
0915       vINTT0r.push_back(rhit);
0916       vINTT0phi.push_back(phi_hit);
0917       vINTT0z.push_back(zhit);
0918     }
0919     else {                       //INTT1
0920       vINTT1r.push_back(rhit);
0921       vINTT1phi.push_back(phi_hit);
0922       vINTT1z.push_back(zhit);
0923     } //if(a_layer)
0924   }// for(i)
0925   //  cout << " nMVTX="<<vMVTXr.size();
0926   //  cout << " nINTT0=" <<vINTT0r.size()<< " nINTT1="<<vINTT1r.size()<<endl;
0927
0928 }
0929
0930 float rz2eta(float r,float z) {
0931   if(z>=0) {
0932     return log((sqrt(r*r+z*z)+z)/r);
0933   } else {
0934     return -log((sqrt(r*r+z*z)-z)/r);
0935   }
0936 }
0937
0938
0939 void FindEmcINTTtriplet(void) {
0940   // constants for cuts
0941   const float dPhi = 0.5;         //dphi from EMC hit to search INTT hit
0942   const float deta_match = 0.03;  //0.03 deta cut between EMC hit and INTT hit
0943   const float deta_match2 = 0.01; //tighter deta cut
0944   const float dphi1_0match = 0.02;   // acceptance range of INTT1phi - INTT0phi
0945   const float dphi1_0match2 = 0.007; // tight cut (after rotational correction)
0946   const float Krot = 0.025;          //correction factor for EMdphi vs INTTdphi
0947
0948   vector<int> INTT0m;
0949   vector<int> INTT1m;
0950   vector<float> INTT0phim;
0951   vector<float> INTT1phim;
0952   vector<float> INTT0detam;
0953   vector<float> INTT1detam;
0954   vector<float> INTT0dphim;
0955   vector<float> INTT1dphim;
0956
0957   hzvtx->Fill(zvtx);
0958
0959   vEmcINTT.clear();
0960   vEmcINTT_Mvtx1.clear();
0961   vEmcINTT_Mvtx2.clear();
0962   vEmcINTT_Mvtx3.clear();
0963   vEmcINTT_R.clear();
0964   vEmcINTT_xO.clear();
0965   vEmcINTT_yO.clear();
0966   vEmcINTT_R_EI.clear();
0967   vEmcINTT_xO_EI.clear();
0968   vEmcINTT_yO_EI.clear();
0969
0970   // search for INTT hits that matches with EMC in phi and eta
0971   for(int iemc=0;iemc<vEMphi.size();iemc++) {
0972     INTT0m.clear();
0973     INTT1m.clear();
0974     INTT0phim.clear();
0975     INTT1phim.clear();
0976     INTT0detam.clear();
0977     INTT1detam.clear();
0978     INTT0dphim.clear();
0979     INTT1dphim.clear();
0980
0981     hEMe->Fill(vEMe.at(iemc));
0982     float EMCeta = rz2eta(vEMr.at(iemc),vEMz.at(iemc)-zvtx);
0983
0984     // search for matched INTT1 hits
0985     int n1=vINTT1phi.size();
0986     for(int i1=0;i1<n1;i1++){
0987       float dphi1=vINTT1phi.at(i1)-vEMphi.at(iemc);
0988       float INTT1eta = rz2eta(vEMr.at(iemc)-vINTT1r.at(i1),vEMz.at(iemc)-vINTT1z.at(i1));
0989       float deta1 = INTT1eta - EMCeta;
0990       if(abs(dphi1)<dPhi && abs(deta1)<deta_match) {
0991     INTT1m.push_back(i1);
0992     INTT1phim.push_back(vINTT1phi.at(i1));
0993     INTT1detam.push_back(deta1);
0994     INTT1dphim.push_back(dphi1);
0995       }
0996     }
0997
0998     // search for matched INTT0 hits
0999     int n0=vINTT0phi.size();
1000     for(int i0=0;i0<n0;i0++){
1001       float dphi0=vINTT0phi.at(i0)-vEMphi.at(iemc);
1002       float INTT0eta = rz2eta(vEMr.at(iemc)-vINTT0r.at(i0),vEMz.at(iemc)-vINTT0z.at(i0));
1003       float deta0 = INTT0eta - EMCeta;
1004       if(abs(dphi0)<dPhi && abs(deta0)<deta_match) {
1005     INTT0m.push_back(i0);
1006     INTT0phim.push_back(vINTT0phi.at(i0));
1007     INTT0detam.push_back(deta0);
1008     INTT0dphim.push_back(dphi0);
1009       }
1010     }
1011
1012     // if there are phi-mathced INTT hit pairs
1013     if(INTT0m.size()>0 && INTT1m.size()>0) {
1014       for(int i1=0;i1<INTT1m.size();i1++) {
1015     for(int i0=0;i0<INTT0m.size();i0++) {
1016       float dphi1_0 = INTT1phim.at(i1)-INTT0phim.at(i0);
1017       hdphi1_0->Fill(dphi1_0);
1018       if(abs(dphi1_0)<0.05) {
1019         float dphiEM_1=vEMphi.at(iemc)-INTT1phim.at(i1);
1020         hdphiM->Fill(dphiEM_1,dphi1_0);
1021         hdphiM2->Fill(dphiEM_1,dphi1_0-Krot*dphiEM_1);
1022
1023         if(abs(dphi1_0-0.03*dphiEM_1)<dphi1_0match2) {//matched!
1024           hdphi0->Fill(INTT0dphim.at(i0));
1025           hdphi1->Fill(INTT1dphim.at(i1));
1026           hdeta0->Fill(INTT0detam.at(i0));
1027           hdeta1->Fill(INTT1detam.at(i1));
1028           hEMC->Fill(vEMz.at(iemc),vEMphi.at(iemc));
1029               hEMeM->Fill(vEMe.at(iemc));
1030         }
1031         // apply tighter deta cut
1032         if(abs(INTT0detam.at(i0))<deta_match2 &&
1033            abs(INTT1detam.at(i1))<deta_match2) {
1034           hdphiM3->Fill(dphiEM_1,dphi1_0-Krot*dphiEM_1);
1035           //apply tighter dphi cut
1036           if(abs(dphi1_0-Krot*dphiEM_1)<dphi1_0match2) {
1037         hdeta1t->Fill(INTT1detam.at(i1));
1038         hdphi1t->Fill(INTT1dphim.at(i1));
1039         hEMCt->Fill(vEMz.at(iemc),vEMphi.at(iemc));
1040         hEMeMt->Fill(vEMe.at(iemc));
1041         // a triplet is found. Store it.
1042         EmcINTTtriplet triplet;
1043         triplet.iEmc = iemc;
1044         triplet.iINTT0 = INTT0m.at(i0);
1045         triplet.iINTT1 = INTT1m.at(i1);
1046         vEmcINTT.push_back(triplet);
1047         if(imode == 2) {
1048           cout <<"one triplet found:"<<iemc<<" "<<INTT0m.at(i0)<< " "<<INTT1m.at(i1)<<endl;
1049         }
1050           }
1051         }
1052       }
1053     }//for(i0)
1054       }//for(i1)
1055     }//if
1056   }// iemc
1057 }
1058
1059 void FindEmcINTT_Mvtx(void) {
1060   // Search for MVTX hits that matches EmcINTTtriplet, and add the
1061   // matched hits
1062   if(imode == 2) {
1063     cout << "# of triplet ="<<vEmcINTT.size()<<endl;
1064   }
1065
1066   const float dphi_MVTX_INTT0 = 0.05;  // rough angle cut
1067   const float dphiMvtx1 = 0.016;
1068   const float dphiMvtx2 = 0.012;
1069   const float dphiMvtx3 = 0.009;
1070
1071   for(int i=0;i<vEmcINTT.size();i++){
1072     vEmcINTT_Mvtx1.push_back(-1);
1073     vEmcINTT_Mvtx2.push_back(-1);
1074     vEmcINTT_Mvtx3.push_back(-1);
1075     vEmcINTT_R.push_back(0);
1076     vEmcINTT_xO.push_back(0);
1077     vEmcINTT_yO.push_back(0);
1078     vEmcINTT_R_EI.push_back(0);
1079     vEmcINTT_xO_EI.push_back(0);
1080     vEmcINTT_yO_EI.push_back(0);
1081   }
1082
1083   const float r_ave = 6.0; // average r of Silicon tracker
1084   for(int i=0;i<vEmcINTT.size();i++) {
1085     float phi0 =vINTT0phi.at((vEmcINTT.at(i)).iINTT0);
1086     float r0   =vINTT0r.at((vEmcINTT.at(i)).iINTT0);
1087     float phi1 =vINTT1phi.at((vEmcINTT.at(i)).iINTT1);
1088     float r1   =vINTT1r.at((vEmcINTT.at(i)).iINTT1);
1089     float rEmc =vEMr.at((vEmcINTT.at(i)).iEmc);
1090     float zEmc =vEMz.at((vEmcINTT.at(i)).iEmc);
1091
1092     float a=(phi1-phi0)/(r1*r1-r0*r0);
1093     float R_orbit = 1/sqrt(abs(2*a));
1094     float phi00 = phi0 - a*r0*r0;
1095     hRorbit->Fill(R_orbit);
1096     hpT->Fill(0.42*R_orbit/100);
1097
1098     //Local orbit. There parameters are good in Si tracker region
1099     float phi00c = phi0 - r_ave/R_orbit;
1100     float phiO = phi00c+PI/2.0;
1101     float xO = R_orbit*cos(phiO);
1102     float yO = R_orbit*sin(phiO);
1103     // store parameters of local orbit
1104     vEmcINTT_R.at(i)=R_orbit;
1105     vEmcINTT_xO.at(i)=xO;
1106     vEmcINTT_yO.at(i)=yO;
1107     cout <<"R(local)="<<R_orbit<<endl;
1108
1109     // Global orbit. A circule that go through
1110     // INTT0, INTT1, Emc
1111     {
1112       int iEmc = (vEmcINTT.at(i)).iEmc;
1113       int iINTT0 = (vEmcINTT.at(i)).iINTT0;
1114       int iINTT1 = (vEmcINTT.at(i)).iINTT1;
1115
1116       float r2   = vEMr.at(iEmc);
1117       float phi2 = vEMphi.at(iEmc);
1118       float x2   = r2*cos(phi2);
1119       float y2   = r2*sin(phi2);
1120       //      cout <<"r2="<<r2<<" phi2="<<phi2<<endl;
1121
1122       float r1   = vINTT1r.at(iINTT1);
1123       float phi1 = vINTT1phi.at(iINTT1);
1124       float x1   = r1*cos(phi1);
1125       float y1   = r1*sin(phi1);
1126       //      cout <<"r1="<<r1<<" phi1="<<phi1<<endl;
1127
1128       float r0   = vINTT0r.at(iINTT0);
1129       float phi0 = vINTT0phi.at(iINTT0);
1130       float x0   = r0*cos(phi0);
1131       float y0   = r0*sin(phi0);
1132       //      cout <<"r0="<<r0<<" phi0="<<phi0<<endl;
1133
1134       float r12 = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
1135       float r01 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
1136       float R_EmcINTT = 0.5*(x2-x0)/((y2-y1)/r12 - (y1-y0)/r01);
1137
1138       float xO_EmcINTT = 0.5*(x1+x0)-R_EmcINTT*(y1-y0)/r01;
1139       float yO_EmcINTT = 0.5*(y1+y0)+R_EmcINTT*(x1-x0)/r01;
1140
1141       vEmcINTT_R_EI.at(i)=abs(R_EmcINTT);
1142       vEmcINTT_xO_EI.at(i)=xO_EmcINTT;
1143       vEmcINTT_yO_EI.at(i)=yO_EmcINTT;
1144
1145       cout << "R_EI = "<<abs(R_EmcINTT)<<endl;
1146     }
1147
1148     //
1149     //
1150     // find the closest Mvtx hits to the tracklet orbit
1151     // for each of the three layers of Mvtx
1152     //
1153     float Mvtx1m_dphi=0.1;
1154     float Mvtx2m_dphi=0.1;
1155     float Mvtx3m_dphi=0.1;
1156     for(int iMvtx=0;iMvtx<vMVTXr.size();iMvtx++) {// vMVTX
1157       float MVTXphi = vMVTXphi.at(iMvtx);
1158       float xMvtx = vMVTXx.at(iMvtx) - xBC;
1159       float yMvtx = vMVTXy.at(iMvtx) - yBC;
1160       float rMvtx = sqrt(xMvtx*xMvtx+yMvtx*yMvtx);
1161
1162       // Mvtx1 (inner-most layer)
1163       if( abs(MVTXphi - phi0) < dphi_MVTX_INTT0) {
1164     float dphiMvtx = MVTXphi-phi00-a*rMvtx*rMvtx;
1165     hdphiMvtx_r->Fill(rMvtx,dphiMvtx);
1166     if(2.5< rMvtx && rMvtx <3.5) {
1167       hdphiMvtx1->Fill(dphiMvtx);
1168       if(abs(dphiMvtx)<dphiMvtx1) {
1169         if(abs(dphiMvtx)<abs(Mvtx1m_dphi)) {
1170           Mvtx1m_dphi = dphiMvtx;
1171           vEmcINTT_Mvtx1.at(i)=iMvtx;
1172         }
1173       }
1174     }
1175     // Mvtx2 (middle layer)
1176     if(3.5< rMvtx && rMvtx <4.3) {
1177       hdphiMvtx2->Fill(dphiMvtx);
1178       if(abs(dphiMvtx)<dphiMvtx2) {
1179         if(abs(dphiMvtx)<abs(Mvtx2m_dphi)) {
1180           Mvtx2m_dphi = dphiMvtx;
1181           vEmcINTT_Mvtx2.at(i)=iMvtx;
1182         }
1183       }
1184     }
1185     // Mvtx3 (outer-most layer)
1186     if(4.3< rMvtx && rMvtx <5.5) {
1187       hdphiMvtx3->Fill(dphiMvtx);
1188       if(abs(dphiMvtx)<abs(Mvtx3m_dphi)) {
1189         Mvtx3m_dphi = dphiMvtx;
1190         vEmcINTT_Mvtx3.at(i)=iMvtx;
1191       }
1192     }
1193       }
1194     }// for(iMVTX);
1195     //
1196     // Now, for each of vEmcINTT triplet, Mvtx hits closest to the track
1197     // orbit are assocated
1198     //
1199     if(imode==2) {
1200       cout << "EMphi = " <<vEMphi.at((vEmcINTT.at(i)).iEmc);
1201       cout <<" INTT1phi= "<<vINTT1phi.at((vEmcINTT.at(i)).iINTT1);
1202       cout <<" INTT0phi= "<<vINTT0phi.at((vEmcINTT.at(i)).iINTT0)<<endl;
1203       cout <<" R_orbit = "<< R_orbit;
1204       cout <<" phi00 = "<< phi00 << endl;
1205       cout << "xO="<<xO<<" yO="<<yO<<" phiO="<<phiO<<endl;
1206       cout << " Mvtx1 = "<<vEmcINTT_Mvtx1.at(i)<<" dphi= "<<Mvtx1m_dphi;
1207       cout << " Mvtx2 = "<<vEmcINTT_Mvtx2.at(i)<<" dphi= "<<Mvtx2m_dphi;
1208       cout << " Mvtx3 = "<<vEmcINTT_Mvtx3.at(i)<<" dphi= "<<Mvtx3m_dphi;
1209       cout <<endl;
1210     }
1211   }// for(i)   loop over vEmcINTT
1212 }
1213 void INTT_zvtx_Finder(float &zvtxINTT){
1214   // 2026/05/23
1215   // 1) form pair of INTT0 hit and INTT1 hit that are close in phi
1216   // 2) assume zvtx position, and calculate the eta range of the hit strip
1217   //    of INTT1 and INTT0 pairs
1218   // 3) calculate the overlap of the eta range of the pairs
1219   // 4) change zvtx and maximize the sum of the eta overlap of the pairs.
1220   // 2026/05/24
1221   // Test results shows that this method doesnt work reasonaly
1222   // well. But the Zvtx resolution is about 1cm or so.
1223   // For more precise ztx, need to combine with MVTX.
1224
1225   const float dphiINTTmax=0.03;
1226   const float Zmax_A = 12.8;  //Max abs(Z) value of A-type INTT (cm)
1227   const float Zhalf_A = 0.8;  // half length of A-type strip (cm)
1228   const float Zhalf_B = 1.0;  // half lenght of B-type strip (cm)
1229
1230
1231   // INTT hits
1232   vector<int>   INTT1hit;  //index of INTT1 hit in a_si_x,a_si_y,a_si_z
1233   vector<float> INTT1phi;  //phi value of INTT1 hit
1234   vector<int>   INTT0hit;  //index of INTT0 hit in a_si_x,a_si_y,a_si_z
1235   vector<float> INTT0phi;  //phi value of INTT1 hit
1236
1237   // pair of (INTT1hit,INTT0hit)
1238   vector<int>   INTTpair1; //index of the pair in INTT1hit
1239   vector<int>   INTTpair0; //index of the pair in INTT0hit
1240   //
1241   vector<float> pair1eta_max;  //max eta of INTT1hit of the pair
1242   vector<float> pair1eta_min;  //min eta of INTT1hit of the pair
1243   vector<float> pair0eta_max;  //max eta of INTT0hit of the pair
1244   vector<float> pair0eta_min;  //min eta of INTT0hit of the pair
1245
1246   //form INTT1hit and INTT0hit
1247   INTT1hit.clear();
1248   INTT0hit.clear();
1249   INTT1phi.clear();
1250   INTT0phi.clear();
1251   for(int i=0;i<a_layer->size();i++) {
1252     if(a_layer->at(i)==3 || a_layer->at(i)==4) {//INTT0
1253       INTT0hit.push_back(i);
1254       INTT0phi.push_back(atan2(a_si_y->at(i),a_si_x->at(i)));
1255     }
1256     if(a_layer->at(i)==5 || a_layer->at(i)==6) {//INTT1
1257       INTT1hit.push_back(i);
1258       INTT1phi.push_back(atan2(a_si_y->at(i),a_si_x->at(i)));
1259     }
1260   }
1261
1262   //form pair of INTT1 and INTT0 hit
1263   INTTpair1.clear();
1264   INTTpair0.clear();
1265   for(int i1=0;i1<INTT1hit.size();i1++) {
1266     for(int i0=0;i0<INTT0hit.size();i0++) {
1267       if(abs(INTT0phi.at(i0)-INTT1phi.at(i1))<dphiINTTmax) {
1268     INTTpair1.push_back(i1);
1269     INTTpair0.push_back(i0);
1270       }
1271     }
1272   }
1273
1274   // print the formed pairs
1275   /*
1276   cout << "INTT pairs:"<<endl;
1277   for(int ipair=0;ipair<INTTpair0.size();ipair++) {
1278     cout << ipair <<": ("<<INTTpair1.at(ipair)<<","<<INTTpair0.at(ipair);
1279     cout <<")";
1280     cout <<"   ("<<INTT1phi.at(INTTpair1.at(ipair))<<",";
1281     cout <<"   ("<<INTT0phi.at(INTTpair0.at(ipair))<<")"<< endl;
1282   }
1283   */
1284
1285   // slide zvtx from -20 to 20. For a given zvtx, calculate the
1286   // sum of (eta overlap) of the pairs.
1287   float max_overlap = 0.;
1288   float zvtx_max_overlap;
1289   float zvtx;
1290   for(zvtx=-20.0;zvtx<20.1;zvtx+=1) {
1291     float sum_overlap = 0.;
1292     for(int ipair=0;ipair<INTTpair0.size();ipair++) {
1293       float z1=a_si_z->at(INTT1hit.at(INTTpair1.at(ipair)));
1294       float z0=a_si_z->at(INTT0hit.at(INTTpair0.at(ipair)));
1295       float z1min, z1max;
1296       float z0min, z0max;
1297       if(abs(z1)<Zmax_A) {  //A-type INTT sensor
1298     z1min=z1-Zhalf_A;
1299     z1max=z1+Zhalf_A;
1300       } else {              //B-type INTT sensor
1301     z1min=z1-Zhalf_B;
1302     z1max=z1+Zhalf_B;
1303       }
1304       if(abs(z0)<Zmax_A) {  //A-type INTT sensor
1305     z0min=z0-Zhalf_A;
1306     z0max=z0+Zhalf_A;
1307       } else {              //B-type INTT sensor
1308     z0min=z1-Zhalf_B;
1309     z0max=z1+Zhalf_B;
1310       }
1311       //convert z to eta
1312       float x1=a_si_x->at(INTT1hit.at(INTTpair1.at(ipair)));
1313       float y1=a_si_y->at(INTT1hit.at(INTTpair1.at(ipair)));
1314       float r1=sqrt(x1*x1+y1*y1);
1315       float x0=a_si_x->at(INTT0hit.at(INTTpair0.at(ipair)));
1316       float y0=a_si_y->at(INTT0hit.at(INTTpair0.at(ipair)));
1317       float r0=sqrt(x0*x0+y0*y0);
1318       float eta1max=rz2eta(r1,z1max-zvtx);
1319       float eta1min=rz2eta(r1,z1min-zvtx);
1320       float eta0max=rz2eta(r0,z0max-zvtx);
1321       float eta0min=rz2eta(r0,z0min-zvtx);
1322
1323       float eta_overlap=0;
1324       if(eta1max > eta0max && eta0max > eta1min) {
1325     eta_overlap = eta0max-eta1min;
1326       } else if(eta0max > eta1max && eta1max > eta0min) {
1327     eta_overlap = eta1max - eta0min;
1328       }
1329       //      cout << ipair <<": ("<< eta1max <<","<< eta1min<<")";
1330       //      cout << " ("<< eta0max <<","<< eta0min<<")";
1331       //      cout << " overlap = "<< eta_overlap << endl;
1332       sum_overlap += eta_overlap;
1333       if(sum_overlap > max_overlap) {
1334     max_overlap = sum_overlap;
1335     zvtx_max_overlap = zvtx;
1336       }
1337     }//for(ipair)
1338     //    cout << "zvtx="<<zvtx<<": sum_overlap = "<<sum_overlap;
1339     //    cout << "  max_overlap="<<max_overlap<<endl;
1340   }//for(zvtx)
1341   zvtxINTT = zvtx_max_overlap;
1342 }