StRoot: StEmcPool/StBarrelCalib2008/BarrelMipCalib.cxx Source File

00001 #include <TH1F.h>
00002 #include <TH2F.h>
00003 #include <TObjArray.h>
00004 #include <TGraph.h>
00005 #include <TVector2.h>
00006 
00007 #include "BarrelMipCalib.h"
00008 #include "StJanBarrelDbMaker.h"
00009 #include "JanBarrelEvent.h"
00010 
00011 #include "StEmcUtil/geometry/StEmcGeom.h"
00012 
00013 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00014 #include "StMuDSTMaker/COMMON/StMuDst.h"
00015 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00016 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
00017 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00018 #include "StThreeVectorF.hh"
00019 #include "StThreeVectorD.hh"
00020 
00021 
00022 //________________________________________________
00023 //________________________________________________
00024 BarrelMipCalib::BarrelMipCalib( TObjArray *HList, StJanBarrelDbMaker* jdb,  StMuDstMaker *mumk){
00025     mJanDbMaker=jdb; assert(mJanDbMaker);
00026     muMaker=mumk;  assert(muMaker);
00027    
00028     assert(HList);
00029     memset(hA,0,sizeof(hA));
00030 
00031 
00032     hA[0]=new TH1F("mipStat","Mip calib statistics; case",30,-0.5,29.5);
00033     hA[1]=new TH1F("mipZver","primary (smalest) vertez Z ; vertx Z (cm)",200,-200,200);
00034     hA[2]=new TH1F("mipTrPt","primTr pT; pT (GeV/c)",100,0,20);
00035     hA[3]=new TH1F("mipTrEta","primTr Eta ; pseudorapidity",100,-2,2);
00036     hA[4]=new TH1F("mipTrNff","primTr nFit/nPos ; nFit/nPos",20,0,1.1);
00037 
00038     int nx=600, ny=720;
00039     hA[5]=new TH2F("mipTrRZ1","primTr punch BPRS  ; Z @ BPRS (cm); phi @ BPRS (rad)",nx,-300,300,ny,0,2*C_PI);
00040     hA[6]=new TH2F("mipTrRZ21","primTr punch BTOW L=21, input; Z @ BTOW layer=21 (cm); phi @ BPRS (rad)",nx,-300,300,ny,0,2*C_PI);
00041     hA[7]=new TH2F("mipTrRZf","primTr in BPRS fiducial ; Z @ BPRS (cm); phi @ BPRS (rad)",nx,-300,300,ny,0,2*C_PI);
00042 
00043     hA[8]=new TH2F("mipDeDx","primTr TPC deDx; track momentum (GeV); TPC dEdX (keV)", 100,0,5,100,0, 20);
00044 
00045     hA[9]=new TH1F("rankZver","highest rank  primary vertez Z ; vertx Z (cm)",200,-200,200);
00046 
00047     hA[10]=new TH1F("mipRxyTr","primTr Rxy of last point on the track; Rxy (cm)",220,0,220);
00048     int nb=200;
00049     float adc1=-50;
00050     float adc2=adc1+nb;
00051  
00052     //free 11
00053     hA[12]=new TH1F("mipTrTw","tower pointed by accepted prim track; BPRS softID",mxBtow,0.5, mxBtow+0.5);  
00054     // accepted events plots
00055     hA[13]=new TH1F("mipZverAc","primary vertez Z , accepted; vertx Z (cm)",200,-200,200);
00056     hA[14]=new TH1F("mipTrPtAc","primTr pT, accepted; pT (GeV/c)",100,0,20);
00057     hA[15]=new TH1F("mipTrPAc","primTr p, accepted; momentum (GeV)",100,0,20);
00058 
00059     // disable to de-activate swapSkan
00060     // hA[16]=new TH2F("swapScan","BTOW tile pointed by MIP vs. tiles with MIP-like ADC ; softID if ADC=MIP;  ID of track-tower",mxBtow,0.5, mxBtow+0.5,mxBtow,0.5, mxBtow+0.5);
00061     // 19 free
00062     hA[17]=new TH1F("mipTrZf","primTr in BPRS fiducial ; Z @ BPRS (cm)",1800,-300,300);
00063     hA[18]=new TH1F("mipTrRf","primTr in BPRS fiducial , z>0; phi @ BPRS (rad)",3600,0,2*C_PI);
00064 
00065     hA[20]=new TH1F("mipMaZ1","Z distance of track from both tower Z-edges;  delZ (cm)", 250,-5,20);
00066     hA[21]=new TH1F("mipMaP1","R*phi distance of track from tower edge;  delR*phi (cm)", 200,-10,10);
00067 
00068     hA[22]=new TH2F("mipBprsTr","ADC for BPRS pointed by TPC MIP track; BPRS softID; rawADC-ped",mxBtow,0.5,mxBtow+0.5, nb,adc1,adc2);
00069     hA[23]=new TH2F("mipBprsTrBt","ADC for BPRS pointed by TPC MIP track + BTOW=MIP; BPRS softID; rawADC-ped",mxBtow,0.5,mxBtow+0.5, nb,adc1,adc2);
00070 
00071     hA[24]=new TH2F("mipBtowTr","ADC for BTOW pointed by TPC track; BTOW softID; rawADC-ped",mxBtow,0.5,mxBtow+0.5, nb,adc1,adc2);
00072     hA[25]=new TH2F("mipBtowTrPr","ADC for BTOW pointed by TPC track+BPRS=MIP; BTOW softID; rawADC-ped",mxBtow,0.5,mxBtow+0.5, nb,adc1,adc2);
00073 
00074 
00075 
00076     for(int i=0;i<mxH;i++)
00077       if(hA[i]) HList->Add(hA[i]);
00078 }
00079 
00080 
00081 
00082 
00083 //________________________________________________
00084 //________________________________________________
00085 void BarrelMipCalib::search(  JanBarrelEvent &fullEve){
00086 
00087   // Access to muDst .......................
00088   StMuEvent* muEve = muMaker->muDst()->event();
00089   int nPrimV=muMaker->muDst()->numberOfPrimaryVertices();
00090   hA[0]->Fill(0);
00091   
00092 #if 0
00093   StMuTriggerIdCollection &tic=muEve->triggerIdCollection();
00094   int trigID=96211;
00095   bool fired=tic.nominal().isTrigger(trigID);
00096 #endif
00097   
00098   printf("\nmipCalib eventID %d nPrimV=%d  =============\n", muEve->eventId(),nPrimV);
00099   if(nPrimV<=0) return;
00100   hA[0]->Fill(1);
00101 
00102   
00103   int nPrimTr =0;
00104   
00105   // find vertex with smalest Z - to compensate for PPV error in FMS-production
00106 
00107   int iVert=-1;
00108   float zVert=99999;
00109   
00110   for(int iv=0;iv<nPrimV;iv++) {
00111     StMuPrimaryVertex* V= muMaker->muDst()->primaryVertex(iv);
00112     assert(V);
00113     if(V->ranking()<0.) continue; //drop vertices with too few matched tracks
00114     muMaker->muDst()->setVertexIndex(iv);
00115     const StThreeVectorF &r=V->position();
00116     const StThreeVectorF &er=V->posError();
00117     if(iv==0) hA[9]->Fill(r.z()); // histo only highe rank vertex
00118 #if 0
00119     cout << "\nPrimary track " << nPrimTr << " momentum " << pr_track->p() << endl;  cout << "\t flag=" << pr_track->flag() << " nHits=" << pr_track->nHits()<< " vertID="<<  pr_track->vertexIndex()<< endl;
00120     cout << "\t primV("<<iv<<")  primDCA=" << pr_track->dca(iv) << ", pT="<< pr_track->pt()<< endl;
00121     if(pr_track->dca(iv).mag()>5)       cout << "^^^^^ 3D DCA magnitude="<<pr_track->dca(iv).mag()<<endl;
00122     cout << "\t first point " << pr_track->firstPoint() << endl;
00123     cout << "\t last point " << pr_track->lastPoint() << endl;
00124 #endif
00125       
00126 
00127     printf("iv=%d   Vz=%.2f +/-%.2f  bestZ=%f\n",iv,r.z(),er.z() ,zVert );
00128   
00129     if( fabs(zVert)< fabs(r.z()))  continue;
00130     zVert=r.z(); 
00131     iVert=iv;
00132   }
00133 
00134   //  printf("iVert=%d\n",iVert);
00135   if(iVert<0) return; // non of vertices has positive rank
00136 
00137   hA[1]->Fill(zVert);
00138   if(fabs(zVert) > cut_zVertex) return; // drop if too large Z 
00139   hA[0]->Fill(2);
00140 
00141   printf("pick zVert=%f iV=%d\n", zVert,iVert);
00142 
00143   muMaker->muDst()->setVertexIndex(iVert);  
00144   Int_t nPrimTrAll=muMaker->muDst()->GetNPrimaryTrack();
00145   for(int itr=0;itr<nPrimTrAll;itr++) {
00146     StMuTrack *pr_track=muMaker->muDst()->primaryTracks(itr);
00147     assert(pr_track->vertexIndex()==iVert);
00148     if(pr_track->flag()<=0) continue;
00149     if(pr_track->flag() != 301) continue; // use TPC-only tracks
00150 
00151     hA[0]->Fill(10);
00152     hA[2]->Fill(pr_track->pt());
00153     if(pr_track->pt()< cut_primPt) continue;    
00154 
00155     hA[0]->Fill(11);
00156     hA[3]->Fill(pr_track->eta());
00157     if(fabs(pr_track->eta())> cut_primEta) continue;
00158 
00159     hA[0]->Fill(12);
00160     float dedx=pr_track->dEdx()*1e6;
00161     float mom= pr_track->p().mag();
00162     // printf("xxx %f %f\n", dedx,mom);
00163     hA[8]->Fill(mom,dedx);
00164     
00165     if(dedx<1.5) continue;
00166     if(dedx>cut_dedx) continue;
00167 
00168     hA[0]->Fill(13);
00169     
00170     const StMuTrack* globTr= pr_track->globalTrack();
00171     
00172     assert(globTr);
00173     assert(globTr->flag()>0);
00174     float nFF=(1.*globTr->nHitsFit())/ globTr->nHitsPoss();
00175     hA[4]->Fill(nFF);
00176     if(nFF< cut_nFitFrac) continue;
00177     hA[0]->Fill(14);
00178 
00179     float Rxy=globTr->lastPoint().perp();
00180     hA[10]->Fill(Rxy);
00181     if(Rxy<cut_primRxy) continue;
00182     hA[0]->Fill(15);
00183  
00184     StPhysicalHelixD TrkHlx=globTr->outerHelix();
00185     
00186     const int mxPL=2; // # of Radii to be tested       
00187     float RxyA[mxPL]={225.4, 248.4};// cm
00188     
00189     int softID=-888;
00190     StThreeVectorD pos3D;
00191     float posPhi=888;
00192     for(int ipl=0;ipl<mxPL;ipl++) { //..... top & botom of the Barrel tile 
00193       float Rbprs= RxyA[ipl];
00194       pairD  d2; 
00195       d2 = TrkHlx.pathLength(Rbprs);
00196       printf(" path ipl=%d =%f, 2=%f, period=%f, trR=%f\n",ipl, d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
00197       if(d2.first>=0 || d2.second<=0) {
00198         LOG_WARN<< Form("extrapolateTrk , unexpected solution for track crossing BPRS\n  d2.firts=%f, second=%f, track ignored",  d2.first, d2.second)<<endm;
00199         continue;
00200       }
00201       
00202       pos3D = TrkHlx.at(d2.second);
00203       double xmagn = ::sqrt( pos3D.x()*pos3D.x() + pos3D.y()*pos3D.y() );
00204       printf(" punchBPRS x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f  eta=%.2f\n",pos3D.x(),pos3D.y(),pos3D.z(),xmagn, globTr->eta());
00205       
00206       int softIDx;
00207       int ierr= mJanDbMaker->mBprsGeom->getId(pos3D.phi(),pos3D.pseudoRapidity(),softIDx);
00208       printf("hit tower id=%d, ierr=%d  ipl=%d\n",softIDx, ierr,ipl);
00209       
00210       posPhi=atan2(pos3D.y(),pos3D.x());
00211       if(posPhi<0) posPhi+=2*C_PI;      
00212       hA[5+ipl]->Fill(pos3D.z(),posPhi);       
00213       if(ierr) break; // do not try another depth for this track
00214       assert(softIDx>0 && softIDx<=mxBtow);
00215 
00216       if(ipl==0) { // BPRS layer
00217         hA[0]->Fill(16);
00218         if(!checkFiducial(pos3D.z(),posPhi,softIDx, Rbprs)) break; // to close to the edge     
00219         softID=softIDx;
00220         hA[0]->Fill(17);
00221       } else {
00222         hA[0]->Fill(18);
00223         if(softID!=softIDx) softID=-999;
00224         hA[0]->Fill(19);
00225         if(!checkFiducial(pos3D.z(),posPhi,softIDx, Rbprs)) softID=-999; // to close to the edge     
00226         // abort this track, it passes through 2 tiles
00227       }
00228     }
00229     
00230     if(softID<1) continue;
00231     nPrimTr++; 
00232     hA[0]->Fill(20);
00233     
00234     hA[14]->Fill(pr_track->pt());
00235     hA[15]->Fill(mom);    
00236     hA[7]->Fill(pos3D.z(),posPhi);       
00237     hA[12]->Fill(softID);
00238 
00239 #if 1 
00240     // ..................try any possible BPRS mapping, if activated
00241     if(hA[16])
00242       for(int j=0; j<mxBtow; j++) {
00243         int ibp=kBPrs;  // switch: BTOW or BPRS - only one!
00244         if(fullEve.statTile[ibp][j]) continue; // skip masked tiles
00245         float adc=fullEve.adcTile[ibp][j];
00246         if(adc<7 ) continue;
00247         if(adc>35 ) continue;
00248         hA[16]->Fill(j+1,softID);
00249       }
00250 #endif
00251 
00252     int id0=softID-1;
00253     uint isMip[mxBTile];
00254     uint badPed[mxBTile];
00255     float adc[mxBTile];
00256     for(int ibp=0;ibp<mxBTile;ibp++) {
00257       isMip[ibp]=0; // default =NO
00258       adc[ibp] =fullEve.adcTile [ibp][id0];
00259       badPed[ibp]=fullEve.statTile[ibp][id0];
00260       if( badPed[ibp]) continue;
00261       float adcL, adcH; // default, to use also tiles w/o gains
00262       mJanDbMaker->cut_mipAdcLH(ibp,softID,adcL,adcH);
00263       if(adc[ibp]<adcL) continue;
00264       if(adc[ibp]>adcH) continue;
00265       isMip[ibp]=1;
00266     }
00267  
00268     if(!badPed[kBPrs]) {  
00269       hA[0]->Fill(21);    
00270       hA[22]->Fill(softID,adc[kBPrs]);
00271       if(isMip[kBTow]) {
00272         hA[0]->Fill(22);    
00273         hA[23]->Fill(softID,adc[kBPrs]);
00274       }
00275     }
00276 
00277    if(!badPed[kBTow]) {  
00278       hA[0]->Fill(23);    
00279       hA[24]->Fill(softID,adc[kBTow]);
00280       if(isMip[kBPrs]) {
00281         hA[0]->Fill(24);    
00282         hA[25]->Fill(softID,adc[kBTow]);
00283       }
00284     }
00285 
00286    if(isMip[kBPrs] && isMip[kBTow])     hA[0]->Fill(25);
00287 
00288 
00289   } // end of loop over tracks
00290   
00291   if(nPrimTr<=0) return;
00292 
00293   hA[0]->Fill(3);
00294   hA[13]->Fill(zVert);
00295   if(nPrimTr>1) hA[0]->Fill(4);
00296 }
00297    
00298 //________________________________________________
00299 //________________________________________________
00300 int BarrelMipCalib::checkFiducial(float zTr, float phiTr, int softID, float Rxy){
00301   
00302   float etaHalfWidth=0.05/2.;
00303   const  float con_phiHalfWidth=(11.1/225.4)/2.; // radians
00304   float  cut_eps=cut_zMargin; // same for all except etaBin=20
00305   if((softID-1)%20==19)  cut_eps=cut_zMargin/2.; // be more forgiving
00306 
00307   float etaTw;
00308   assert( mJanDbMaker->mBprsGeom->getEta(softID,etaTw)==0);
00309  
00310   if(fabs(etaTw)<0.050) etaHalfWidth=(0.05-0.0035)/2.;
00311   if(fabs(etaTw)>0.95) etaHalfWidth=(0.9835-0.95)/2.;
00312  
00313   // those are edges of the tower in theta
00314   float thetaL=2.*atan(exp(-etaTw+etaHalfWidth));
00315   float thetaH=2.*atan(exp(-etaTw-etaHalfWidth));
00316  
00317   float zL=Rxy/tan(thetaL);
00318   float zH=Rxy/tan(thetaH);
00319                       
00320   float delZ1=zTr-zL;
00321   float delZ2=zH-zTr;
00322   bool isInZ= (delZ1>cut_eps) && (delZ2>cut_eps);
00323   printf("checkFiducial(zTr=%f phiTr=%f, id=%d  Rxy=%f\n",zTr,phiTr,softID,Rxy);
00324   printf(" theta range (deg)=%.1f %.1f\n",thetaL/3.1416*180,thetaH/3.1416*180); 
00325   printf("    etaTw=%f w2=%f  zL=%f zH=%f  delZ1=%f delZ2=%f isInZ=%d\n",etaTw,etaHalfWidth,zL,zH,delZ1, delZ2,isInZ);
00326 
00327  if(  isInZ) {
00328    hA[20]->Fill(delZ1);
00329    hA[20]->Fill(delZ2);
00330  }
00331 
00332   float phiTw;
00333   assert( mJanDbMaker->mBprsGeom->getPhi(softID,phiTw)==0);
00334   TVector2  uniTw; uniTw.SetMagPhi(1.,phiTw);
00335   TVector2  uniTr; uniTr.SetMagPhi(1.,phiTr);
00336   float delPhi=uniTw.DeltaPhi(uniTr);
00337 
00338   bool isInPhi=(con_phiHalfWidth- fabs(delPhi)) *Rxy> cut_eps;
00339                   
00340   printf("  phiTr=%f phiTw=%f delPhi=%f  delRphi=%f isInPhi=%d\n", phiTr, phiTw,delPhi, delPhi*Rxy, isInPhi);
00341  if(isInPhi) hA[21]->Fill(delPhi *Rxy);
00342 
00343   return isInZ && isInPhi ;
00344 
00345 }
00346   
00347 
00348 // test code to find where in eta are towers at eat=1.0
00349 //________________________________________________
00350 //________________________________________________
00351 void BarrelMipCalib::searchEtaBin20(  JanBarrelEvent &fullEve){
00352 
00353   // Access to muDst .......................
00354   StMuEvent* muEve = muMaker->muDst()->event();
00355   int nPrimV=muMaker->muDst()->numberOfPrimaryVertices();
00356   hA[0]->Fill(0);
00357   
00358   printf("\nWARN-ALTERNATIVE CODE: mipCalib eventID %d nPrimV=%d  =============\n", muEve->eventId(),nPrimV);
00359 
00360   // ..............filter events , keep only those with MIP tower at eta=1.0
00361 
00362   int pickBin=17;
00363   int myId0=-1;
00364   for(int j=0; j<mxBtow; j++) {
00365     if((j%20)!=pickBin-1) continue; // fixed eta-bin
00366     if((j/20)%3) continue; // every 3rd phi- row
00367 
00368     int ibp=kBTow;  // switch: BTOW or BPRS - only one!
00369     if(fullEve.statTile[ibp][j]) continue; // skip masked tiles
00370     float adc=fullEve.adcTile[ibp][j];
00371     if(adc<10 ) continue;
00372     if(adc>30 ) continue;
00373     myId0=j;
00374     break;
00375   }
00376   if(myId0<0) return;
00377   hA[0]->Fill(1);
00378   float phiTw;
00379   assert( mJanDbMaker->mBprsGeom->getPhi(myId0+1,phiTw)==0);
00380   TVector2  uniTw; uniTw.SetMagPhi(1.,phiTw);
00381 
00382   if(nPrimV<=0) return;
00383   hA[0]->Fill(2);
00384 
00385   
00386   int nPrimTr =0;
00387   
00388   // find vertex with smalest Z - to compensate for PPV error in FMS-production
00389 
00390   int iVert=-1;
00391   float zVert=99999;
00392   
00393   for(int iv=0;iv<nPrimV;iv++) {
00394     StMuPrimaryVertex* V= muMaker->muDst()->primaryVertex(iv);
00395     assert(V);
00396     if(V->ranking()<0.) continue; //drop vertices with too few matched tracks
00397     muMaker->muDst()->setVertexIndex(iv);
00398     const StThreeVectorF &r=V->position();
00399     const StThreeVectorF &er=V->posError();
00400     if(iv==0) hA[9]->Fill(r.z()); // histo only highe rank vertex
00401 
00402     printf("iv=%d   Vz=%.2f +/-%.2f  bestZ=%f\n",iv,r.z(),er.z() ,zVert );
00403   
00404     if( fabs(zVert)< fabs(r.z()))  continue;
00405     zVert=r.z(); 
00406     iVert=iv;
00407   }
00408 
00409   //  printf("iVert=%d\n",iVert);
00410   if(iVert<0) return; // non of vertices has positive rank
00411 
00412   hA[1]->Fill(zVert);
00413   if(fabs(zVert) > cut_zVertex) return; // drop if too large Z 
00414   hA[0]->Fill(3);
00415 
00416   printf("pick zVert=%f iV=%d\n", zVert,iVert);
00417 
00418   muMaker->muDst()->setVertexIndex(iVert);  
00419   Int_t nPrimTrAll=muMaker->muDst()->GetNPrimaryTrack();
00420   for(int itr=0;itr<nPrimTrAll;itr++) {
00421     StMuTrack *pr_track=muMaker->muDst()->primaryTracks(itr);
00422     assert(pr_track->vertexIndex()==iVert);
00423     if(pr_track->flag()<=0) continue;
00424 
00425     hA[0]->Fill(10);
00426     hA[2]->Fill(pr_track->pt());
00427     if(pr_track->pt()< cut_primPt) continue;    
00428 
00429     hA[0]->Fill(11);
00430     hA[3]->Fill(pr_track->eta());
00431     if(fabs(pr_track->eta())> cut_primEta) continue;
00432 
00433     hA[0]->Fill(12);
00434     float dedx=pr_track->dEdx()*1e6;
00435     float mom= pr_track->p().mag();
00436     // printf("xxx %f %f\n", dedx,mom);
00437     hA[8]->Fill(mom,dedx);
00438     
00439     if(dedx<1.5) continue;
00440     if(dedx>cut_dedx) continue;
00441 
00442     hA[0]->Fill(13);
00443     
00444     const StMuTrack* globTr= pr_track->globalTrack();
00445     
00446     assert(globTr);
00447     assert(globTr->flag()>0);
00448     float nFF=(1.*globTr->nHitsFit())/ globTr->nHitsPoss();
00449     hA[4]->Fill(nFF);
00450     if(nFF< cut_nFitFrac) continue;
00451     
00452     hA[0]->Fill(14);
00453     StPhysicalHelixD TrkHlx=globTr->outerHelix();
00454     
00455     const int mxPL=2; // # of Radii to be tested       
00456     float RxyA[mxPL]={225.4, 248.4};// cm
00457     
00458     int softID=-888;
00459     StThreeVectorD pos3D;
00460     float posPhi=888;
00461     for(int ipl=0;ipl<1;ipl++) { //..... top & botom of the Barrel tile 
00462       float Rbprs= RxyA[ipl];
00463       pairD  d2; 
00464       d2 = TrkHlx.pathLength(Rbprs);
00465       printf(" path ipl=%d =%f, 2=%f, period=%f, trR=%f\n",ipl, d2.first ,d2.second,TrkHlx.period(),1./TrkHlx.curvature());
00466       if(d2.first>=0 || d2.second<=0) {
00467         LOG_WARN<< Form("extrapolateTrk , unexpected solution for track crossing BPRS\n  d2.firts=%f, second=%f, track ignored",  d2.first, d2.second)<<endm;
00468         continue;
00469       }
00470       
00471       pos3D = TrkHlx.at(d2.second);
00472       double xmagn = ::sqrt( pos3D.x()*pos3D.x() + pos3D.y()*pos3D.y() );
00473       printf(" punchBPRS x,y,z=%.1f, %.1f, %.1f, Rxy=%.1f  eta=%.2f\n",pos3D.x(),pos3D.y(),pos3D.z(),xmagn, globTr->eta());
00474       
00475       posPhi=atan2(pos3D.y(),pos3D.x());
00476       if(posPhi<0) posPhi+=2*C_PI;      
00477       hA[5+ipl]->Fill(pos3D.z(),posPhi);    
00478    
00479       TVector2  uniTr; uniTr.SetMagPhi(1.,posPhi);
00480       float delPhi=uniTw.DeltaPhi(uniTr);
00481 
00482       if(fabs(delPhi)>0.15) continue;
00483       softID=111; // take every track
00484     }
00485     
00486     if(softID<1) continue;
00487     nPrimTr++; 
00488     hA[0]->Fill(19);
00489     
00490     hA[14]->Fill(pr_track->pt());
00491     hA[15]->Fill(mom);    
00492     hA[7]->Fill(pos3D.z(),posPhi);       
00493     hA[17]->Fill(pos3D.z());
00494     if(pos3D.z()>0) hA[18]->Fill(posPhi);       
00495     //  hA[12]->Fill(softID);
00496     
00497   } // end of loop over tracks
00498   
00499   if(nPrimTr<=0) return;
00500   
00501   hA[0]->Fill(4);
00502   hA[13]->Fill(zVert);
00503   if(nPrimTr>1) hA[0]->Fill(5);
00504 }
00505    
00506   
00507