From 97a4d5386a1e8114730e7a1af54cd9f135ee8ccd Mon Sep 17 00:00:00 2001 From: dibari Date: Thu, 19 Mar 2009 12:02:14 +0000 Subject: [PATCH] Bug in magnetic field solved. IMPORTANT change: the first two fields of GetHMPIDtrk are the x an y of the extrapolated track to the PC (before was at the RADIATOR). This change is to minimize the infos in ESD and to optimize the code for new mag field. --- HMPID/AliHMPIDRecon.cxx | 6 +++--- HMPID/AliHMPIDRecon.h | 2 +- HMPID/AliHMPIDTracker.cxx | 8 ++++---- HMPID/AliHMPIDtrack.cxx | 41 ++++++++------------------------------- HMPID/AliHMPIDtrack.h | 5 ++--- HMPID/HESDfromKin.C | 7 +++++-- 6 files changed, 23 insertions(+), 46 deletions(-) diff --git a/HMPID/AliHMPIDRecon.cxx b/HMPID/AliHMPIDRecon.cxx index 3bbaf6d63f5..0e832ebba7e 100644 --- a/HMPID/AliHMPIDRecon.cxx +++ b/HMPID/AliHMPIDRecon.cxx @@ -79,7 +79,7 @@ void AliHMPIDRecon::DeleteVars()const delete [] fPhotWei; } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean) +void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean,Float_t xRa,Float_t yRa) { // Pattern recognition method based on Hough transform // Arguments: pTrk - track for which Ckov angle is to be found @@ -95,8 +95,8 @@ void AliHMPIDRecon::CkovAngle(AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t inde InitVars(nClusTot); - Float_t xRa,yRa,th,ph; - pTrk->GetHMPIDtrk(xRa,yRa,th,ph); //initialize this track: th and ph angles at middle of RAD + Float_t xPc,yPc,th,ph; + pTrk->GetHMPIDtrk(xPc,yPc,th,ph); //initialize this track: th and ph angles at middle of RAD SetTrack(xRa,yRa,th,ph); fParam->SetRefIdx(nmean); diff --git a/HMPID/AliHMPIDRecon.h b/HMPID/AliHMPIDRecon.h index 705a10aab84..986507a2181 100644 --- a/HMPID/AliHMPIDRecon.h +++ b/HMPID/AliHMPIDRecon.h @@ -28,7 +28,7 @@ public : void InitVars (Int_t n); //init space for variables void DeleteVars ()const; //delete variables - void CkovAngle (AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean ); //reconstructed Theta Cerenkov + void CkovAngle (AliESDtrack *pTrk,TClonesArray *pCluLst,Int_t index,Double_t nmean,Float_t xRa,Float_t yRa );//reconstructed Theta Cerenkov Bool_t FindPhotCkov (Double_t cluX,Double_t cluY,Double_t &thetaCer,Double_t &phiCer ); //find ckov angle for single photon candidate Double_t FindRingCkov (Int_t iNclus ); //best ckov for ring formed by found photon candidates void FindRingGeom (Double_t ckovAng,Int_t level=1 ); //estimated area of ring in cm^2 and portion accepted by geometry diff --git a/HMPID/AliHMPIDTracker.cxx b/HMPID/AliHMPIDTracker.cxx index 4254a582088..29ce92554a3 100644 --- a/HMPID/AliHMPIDTracker.cxx +++ b/HMPID/AliHMPIDTracker.cxx @@ -77,8 +77,8 @@ Int_t AliHMPIDTracker::IntTrkCha(Int_t ch,AliHMPIDtrack *pTrk,Float_t &xPc,Float Double_t p2[3],n2[3]; pParam->Norm(ch,n2); pParam->Point(ch,p2,AliHMPIDParam::kPc); //point & norm for entrance to PC plane - if(pTrk->Intersect(pTrk,p1,n1)==kFALSE) return -1; //try to intersect track with the middle of radiator - if(pTrk->Intersect(pTrk,p2,n2)==kFALSE) return -1; + if(pTrk->Intersect(p1,n1)==kFALSE) return -1; //try to intersect track with the middle of radiator + if(pTrk->Intersect(p2,n2)==kFALSE) return -1; pParam->Mars2LorsVec(ch,n1,theta,phi); //track angles at RAD pParam->Mars2Lors (ch,p1,xRa,yRa); //TRKxRAD position pParam->Mars2Lors (ch,p2,xPc,yPc); //TRKxPC position @@ -151,7 +151,7 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea if(ipCh<0) continue; //no intersection at all, go after next track - pTrk->SetHMPIDtrk(xRa,yRa,theta,phi); //store initial infos + pTrk->SetHMPIDtrk(xPc,yPc,theta,phi); //store initial infos pTrk->SetHMPIDcluIdx(ipCh,9999); //set chamber, index of cluster + cluster size // track intersects the chamber ipCh: find the MIP @@ -273,7 +273,7 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea } // recon.SetImpPC(xPc,yPc); //store track impact to PC - recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(ipCh),index,nmean); //search for Cerenkov angle of this track + recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(ipCh),index,nmean,xRa,yRa); //search for Cerenkov angle of this track AliHMPIDPid pID; Double_t prob[5]; diff --git a/HMPID/AliHMPIDtrack.cxx b/HMPID/AliHMPIDtrack.cxx index 16ff9751213..0168713794a 100644 --- a/HMPID/AliHMPIDtrack.cxx +++ b/HMPID/AliHMPIDtrack.cxx @@ -201,8 +201,8 @@ Bool_t AliHMPIDtrack::PropagateTo(const AliCluster3D *c) { } return kTRUE; }//PropagateTo() - -Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) const { +//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3]) const { //+++++++++++++++++++++++++++++++++++++++++ // Origin: K. Shileev (Kirill.Shileev@cern.ch) // Finds point of intersection (if exists) of the helix with the plane. @@ -216,10 +216,10 @@ Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) //estimates initial helix length up to plane Double_t s=(pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2]; Double_t dist=99999,distPrev=dist; - Double_t x[3],p[3]; + Double_t p[3],x[3]; while(TMath::Abs(dist)>0.00001){ //calculates helix at the distance s from x0 ALONG the helix - Propagate(s,x,p,bz); + Propagate(s,x,p); //distance between current helix position and plane dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2]; if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;} @@ -232,35 +232,7 @@ Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) return kTRUE; }//Intersect() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -Bool_t AliHMPIDtrack::Intersect(AliHMPIDtrack *pTrk,Double_t pnt[3], Double_t norm[3]) { - // Finds point of intersection (if exists) of the helix with the plane. - // Stores result in fX and fP. - // Arguments: planePoint,planeNorm - the plane defined by any plane's point - // and vector, normal to the plane - // Returns: kTrue if helix intersects the plane, kFALSE otherwise. - - Double_t x0[3]; pTrk->GetXYZ(x0); //get track position in MARS - Double_t dist=99999,distPrev=dist; - Double_t x[3],p[3], - pntrad= TMath::Sqrt(pnt[0]*pnt[0]+pnt[1]*pnt[1]); - while(TMath::Abs(dist)> 0.000001){//0.00001){ - //calculates helix at the distance s from x0 ALONG the helix - pTrk->PropagateTo(pntrad);pTrk->GetXYZ(x);pTrk->GetPxPyPz(p); - //distance between current helix position and plane - dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2]; - pntrad=pntrad-dist*0.7; - //Printf("--- 111 --- dist %lf",dist); - if(TMath::Abs(2.0*dist) >= TMath::Abs(distPrev)) {return kFALSE;} - distPrev=dist; - } - //on exit pnt is intersection point,norm is track vector at that point, - //Printf("--- 222 --- dist %lf",dist); -// Printf(""); - for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];} - return kTRUE; -}//Intersect() -//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -void AliHMPIDtrack::Propagate(Double_t len, Double_t x[3],Double_t p[3], Double_t bz) const { +void AliHMPIDtrack::Propagate(Double_t len, Double_t x[3],Double_t p[3]) const { //+++++++++++++++++++++++++++++++++++++++++ // Origin: K. Shileev (Kirill.Shileev@cern.ch) // Extrapolate track along simple helix in magnetic field @@ -270,6 +242,9 @@ void AliHMPIDtrack::Propagate(Double_t len, Double_t x[3],Double_t p[3], Double_ // The momentum returned for straight-line tracks is meaningless ! //+++++++++++++++++++++++++++++++++++++++++ GetXYZ(x); + Double_t bField[3]; + TGeoGlobalMagField::Instance()->Field(x,bField); + Double_t bz = bField[2]; if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks Double_t unit[3]; GetDirection(unit); x[0]+=unit[0]*len; diff --git a/HMPID/AliHMPIDtrack.h b/HMPID/AliHMPIDtrack.h index 2ae5d483689..bf115eb617f 100644 --- a/HMPID/AliHMPIDtrack.h +++ b/HMPID/AliHMPIDtrack.h @@ -29,12 +29,11 @@ public: Double_t GetPredictedChi2(const AliCluster3D *c) const; Bool_t PropagateTo(const AliCluster3D *c); Bool_t PropagateTo(Double_t xr, Double_t x0 = 8.72, Double_t rho = 5.86e-3); //Use material definition as for TOF??? - void Propagate(Double_t len,Double_t x[3],Double_t p[3],Double_t bz) const; //HMPID method moved from AliExternalTrackParam + void Propagate(Double_t len,Double_t x[3],Double_t p[3]) const; //HMPID method moved from AliExternalTrackParam Bool_t PropagateToR(Double_t r,Double_t step); Bool_t Rotate(Double_t alpha, Bool_t absolute); Int_t GetProlongation(Double_t xk, Double_t &y, Double_t &z); - Bool_t Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) const; //HMPID method moved from AliExternalTrackParam - Bool_t Intersect(AliHMPIDtrack *pTrk,Double_t pnt[3], Double_t norm[3]) ; //just for test + Bool_t Intersect(Double_t pnt[3], Double_t norm[3]) const; //HMPID method moved from AliExternalTrackParam Bool_t Update(const AliHMPIDCluster *pClu, Double_t chi2, Int_t index); protected: diff --git a/HMPID/HESDfromKin.C b/HMPID/HESDfromKin.C index effd7d04410..4e80cce8143 100644 --- a/HMPID/HESDfromKin.C +++ b/HMPID/HESDfromKin.C @@ -73,6 +73,7 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck) if(!pStack->IsPhysicalPrimary(i)) continue; TParticle *pTrack=pStack->Particle(i); if(pTrack->GetPDG()->Charge()==0) continue; + Printf("track n. %i",i); AliESDtrack trk(pTrack); Float_t xPc,yPc,xRa,yRa,thRa,phRa; Int_t iCh=pTracker.IntTrkCha(&trk,xPc,yPc,xRa,yRa,thRa,phRa); //get chamber intersected by this track @@ -87,7 +88,7 @@ void SimEsd(AliLoader *pHL,AliESDEvent *pEsd,Bool_t htaCheck) if(phRa<0) phRa += TMath::TwoPi(); // to be verified - trk.SetHMPIDtrk(xRa,yRa,thRa,phRa); //store initial infos + trk.SetHMPIDtrk(xPc,yPc,thRa,phRa); //store initial infos pEsd->AddTrack(&trk); Int_t status; @@ -322,7 +323,9 @@ Bool_t InitGRP() { else { Printf(Form("Cannot determine the beam type from %s, assume no LHC magnet field",beamtype)); } - + Printf("------------------------------"); + Printf(" Summary for B: %s",s.Data()); + Printf("------------------------------"); AliMagF* fld = new AliMagF("MagneticFieldMap", s.Data(), 2, fcL3, fcDip, 10., map, path, btype,beamenergy); TGeoGlobalMagField::Instance()->SetField( fld ); -- 2.43.0