#include "AliITSPid.h" #include "TMath.h" #include "AliITSIOTrack.h" #include ClassImp(AliITSPid) Float_t AliITSPid::qcorr(Float_t xc) { assert(0); Float_t fcorr; fcorr=( 0.766 +0.9692*xc -1.267*xc*xc )*( 1.-TMath::Exp(-xc*64.75) ); if(fcorr<=0.1)fcorr=0.1; return qtot/fcorr; } Float_t AliITSPid::qtrm(Int_t track) { TVector q(*( this->GetVec(track) )); Int_t ml=(Int_t)q(0); if(ml<1)return 0.; if(ml>6)ml=6; float vf[6]; Int_t nl=0; for(Int_t i=0;ifSigmin){vf[nl]=q(i+1);nl++;}} if(nl==0)return 0.; switch(nl){ case 1:q(6)=q(1); break; case 2:q(6)=(q(1)+q(2))/2.; break; default: for(int fi=0;fi<2;fi++){ Int_t swap; do{ swap=0; float qmin=vf[fi]; for(int j=fi+1;jvf[j]){qmin=vf[j];vf[j]=vf[fi];vf[fi]=qmin;swap=1;}; }while(swap==1); } q(6)= (vf[0]+vf[1])/2.; break; } for(Int_t i=0;iSetVec(track,q); return (q(6)); } Float_t AliITSPid::qtrm(Float_t qarr[6],Int_t narr) { Float_t q[6],qm,qmin; Int_t nl,ml; if(narr>0&&narr<7){ml=narr;}else{return 0;}; nl=0; for(Int_t i=0;ifSigmin){q[nl]=qarr[i];nl++;}} if(nl==0)return 0.; switch(nl){ case 1:qm=q[0]; break; case 2:qm=(q[0]+q[1])/2.; break; default: Int_t swap; for(int fi=0;fi<2;fi++){ do{ swap=0; qmin=q[fi]; for(int j=fi+1;jq[j]){qmin=q[j];q[j]=q[fi];q[fi]=qmin;swap=1;}; }while(swap==1); } qm= (q[0]+q[1])/2.; break; } return qm; } //----------------------------------------------------------- //inline Int_t AliITSPid::wpik(Int_t nc,Float_t q) { // return pion(); Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk; Float_t appi,apk; qmpi =cut[nc][1]; sigpi=cut[nc][2]; qmk =cut[nc][3]; sigk =cut[nc][4]; appi = aprob[0][nc-5]; apk = aprob[1][nc-5]; // cout<<"qmpi,sigpi,qmk,sigk="<ppi){return kaon();}else{return pion();} } //----------------------------------------------------------- //inline Int_t AliITSPid::wpikp(Int_t nc,Float_t q) { Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp; Float_t appi,apk,app; qmpi =cut[nc][1]; sigpi=cut[nc][2]; qmk =cut[nc][3]; sigk =cut[nc][4]; qmp =cut[nc][5]; sigp =cut[nc][6]; //appi = apk = app = 1.; appi = aprob[0][nc-5]; apk = aprob[1][nc-5]; app = aprob[2][nc-5]; //ppi =TMath::Erf( TMath::Abs((q-qmpi)/sigpi) )+0.5; //pka =TMath::Erf( TMath::Abs((q-qmk )/sigk) )+0.5; //ppr =TMath::Erf( TMath::Abs((q-qmp )/sigp) )+0.5; ppi=appi*TMath::Gaus(q,qmpi,sigpi) /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) ); pk = apk*TMath::Gaus(q,qmk, sigk ) /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) ); pp = app*TMath::Gaus(q,qmp, sigp ) /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) ); fWp=pp; fWpi=ppi; fWk=pk; //cout<<" wpikp: mid,sig pi,k,p="<pk&&ppi>pp ) { return pion(); } if(pk>pp){return kaon();}else{return proton();} } //----------------------------------------------------------- Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm) { return 0; } //----------------------------------------------------------- Int_t AliITSPid::GetPcode(AliTPCtrack *track) { Double_t xk,par[5]; track->GetExternalParameters(xk,par); Float_t phi=TMath::ASin(par[2]) + track->GetAlpha(); if (phi<-TMath::Pi()) phi+=2*TMath::Pi(); if (phi>=TMath::Pi()) phi-=2*TMath::Pi(); Float_t lam=TMath::ATan(par[3]); Float_t pt_1=TMath::Abs(par[4]); Float_t mom=1./(pt_1*TMath::Cos(lam)); Float_t dedx=track->GetdEdx(); Int_t pcode=GetPcode(dedx/40.,mom); // cout<<"TPCtrack dedx,mom,pcode="<GetPx(); py=track->GetPy(); pz=track->GetPz(); Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz); //??????????????????? // Float_t dedx=1.0; Float_t dedx=track->GetdEdx(); //??????????????????? Int_t pcode=GetPcode(dedx,mom); // cout<<"ITSV1 dedx,mom,pcode="<Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848); track->PropagateTo(3.,0.0028,65.19); track->PropagateToVertex(); Double_t xk,par[5]; track->GetExternalParameters(xk,par); Float_t lam=TMath::ATan(par[3]); Float_t pt_1=TMath::Abs(par[4]); Float_t mom=0.; if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;}; Float_t dedx=track->GetdEdx(); // cout<<"lam,pt_1,mom,dedx="<cut[5][5] ) {return proton();} else {return wpik(5,q);}; //6)----------------------470-530 Mev/c ------------- if ( pm<=cut[6][0] ) if ( q>cut[6][5] ) {return proton();} else {return wpik(6,q);}; //7)----------------------530-590 Mev/c ------------- if ( pm<=cut[7][0] ) if ( q<=cut[7][5] ) {return wpik(7,q);} else {return proton();}; //8)----------------------590-650 Mev/c ------------- if ( pm<=cut[8][0] ) if ( q<=cut[8][5] ) {return wpik(8,q);} else {return proton();}; //9)----------------------650-730 Mev/c ------------- if ( pm<=cut[9][0] ) if ( q<=cut[9][5] ) {return wpik(9,q);} else {return proton();}; //10)----------------------730-830 Mev/c ------------- if( pm<=cut[10][0] ){ return wpikp(10,q); } //11)----------------------830-930 Mev/c ------------- if( pm<=cut[11][0] ){ return wpikp(11,q); } //12)----------------------930-1030 Mev/c ------------- if( pm<=cut[12][0] ){ return wpikp(12,q); } return fPcode; } //----------------------------------------------------------- void AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi, Float_t klo,Float_t khi,Float_t plo,Float_t phi) { cut[n][0]=pm; cut[n][1]=pilo; cut[n][2]=pihi; cut[n][3]=klo; cut[n][4]=khi; cut[n][5]=plo; cut[n][6]=phi; return ; } //------------------------------------------------------------ void AliITSPid::SetVec(Int_t ntrack,TVector info) { TClonesArray& arr=*trs; new( arr[ntrack] ) TVector(info); } //----------------------------------------------------------- TVector* AliITSPid::GetVec(Int_t ntrack) { TClonesArray& arr=*trs; return (TVector*)arr[ntrack]; } //----------------------------------------------------------- void AliITSPid::SetEdep(Int_t track,Float_t Edep) { TVector xx(0,11); if( ((TVector*)trs->At(track))->IsValid() ) {TVector yy( *((TVector*)trs->At(track)) );xx=yy; } Int_t j=(Int_t)xx(0); if(j>4)return; xx(++j)=Edep;xx(0)=j; TClonesArray &arr=*trs; new(arr[track])TVector(xx); } //----------------------------------------------------------- void AliITSPid::SetPmom(Int_t track,Float_t Pmom) { TVector xx(0,11); if( ((TVector*)trs->At(track))->IsValid() ) {TVector yy( *((TVector*)trs->At(track)) );xx=yy; } xx(10)=Pmom; TClonesArray &arr=*trs; new(arr[track])TVector(xx); } //----------------------------------------------------------- void AliITSPid::SetPcod(Int_t track,Int_t partcode) { TVector xx(0,11); if( ((TVector*)trs->At(track))->IsValid() ) {TVector yy( *((TVector*)trs->At(track)) );xx=yy; } if(xx(11)==0) {xx(11)=partcode; mxtrs++; TClonesArray &arr=*trs; new(arr[track])TVector(xx); } } //----------------------------------------------------------- void AliITSPid::Print(Int_t track) {cout<At(track))->IsValid() ) {TVector xx( *((TVector*)trs->At(track)) ); xx.Print(); } else {cout<<"No data for track "<GetEntries()==0){cout<<"No entries in TAB"<GetEntries();i++) { TVector xx( *((TVector*)trs->At(i)) ); if( xx.IsValid() && xx(0)>0 ) { TVector xx( *((TVector*)trs->At(i)) ); if(xx(0)>=2) { // 1)Calculate Qtrm xx(6)=(this->qtrm(i)); }else{ xx(6)=xx(1); } // 2)Calculate Wpi,Wk,Wp this->GetPcode(xx(6),xx(10)/1000.); xx(7)=GetWpi(); xx(8)=GetWk(); xx(9)=GetWp(); // 3)Print table if(xx(0)>0){ // cout<7){ cout.width(7);cout.precision(5);cout<GetEntries();i++){ TVector xx(0,11); TClonesArray &arr=*trs; new(arr[i])TVector(xx); } } //----------------------------------------------------------- AliITSPid::AliITSPid(Int_t ntrack) { fSigmin=0.01; trs = new TClonesArray("TVector",ntrack); TClonesArray &arr=*trs; for(Int_t i=0;i