Bug in magnetic field solved.
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Mar 2009 12:02:14 +0000 (12:02 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Mar 2009 12:02:14 +0000 (12:02 +0000)
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
HMPID/AliHMPIDRecon.h
HMPID/AliHMPIDTracker.cxx
HMPID/AliHMPIDtrack.cxx
HMPID/AliHMPIDtrack.h
HMPID/HESDfromKin.C

index 3bbaf6d..0e832eb 100644 (file)
@@ -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);
index 705a10a..986507a 100644 (file)
@@ -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
index 4254a58..29ce925 100644 (file)
@@ -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];
index 16ff975..0168713 100644 (file)
@@ -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;   
index 2ae5d48..bf115eb 100644 (file)
@@ -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:
index effd7d0..4e80cce 100644 (file)
@@ -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 );