// $Id$
+/**
+ trigger_pp -> Calculate efficiency of trigger and fake tracks in pp pileup events
+ plot_pp -> Plot vertex resolution
+ plot_pp_vert -> Plot vertex resolution
+ If you are looking for the old trigger macro search for
+ OLD_TRIGGER_MACRO.
+
+ Author: C. Loizides
+*/
+
+
+/* ---trigger_pp---
+ Calculate the efficiency of found tracks in pileup event
+ and found fake tracks after applying different DCA cuts to
+ the assumed vertex of (0,0,0). Tracks are stored in the ntuple
+ root file (see fill_pp.C). Good tracks can either be used from
+ the single pp event or the good_tracks_tpc event (the macro
+ will detect which ntuple file you specified).
+
+ Script: Track single pp events and store in ntuple:
+--------------------------------------------------------
+#!/bin/bash
+
+dir=/mnt/data/loizides/alidata/head
+wdir=$dir/pp/getrackted-0-375/
+mdir=~/work/l3/level3code/exa/
+sfile=~/pp-tracks.root
+rfile=$dir/PythiaPP_1000Events_TPConly_ev0-375_digits.root
+
+cd $dir
+ln -s -f $rfile digitfile.root
+cd $mdir;
+
+for i in `seq 0 365`; do
+ echo $i
+ trigev=$i
+ cd $mdir
+ aliroot -b -q runtracker_pp.C\($trigev,0,\"$rfile\",\"$wdir\"\)
+ aliroot -b -q fill_pp.C\($trigev,\"$wdir\",\"$sfile\"\)
+done
+--------------------------------------------------------
+
+ Script: Track pileup events and store in ntuple:
+-----------------------------------------------------
+#!/bin/bash
+
+dir=/mnt/data/loizides/alidata/head
+wdir=/mnt/data/loizides/alidata/head/pp/pileup-25
+mdir=~/l3/level3code/exa
+sfile=~/pileup25-tracks.root
+
+cd $mdir;
+
+for i in `seq 0 100`; do
+ echo $i
+ pdir=$wdir/pileup-$i/
+ cd $pdir
+ trigev=`ls digits_*_0_-1.raw | cut -d_ -f 2`
+ cd $mdir
+ aliroot -b -q runtracker_pp.C\($trigev,\"$pdir\",0,\"$pdir\"\)
+ aliroot -b -q fill_pp.C\($trigev,\"$pdir\",\"$sfile\",$i\)
+done
+-----------------------------------------------------
+
+ Parameter to trigger_pp are:
+ Pileupfile is the ntuple root file containing the information of the pileup event
+ Ppfile gives the ntuple root file containing the reference information of the single pp event (can either be good_particles or single tracked pp events)
+ Evi and Eve give start and end event number to process.
+*/
+
+void trigger_pp(Char_t *pileupfile, Char_t *ppfile, Int_t evi=0, Int_t eve=100)
+{
+
+ //-----------------------------------------------------------
+ const Char_t *LTEXT_="for 25 piles, min. 2 trigger tracks";
+ const Int_t MINTRACKS_=2; //min. tracks of the triggered
+ const Int_t N_=100; //number of data points
+ const Float_t WIDTH_=1; //the width of the xbin
+ //-----------------------------------------------------------
+
+
+ //"pt:phi:eta:xvert:yvert:zvert:imp:nhits:px:py:pz:event:mc")
+ Float_t pt1,phi1,eta1,xvert1,yvert1,zvert1,imp1,nhits1,px1,py1,pz1,event1,mc1; //pileup
+ Float_t pt2,phi2,eta2,xvert2,yvert2,zvert2,imp2,nhits2,px2,py2,pz2,event2,mc2; //pp
+
+ Int_t countgoodevents=0;
+ Int_t nent1,nent2;
+ Char_t dummy[1024];
+
+ Int_t tnorm;
+ Int_t fnorm;
+ Int_t tcounts[N_];
+ Int_t fcounts[N_];
+
+ Float_t tallcounts[N_];
+ Float_t fallcounts[N_];
+ Int_t totnorm=0;
+
+ Float_t thcounts[N_];
+ Float_t fhcounts[N_];
+ Float_t thmeancounts[N_];
+ Float_t fhmeancounts[N_];
+ Float_t thsigmacounts[N_];
+ Float_t fhsigmacounts[N_];
+ Float_t xcut[N_];
+ for(Int_t k=0;k<N_;k++){
+ tallcounts[k]=0;
+ fallcounts[k]=0;
+ thmeancounts[k]=0;
+ fhmeancounts[k]=0;
+ thsigmacounts[k]=0;
+ fhsigmacounts[k]=0;
+ xcut[k]=(k+1)*WIDTH_;
+ }
+
+ TFile file1=TFile(pileupfile,"READ");
+ TFile file2=TFile(ppfile,"READ");
+ if(!file1.IsOpen() || !file2.IsOpen()) return;
+
+ for(Int_t i=evi;i<eve;i++){ //loop over pileup-events
+
+ sprintf(dummy,"good-pptracks-%d",i);
+ TNtuple *ntuppel1 = (TNtuple*)file1.Get(dummy);
+ Float_t nhits1=0,pdg1=0,sector1=0;
+ if(ntuppel1){
+ ntuppel1->SetBranchAddress("mc",&mc1);
+ ntuppel1->SetBranchAddress("pdg",&pdg1);
+ ntuppel1->SetBranchAddress("px",&px1);
+ ntuppel1->SetBranchAddress("py",&py1);
+ ntuppel1->SetBranchAddress("pz",&pz1);
+ ntuppel1->SetBranchAddress("xvert",&xvert1);
+ ntuppel1->SetBranchAddress("yvert",&yvert1);
+ ntuppel1->SetBranchAddress("zvert",&zvert1);
+ ntuppel1->SetBranchAddress("nhits",&nhits1);
+ ntuppel1->SetBranchAddress("sector",§or1);
+ event1=i;
+ } else {
+ sprintf(dummy,"pptracks-%d",i);
+ TNtuple *ntuppel1 = (TNtuple*)file1.Get(dummy);
+ if(!ntuppel1) continue;
+ ntuppel1->SetBranchAddress("pt",&pt1);
+ ntuppel1->SetBranchAddress("phi",&phi1);
+ ntuppel1->SetBranchAddress("eta",&eta1);
+ ntuppel1->SetBranchAddress("xvert",&xvert1);
+ ntuppel1->SetBranchAddress("yvert",&yvert1);
+ ntuppel1->SetBranchAddress("zvert",&zvert1);
+ ntuppel1->SetBranchAddress("imp",&imp1);
+ ntuppel1->SetBranchAddress("px",&px1);
+ ntuppel1->SetBranchAddress("py",&py1);
+ ntuppel1->SetBranchAddress("pz",&pz1);
+ ntuppel1->SetBranchAddress("event",&event1);
+ ntuppel1->SetBranchAddress("mc",&mc1);
+ }
+
+ tnorm=0;
+ fnorm=0;
+ for(Int_t j=0;j<N_;j++){
+ tcounts[j]=0;
+ fcounts[j]=0;
+ thcounts[j]=0;
+ fhcounts[j]=0;
+ }
+
+ nent1=ntuppel1->GetEntries();
+ for(Int_t j=0;j<nent1;j++){
+ ntuppel1->GetEvent(j);
+ if(nhits1!=0){ // file1 is good-particle file
+ pt1=TMath::Sqrt(px1*px1+py1*py1);
+ Float_t p=TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
+ eta1=0.5*TMath::log((p+pz1)/(p-pz1));
+ phi1=TMath::atan(py1/px1);
+ imp1=TMath::Sqrt(xvert1*xvert1+yvert1*yvert1+zvert1*zvert1);
+ }
+ if(j==0){ //get info of triggered event from file2
+ sprintf(dummy,"good-pptracks-%d",event1);
+ TNtuple *ntuppel2 = (TNtuple*)file2.Get(dummy);
+ if(ntuppel2){
+ Float_t nhits2,pdg2,sector2;
+ ntuppel2->SetBranchAddress("mc",&mc2);
+ ntuppel2->SetBranchAddress("pdg",&pdg2);
+ ntuppel2->SetBranchAddress("px",&px2);
+ ntuppel2->SetBranchAddress("py",&py2);
+ ntuppel2->SetBranchAddress("pz",&pz2);
+ ntuppel2->SetBranchAddress("xvert",&xvert2);
+ ntuppel2->SetBranchAddress("yvert",&yvert2);
+ ntuppel2->SetBranchAddress("zvert",&zvert2);
+ ntuppel2->SetBranchAddress("nhits",&nhits2);
+ ntuppel2->SetBranchAddress("sector",§or2);
+ event2=event1;
+ nent2=ntuppel2->GetEntries();
+ for(Int_t k=0;k<nent2;k++){
+ ntuppel2->GetEvent(k);
+ if(mc2<0)continue; //dont count fake tracks
+ //if(TMath::fabs(zvert2)>30) continue;
+ pt2=TMath::Sqrt(px2*px2+py2*py2);
+ if(nhits2<63) continue;
+ if(pt2<0.1) continue;
+ if((px2==0)&&(py2==0)&&(pz2==0)) continue;
+ Float_t p=TMath::Sqrt(px2*px2+py2*py2+pz2*pz2);
+ eta2=0.5*TMath::log((p+pz2)/(p-pz2));
+ if(TMath::fabs(eta2)>0.9) continue;
+ phi2=TMath::atan(py2/px2);
+ imp2=TMath::Sqrt(xvert2*xvert2+yvert2*yvert2+zvert2*zvert2);
+ tnorm++;
+ }
+ } else {
+ sprintf(dummy,"pptracks-%d",event1);
+ TNtuple *ntuppel2 = (TNtuple*)file2.Get(dummy);
+ if(ntuppel2){
+ ntuppel2->SetBranchAddress("pt",&pt2);
+ ntuppel2->SetBranchAddress("phi",&phi2);
+ ntuppel2->SetBranchAddress("eta",&eta2);
+ ntuppel2->SetBranchAddress("xvert",&xvert2);
+ ntuppel2->SetBranchAddress("yvert",&yvert2);
+ ntuppel2->SetBranchAddress("zvert",&zvert2);
+ ntuppel2->SetBranchAddress("imp",&imp2);
+ ntuppel2->SetBranchAddress("px",&px2);
+ ntuppel2->SetBranchAddress("py",&py2);
+ ntuppel2->SetBranchAddress("pz",&pz2);
+ ntuppel2->SetBranchAddress("event",&event2);
+ ntuppel2->SetBranchAddress("mc",&mc2);
+
+ nent2=ntuppel2->GetEntries();
+ for(Int_t k=0;k<nent2;k++){
+ ntuppel2->GetEvent(k);
+ if(mc2<0)continue; //dont count fake tracks
+ if(pt2<0.1) continue;
+ //if(TMath::fabs(zvert2)>30) continue;
+ if(TMath::fabs(eta2)>0.9) continue;
+ tnorm++;
+ //cout << k << ": " << pt2 << " " << mc2 <<endl;
+ }
+ }
+ }
+ }
+ if(tnorm<MINTRACKS_){
+ tnorm=0;
+ break; //triggered event has to have one track at least
+ }
+
+ if(pt1<0.1) continue;
+ if(TMath::fabs(eta1>0.9)) continue;
+
+ if(mc1<-1) continue; //dont count fake tracks
+ else if(mc1==-1) fnorm++;
+
+ for(Int_t k=0;k<N_;k++){ //do the counting
+ Float_t cut=xcut[k];
+ //Float_t impt=TMath::Sqrt(xvert1*xvert1+yvert1*yvert1);
+ Float_t impz=TMath::abs(zvert1);
+ if(impz>cut) continue;
+ if(mc1==-1) fcounts[k]++;
+ else tcounts[k]++;
+ }
+
+ }//pp track loop
+
+ if(tnorm==0) continue; //break from above
+ countgoodevents++; //for normalization
+
+ totnorm+=tnorm;
+
+ for(Int_t k=0;k<N_;k++){ //do the counting
+ if(tcounts[k]> tnorm)tcounts[k]=tnorm;
+ thcounts[k]=(Float_t)tcounts[k]/tnorm;
+ fhcounts[k]=(Float_t)fcounts[k]/tnorm;
+ tallcounts[k]+=tcounts[k];
+ fallcounts[k]+=fcounts[k];
+ thmeancounts[k]+=thcounts[k];
+ fhmeancounts[k]+=fhcounts[k];
+ thsigmacounts[k]+=thcounts[k]*thcounts[k];
+ fhsigmacounts[k]+=fhcounts[k]*fhcounts[k];
+ if((k==N_-1)&&(thcounts[k]< 0.99)){
+ //cout << "Warning " << i << " " << event2 << " " << tcounts[k] << " " << tnorm << endl;
+ }
+ }
+ }//pileup loop
+
+ file1.Close();
+ file2.Close();
+
+ cout << "Events used: " << countgoodevents << endl;
+ for(Int_t k=0;k<N_;k++){ //do the counting
+ thmeancounts[k]/=countgoodevents;
+ fhmeancounts[k]/=countgoodevents;
+
+ tallcounts[k]/=totnorm;
+ fallcounts[k]/=totnorm;
+
+ if(countgoodevents>1){
+ Float_t stemp=thsigmacounts[k]/countgoodevents;
+ Float_t s2temp=fhsigmacounts[k]/countgoodevents;
+ Int_t N=countgoodevents-1;
+ thsigmacounts[k]=TMath::sqrt((stemp -thmeancounts[k]*thmeancounts[k])/N);
+ fhsigmacounts[k]=TMath::sqrt((s2temp-fhmeancounts[k]*fhmeancounts[k])/N);
+ }else{
+ thsigmacounts[k]=0;
+ fhsigmacounts[k]=0;
+ }
+ cout << k << ": " << thmeancounts[k] << "+-"<< thsigmacounts[k] << " " << fhmeancounts[k] << "+-" << fhsigmacounts[k] << endl;
+ }
+
+ TCanvas *c1 = new TCanvas("c1","PileUp",1000,500);
+ //c1->SetFillColor(42);
+ //c1->SetGrid();
+ //c1->GetFrame()->SetFillColor(21);
+ //c1->GetFrame()->SetBorderSize(12);
+
+ TGraphErrors *g1=new TGraphErrors(N_,xcut,thmeancounts,0,&thsigmacounts[0]);
+ //TGraphErrors *g1=new TGraphErrors(N_,xcut,tallcounts);
+ sprintf(dummy,"Good tracks %s",LTEXT_);
+ g1->SetTitle(dummy);
+ g1->SetMarkerColor(4);
+ g1->GetHistogram()->SetXTitle("Z_{cut} [cm]");
+ //g1->GetHistogram()->SetYTitle("Efficiency");
+ g1->SetMarkerStyle(21);
+
+ TGraphErrors *g2=new TGraphErrors(N_,xcut,fhmeancounts,0,fhsigmacounts);
+ //TGraphErrors *g2=new TGraphErrors(N_,xcut,fallcounts);
+ g2->SetTitle("Fake tracks (relative to good triggered tracks)");
+ g2->GetHistogram()->SetXTitle("Z_{cut} [cm]");
+ g2->SetMarkerColor(4);
+ g2->SetMarkerStyle(21);
+
+ c1->Divide(2,1);
+ c1->cd(1);
+ g1->Draw("AP");
+ c1->cd(2);
+ g2->Draw("AP");
+ c1->Update();
+}
+
+/* Plot the vertex reconstruction resolution */
+
+void plot_pp_vert(Int_t evi=0,Int_t evs=100,Char_t *infile="pp-tracks.root")
+{
+ const Int_t zcut=3;
+ const Int_t zcutplot=5;
+
+ Char_t dummy[1000];
+ Char_t fname[1000];
+ Char_t fname2[1000];
+
+ Float_t pt,phi,eta,xvert,yvert,zvert,imp,nhits,px,py,pz,event,mc;
+ TH1F *vhist = new TH1F("zhist","Vertex reconstruction resolution",50,-zcutplot,zcutplot);
+
+ TFile file = TFile(infile,"READ");
+ if(!file.IsOpen()) return;
+
+ TNtuple *ntuppel;
+
+ for(int ev=evi; ev<evs; ev++){
+ //sprintf(dummy,"good-pptracks-%d",ev); /*dont use this*/
+ sprintf(dummy,"pptracks-%d",ev);
+
+ if(ntuppel) delete ntuppel;
+ ntuppel = (TNtuple*)file.Get(dummy);
+ if(!ntuppel) continue;
+ if(ntuppel.GetEntries()==0) continue;
+
+ sprintf(fname,"(eta >-0.9) && (eta<0.9) && (pt>0.1) && (nhits>63) && (zvert<%d) && (zvert>-%d)",zcut,zcut);
+ ntuppel->Draw("zvert>>histo",fname,"groff");
+ Float_t mean = histo->GetMean();
+ //cout << mean << endl;
+ vhist->Fill(mean);
+ }
+
+ gStyle->SetStatColor(10);
+ gStyle->SetOptFit(1);
+
+ TF1 *f1 = new TF1("f1","gaus",-(Float_t)zcutplot/2.,(Float_t)zcutplot/2.);
+ vhist->Fit("f1","R");
+
+ vhist->SetXTitle("Z* [cm]");
+ vhist->Draw("E1P");
+}
+
+/* Plot the impact parameter resolution */
+
+void plot_pp(Int_t evi=0,Int_t evs=100,Char_t *infile="pp-tracks.root")
+{
+ const Int_t zcut=3;
+ const Int_t zcutplot=5;
+
+ Char_t dummy[1000];
+ Char_t fname[1000];
+ Char_t fname2[1000];
+
+ TFile file = TFile(infile,"READ");
+ if(!file.IsOpen()) return;
+
+ TH1F *vhist = new TH1F("zhist","Impact parameter resolution",50,-zcut,zcut);
+
+ TNtuple *ntuppel;
+ Int_t a=0;
+ for(int ev=evi; ev<evs; ev++){
+ //sprintf(dummy,"good-pptracks-%d",ev); /*dont use this*/
+ sprintf(dummy,"pptracks-%d",ev);
+
+ if(ntuppel) delete ntuppel;
+ ntuppel = (TNtuple*)file.Get(dummy);
+ if(!ntuppel) continue;
+ if(ntuppel.GetEntries()==0) continue;
+
+ sprintf(fname,"(eta>-0.9) && (eta<0.9) && (pt>0.1) && (nhits>63) && (zvert<%d) && (zvert>-%d)",zcut,zcut);
+ ntuppel->Draw("zvert>>histo",fname,"groff");
+ Float_t mean = histo->GetMean();
+ a++;
+ //cout << mean << endl;
+
+ sprintf(fname2,"(zvert-%f)>>+zhist",mean);
+ ntuppel->Draw(fname2,fname);
+ }
+
+ gStyle->SetStatColor(10);
+ gStyle->SetOptFit(1);
+ cout << a <<endl;
+ TF1 *f1 = new TF1("f1","gaus",-(Float_t)zcutplot/2.,(Float_t)zcutplot/2.);
+ vhist->Fit("f1","R");
+
+ vhist->SetXTitle("Z_{DCA}-Z* [cm]");
+ vhist->DrawCopy("E1P");
+ file->Close();
+}
+
+
+#ifdef OLD_TRIGGER_MACRO
void trigger_pp(char *outfile="results.root")
{
for(int event=0; event<1; event++)
{
char fname[256];
- sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/tracks.raw",event);
- //sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/tracks.raw",event);
+ sprintf(fname,"/data1/AliRoot/pp/pileup/tracks_0.raw");
+ //sprintf(fname,"aliruntfile.root");
//Get the tracks:
- AliL3TrackArray *tracks = new AliL3TrackArray();
+ /*AliL3TrackArray *tracks = new AliL3TrackArray();
AliL3FileHandler *file = new AliL3FileHandler();
file->SetBinaryInput(fname);
file->Binary2TrackArray(tracks);
file->CloseBinaryInput();
- delete file;
+ delete file;*/
+
+ sprintf(fname,"/data1/AliRoot/pp/pileup/");
+
+ AliL3Evaluate *eval=new AliL3Evaluate(fname,63,63);
+ eval->LoadData(event,-1);
+ eval->AssignIDs();
- sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/pileups/recon_%d/",event);
+ AliL3TrackArray *tracks=eval->GetTracks();
+
+ sprintf(fname,"/data1/AliRoot/pp/pileup/");
//sprintf(fname,"/prog/alice/data/Rawdata/1_patch/pp/recon_%d/",event);
Int_t ntracks=0;
Double_t xc,yc,zc;
Double_t impact;
AliL3Vertex vertex;
- AliL3Fitter *fitter = new AliL3Fitter(&vertex);
- fitter->LoadClusters(fname);
- //fitter->NoVertex();
-
+
+ Int_t mcid = 0;
+
AliL3TrackArray *ftracks = new AliL3TrackArray();
for(int i=0; i<tracks->GetNTracks(); i++)
{
track = (AliL3Track*)tracks->GetCheckedTrack(i);
if(!track) continue;
+
+ //Assign MCid
+ mcid = track->GetMCid();
+ //if (mcid != -1)
+ //cout << "MCid " << mcid << endl;
+
track->CalculateHelix();
//cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
//cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
- fitter->FitHelix(track);//refit the tracks
+
track->CalculateHelix();
track->GetClosestPoint(&vertex,xc,yc,zc);
meas[0]=track->GetPt();
if(fabs(track->GetPseudoRapidity())>0.9) continue;
if(track->GetNHits() < 100) continue;
if(track->GetPt()<0.2) continue;
+
impact = sqrt(xc*xc+yc*yc+zc*zc);
if(fabs(zc)>3) continue;
ftracks->AddLast(track);
}
cout<<endl<<"There was "<<ntracks<<" accepted tracks, out of total "<<tracks->GetNTracks()<<" found tracks"<<endl;
- display(ftracks,fname);
+ //display(ftracks,fname);
delete tracks;
- delete fitter;
delete ftracks;
}
ntuppel->Write();
rfile->Close();
delete ntuppel;
-
}
void display(AliL3TrackArray *tracks,char *path)
{
int slice[2]={0,35};
d = new AliL3Display(slice);
- d->Setup("tracks.raw",path);
+ d->Setup("tracks_0.raw",path);
d->SetTracks(tracks);
//d->DisplayClusters();
- d->DisplayAll();
- //d->DisplayTracks();
+ d->DisplayAll();
+ //d->DisplayTracks();
}
void ploteff()
{
- gROOT->LoadMacro("XFunct.C");
//double x[6]={1,2,4,6,8,10};
//double y[6]={0.873684,1.01379,1.17751,1.28614,1.31638,1.32022};
double geteff(char *infile,char *singlefile,double cut)
{
- gROOT->LoadMacro("XFunct.C");
gStyle->SetStatColor(10);
gStyle->SetOptFit(1);
Int_t pileups[25],single[25];
void plot(char *infile)
{
- gROOT->LoadMacro("XFunct.C");
gStyle->SetStatColor(10);
gStyle->SetOptFit(1);
file = TFile::Open(infile,"READ");
gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
loadlibs();
- TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
+ //TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
+ TFile *rootfile = TFile::Open("/prog/alice/pp/pileup/100pileup/alirunfile.root","READ");
gAlice = (AliRun*)rootfile->Get("gAlice");
// TNtuple *ntup = new TNtuple(
// cuts:
- Double_t vertcut = 0.001; // cut on the vertex position
+ //Double_t vertcut = 0.001; // cut on the vertex position
+ Double_t vertcut = 2.0;
Double_t decacut = 3.; // cut if the part. decays close to the vert.
Double_t timecut = 0.;
if(xxx)continue;
TParticlePDG *ppdg = part->GetPDG();
- if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
+ //if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
+ if(TMath::Abs(ppdg->Charge())<1)continue; // only charged (no quarks)
Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
if(dist>vertcut)continue; // cut on the vertex
return nprch1-nprch2;
}
+#endif