X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCPid.cxx;h=6c713df7392dde82f647b33809eec0e6340eb714;hb=c0b978f01970d3f0f0fa200479c81ed319e59b64;hp=d9394395283790e3a2d3a26ca5490e585c67ad83;hpb=9dec0e8165300453ea4ccc7238974eff8c414bee;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCPid.cxx b/TPC/AliTPCPid.cxx index d9394395283..6c713df7392 100644 --- a/TPC/AliTPCPid.cxx +++ b/TPC/AliTPCPid.cxx @@ -13,32 +13,45 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.3 2003/03/04 07:36:15 kowal2 -Logs added - -*/ +/* $Id$ */ #include "AliTPCPid.h" #include "TMath.h" +//#include +#include +#include //#include "AliITSIOTrack.h" +#include "AliKalmanTrack.h" #include "Riostream.h" ClassImp(AliTPCPid) // Correction 13.01.2003 Z.S.,Dubna // 22.01.2003 //------------------------------------------------------------ -Float_t AliTPCPid::qcorr(Float_t xc) +// pid class by B. Batyunya +// stupid corrections by M.K. +// +//________________________________________________________ + AliTPCPid::AliTPCPid( const AliTPCPid& r):TObject(r) +{ + // dummy copy constructor +} +Float_t AliTPCPid::Qcorr(Float_t xc) { + // + // charge correction + // 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; +return fqtot/fcorr; } //----------------------------------------------------------- -Float_t AliTPCPid::qtrm(Int_t track) +Float_t AliTPCPid::Qtrm(Int_t track) const { + // + // dummy comment (Boris!!!) + // TVector q(*( this->GetVec(track) )); Int_t ml=(Int_t)q(0); if(ml<1)return 0.; @@ -62,9 +75,9 @@ Float_t AliTPCPid::qtrm(Int_t track) } for(Int_t i=0;iSetVec(track,q); return (q(6)); -}// --- End AliTPCPid::qtrm --- +}// --- End AliTPCPid::Qtrm --- -Float_t AliTPCPid::qtrm(Float_t qarr[6],Int_t narr) +Float_t AliTPCPid::Qtrm(Float_t qarr[6],Int_t narr) { //.................. Float_t q[6],qm,qmin; @@ -87,18 +100,21 @@ Float_t AliTPCPid::qtrm(Float_t qarr[6],Int_t narr) break; } return qm; -}// --- End AliTPCPid::qtrm --- +}// --- End AliTPCPid::Qtrm --- -Int_t AliTPCPid::wpik(Int_t nc,Float_t q) +Int_t AliTPCPid::Wpik(Int_t nc,Float_t q) { + // + // pi-k + // 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]; + qmpi =fcut[nc][1]; + sigpi=fcut[nc][2]; + qmk =fcut[nc][3]; + sigk =fcut[nc][4]; + appi = faprob[0][nc-5]; + apk = faprob[1][nc-5]; if( !fSilent ){ cout<<"qmpi,sigpi,qmk,sigk="<0.){ ppi=appi*TMath::Gaus(q,qmpi,sigpi)/dn; pk = apk*TMath::Gaus(q,qmk, sigk )/dn; - }else{fWpi=1;return pion();} + }else{fWpi=1;return Pion();} Float_t rpik=ppi/(pk+0.0000001); if( !fSilent ) - cout<<"q,dqpi,dqk, wpik: ppi,pk,rpik=" + cout<<"q,dqpi,dqk, Wpik: ppi,pk,rpik=" <ppi){return kaon();}else{return pion();} + if( pk>ppi){return Kaon();}else{return Pion();} } //----------------------------------------------------------- -//inline -Int_t AliTPCPid::wpikp(Int_t nc,Float_t q) +Int_t AliTPCPid::Wpikp(Int_t nc,Float_t q) { + // + // pi-k-p + // 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]; + qmpi =fcut[nc][1]; + sigpi=fcut[nc][2]; + qmk =fcut[nc][3]; + sigk =fcut[nc][4]; + qmp =fcut[nc][5]; + sigp =fcut[nc][6]; //appi = apk = app = 1.; - appi = aprob[0][nc-5]; - apk = aprob[1][nc-5]; - app = aprob[2][nc-5]; + appi = faprob[0][nc-5]; + apk = faprob[1][nc-5]; + app = faprob[2][nc-5]; Double_t dn= appi*TMath::Gaus(q,qmpi,sigpi) +apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp); @@ -146,112 +164,126 @@ Int_t AliTPCPid::wpikp(Int_t nc,Float_t q) pk = apk*TMath::Gaus(q,qmk, sigk )/dn; pp = app*TMath::Gaus(q,qmp, sigp )/dn; } - else{fWpi=1;return pion();} + else{fWpi=1;return Pion();} fWp=pp; fWpi=ppi; fWk=pk; if( !fSilent ){ -cout<<" wpikp: mid,sig pi,k,p="<pk&&ppi>pp ) { return pion(); } - if(pk>pp){return kaon();}else{return proton();} + if( ppi>pk&&ppi>pp ) { return Pion(); } + if(pk>pp){return Kaon();}else{return Proton();} } //----------------------------------------------------------- -Int_t AliTPCPid::GetPcode(TClonesArray* rps,Float_t pm) +Int_t AliTPCPid::GetPcode(TClonesArray* /*rps*/,Float_t /*pm*/) const { + // + // dummy ??? + // return 0; } //----------------------------------------------------------- Int_t AliTPCPid::GetPcode(AliTPCtrack *track) { + // + // get particle code + // 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 pt1=TMath::Abs(par[4]); + Float_t mom=1./(pt1*TMath::Cos(lam)); Float_t dedx=track->GetdEdx(); Int_t pcode=GetPcode(dedx/40.,mom); cout<<"TPCtrack 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 pt1=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="<GetPIDsignal(); + cout<<"lam,pt1,mom,dedx="<Eval(pm) ){fWpi=1.; return pion(); } - else {fWk =1.; return kaon();} } + if ( pm<=fcut[2][0] ) + { if( qEval(pm) ){fWpi=1.; return Pion(); } + else {fWk =1.; return Kaon();} } //3)----------------------200-400 Mev/c ------------- - if ( pm<=cut[3][0] ) + if ( pm<=fcut[3][0] ) if( qEval(pm) ) - { fWpi=1.;return pion(); } + { fWpi=1.;return Pion(); } else { if ( q<=fCutPr->Eval(pm) ) - {fWk=1.;return kaon();} - else {fWp=1.;return proton();} + {fWk=1.;return Kaon();} + else {fWp=1.;return Proton();} } //4)----------------------400-450 Mev/c ------------- - if ( pm<=cut[4][0] ) + if ( pm<=fcut[4][0] ) if( qEval(pm) ) - { fWpi=1.;return pion(); } + { fWpi=1.;return Pion(); } else { if( qEval(pm) ) - {fWk=1.;return kaon();} else {fWp=1.;return proton(); } + {fWk=1.;return Kaon();} else {fWp=1.;return Proton(); } } //5)----------------------450-500 Mev/c ------------- - if ( pm<=cut[5][0] ) + if ( pm<=fcut[5][0] ) if ( q>fCutPr->Eval(pm) ) - {fWp=1.;return proton();} else {return wpik(5,q);}; + {fWp=1.;return Proton();} else {return Wpik(5,q);}; //6)----------------------500-550 Mev/c ------------- - if ( pm<=cut[6][0] ) + if ( pm<=fcut[6][0] ) if ( q>fCutPr->Eval(pm) ) - {fWp=1.;return proton();} else {return wpik(6,q);}; + {fWp=1.;return Proton();} else {return Wpik(6,q);}; //7)----------------------550-600 Mev/c ------------- - if ( pm<=cut[7][0] ) + if ( pm<=fcut[7][0] ) if ( q>fCutPr->Eval(pm) ) - {fWp=1.;return proton();} else {return wpik(7,q);}; + {fWp=1.;return Proton();} else {return Wpik(7,q);}; //8)----------------------600-650 Mev/c ------------- - if ( pm<=cut[8][0] ) - if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return proton();} - else {return wpik(8,q);}; + if ( pm<=fcut[8][0] ) + if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();} + else {return Wpik(8,q);}; //9)----------------------650-730 Mev/c ------------- - if ( pm<=cut[9][0] ) - if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return proton();} - else {return wpik(9,q);}; + if ( pm<=fcut[9][0] ) + if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();} + else {return Wpik(9,q);}; //10)----------------------730-830 Mev/c ------------- - if( pm<=cut[10][0] ) - if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return proton();} - else {return wpik(10,q);}; + if( pm<=fcut[10][0] ) + if ( q>fCutPrTune*fCutPr->Eval(pm) ){fWp=1.;return Proton();} + else {return Wpik(10,q);}; //11)----------------------830-930 Mev/c ------------- - if( pm<=cut[11][0] ){ return wpikp(11,q); } + if( pm<=fcut[11][0] ){ return Wpikp(11,q); } //12)----------------------930-1030 Mev/c ------------- - if( pm<=cut[12][0] ) - { return wpikp(12,q); }; + if( pm<=fcut[12][0] ) + { return Wpikp(12,q); }; return fPcode; } @@ -259,30 +291,42 @@ Int_t AliTPCPid::GetPcode(Float_t q,Float_t pm) void AliTPCPid::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; + // + // set cuts + // + fcut[n][0]=pm; + fcut[n][1]=pilo; + fcut[n][2]=pihi; + fcut[n][3]=klo; + fcut[n][4]=khi; + fcut[n][5]=plo; + fcut[n][6]=phi; return ; } //------------------------------------------------------------ -void AliTPCPid::SetVec(Int_t ntrack,TVector info) +void AliTPCPid::SetVec(Int_t ntrack,TVector info) const { + // + // new track vector + // TClonesArray& arr=*trs; new( arr[ntrack] ) TVector(info); } //----------------------------------------------------------- -TVector* AliTPCPid::GetVec(Int_t ntrack) +TVector* AliTPCPid::GetVec(Int_t ntrack) const { + // + // get track vector + // TClonesArray& arr=*trs; return (TVector*)arr[ntrack]; } //----------------------------------------------------------- void AliTPCPid::SetEdep(Int_t track,Float_t Edep) { + // + // energy deposit + // TVector xx(0,11); if( ((TVector*)trs->At(track))->IsValid() ) {TVector yy( *((TVector*)trs->At(track)) );xx=yy; } @@ -294,6 +338,9 @@ void AliTPCPid::SetEdep(Int_t track,Float_t Edep) //----------------------------------------------------------- void AliTPCPid::SetPmom(Int_t track,Float_t Pmom) { + // + // set part. momentum + // TVector xx(0,11); if( ((TVector*)trs->At(track))->IsValid() ) {TVector yy( *((TVector*)trs->At(track)) );xx=yy; } @@ -304,18 +351,25 @@ void AliTPCPid::SetPmom(Int_t track,Float_t Pmom) //----------------------------------------------------------- void AliTPCPid::SetPcod(Int_t track,Int_t partcode) { + // + // set part. code + // 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++; + {xx(11)=partcode; fmxtrs++; TClonesArray &arr=*trs; new(arr[track])TVector(xx); } } //----------------------------------------------------------- void AliTPCPid::Print(Int_t track) -{cout<At(track))->IsValid() ) {TVector xx( *((TVector*)trs->At(track)) ); xx.Print(); @@ -326,6 +380,9 @@ void AliTPCPid::Print(Int_t track) //----------------------------------------------------------- void AliTPCPid::Tab(void) { + // + // fill table + // if(trs->GetEntries()==0){cout<<"No entries in TAB"<GetEntries();i++) if(xx(0)>=2) { // 1)Calculate Qtrm - xx(6)=(this->qtrm(i)); + xx(6)=(this->Qtrm(i)); }else{ xx(6)=xx(1); @@ -368,6 +425,9 @@ for(Int_t i=0;iGetEntries();i++) } void AliTPCPid::Reset(void) { + // + // reset + // for(Int_t i=0;iGetEntries();i++){ TVector xx(0,11); TClonesArray &arr=*trs; @@ -380,11 +440,13 @@ AliTPCPid::AliTPCPid(Int_t ntrack) trs = new TClonesArray("TVector",ntrack); TClonesArray &arr=*trs; for(Int_t i=0;i