]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/TPCLib/AliHLTTPCConfMapFit.cxx
- fixing small bug in the cluster analyser
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCConfMapFit.cxx
index 2ddeefc50f53ea3b1c594343e732e813d28384a8..dc4f1ac9761fc53e37ed25a13436dfb44ecd9ff2 100644 (file)
@@ -2,11 +2,12 @@
 // Original: AliHLTConfMapFit.cxx,v 1.14 2005/06/14 10:55:21 cvetan 
 
 /**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * This file is property of and copyright by the ALICE HLT Project        * 
+ * ALICE Experiment at CERN, All rights reserved.                         *
  *                                                                        *
- * Authors: Anders Vestbo                                                 *
- *          Matthias Richter <Matthias.Richter@ift.uib.no>                *
- *          for The ALICE Off-line Project.                               *
+ * Primary Authors: Anders Vestbo, maintained by                          *
+ *                  Matthias Richter <Matthias.Richter@ift.uib.no>        *
+ *                  for The ALICE HLT Project.                            *
  *                                                                        *
  * Permission to use, copy, modify and distribute this software and its   *
  * documentation strictly for non-commercial purposes is hereby granted   *
     @brief  Fit class for conformal mapping tracking.
 */
 
+// see header file for class documentation                                   //
+// or                                                                        //
+// refer to README to build package                                          //
+// or                                                                        //
+// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt                          //
+
 #include "AliHLTTPCRootTypes.h"
 #include "AliHLTTPCLogging.h"
 #include "AliHLTTPCVertex.h"
@@ -35,8 +42,7 @@
 using namespace std;
 #endif
 
-ClassImp(AliHLTTPCConfMapFit)
-
+ClassImp(AliHLTTPCConfMapFit);
 
 AliHLTTPCConfMapFit::AliHLTTPCConfMapFit()
   :
@@ -55,22 +61,6 @@ AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(AliHLTTPCConfMapTrack *track,AliHLTTPCV
   //constructor
 }
 
-AliHLTTPCConfMapFit::AliHLTTPCConfMapFit(const AliHLTTPCConfMapFit&)
-  :
-  fTrack(NULL),
-  fVertex(NULL)
-{
-  // dummy copy constructor
-  //HLTFatal("copy constructor untested");
-}
-
-AliHLTTPCConfMapFit& AliHLTTPCConfMapFit::operator=(const AliHLTTPCConfMapFit&)
-{ 
-  // dummy assignment operator
-  //HLTFatal("assignment operator untested");
-  return *this;
-}
-
 AliHLTTPCConfMapFit::~AliHLTTPCConfMapFit()
 {
   // destructor
@@ -94,8 +84,6 @@ Int_t AliHLTTPCConfMapFit::FitHelix()
   return 0;
 }
 
-// #### -B0-CHANGE-START == JMT
-
 Int_t AliHLTTPCConfMapFit::FitStraightLine() {
     //fit the straight line 
     if(FitLineXY()) {
@@ -111,8 +99,6 @@ Int_t AliHLTTPCConfMapFit::FitStraightLine() {
     return 0;
 }
 
-// #### -B0-CHANGE-END == JMT
-
 Int_t AliHLTTPCConfMapFit::FitCircle()
 {
   //-----------------------------------------------------------------
@@ -129,7 +115,7 @@ Int_t AliHLTTPCConfMapFit::FitCircle()
   Double_t xav   = 0.0 ;
   Double_t yav   = 0.0 ;
   
-  Int_t num_of_hits = fTrack->GetNumberOfPoints();
+  Int_t numOfHits = fTrack->GetNumberOfPoints();
   //
   //     Loop over hits calculating average
   Int_t co=0;
@@ -143,9 +129,9 @@ Int_t AliHLTTPCConfMapFit::FitCircle()
       xav       += cHit->GetXYWeight() * cHit->GetX() ;
       yav       += cHit->GetXYWeight() * cHit->GetY() ;
     }
-  if(co!=num_of_hits) 
+  if(co!=numOfHits) 
     LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitCircle","TrackFit")<<AliHLTTPCLog::kDec<<
-      "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
+      "Mismatch of hits. Counter: "<<co<<" nHits: "<<numOfHits<<ENDLOG;
   if (fTrack->ComesFromMainVertex() == true)
     {    
       wsum += fVertex->GetXYWeight() ;
@@ -246,7 +232,7 @@ Int_t AliHLTTPCConfMapFit::FitCircle()
 
   Double_t xixi, yiyi, riri, wiriri, xold, yold ;
   
-  //for (hit_counter=0; hit_counter<num_of_hits; hit_counter++) 
+  //for (hit_counter=0; hit_counter<numOfHits; hit_counter++) 
   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
     { 
       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);  
@@ -311,7 +297,7 @@ Int_t AliHLTTPCConfMapFit::FitCircle()
   rrrrav  = rrrrav / wsum ;
   xyav    = xyav   / wsum ;
 
-  Int_t const ntry = 5 ;
+  const Int_t ntry = 5 ;
 //
 //-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
 //-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
@@ -428,6 +414,9 @@ Int_t AliHLTTPCConfMapFit::FitCircle()
   fTrack->SetCenterX(acent);
   fTrack->SetCenterY(bcent);
 
+  fTrack->SetY0err(0.03);
+  //set error for pT and Y. psi, Z and Tgl are set.
+
   //
 //    Get errors from fast fit
 //
@@ -457,7 +446,7 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
   Double_t radius = (Double_t)(fTrack->GetPt() / AliHLTTPCTransform::GetBFieldValue() ) ;
 
   //TObjArray *hits = fTrack->GetHits();
-  //Int_t num_of_hits = fTrack->GetNumberOfPoints();
+  //Int_t numOfHits = fTrack->GetNumberOfPoints();
 
   if (0)// fTrack->ComesFromMainVertex() == true ) 
     {
@@ -473,15 +462,15 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
     }
   
   Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
-  Double_t total_s ;
+  Double_t totalS ;
   
   if ( fabs(localPsi) < 1. ) 
     {
-      total_s = 2.0 * radius * asin ( localPsi ) ;
+      totalS = 2.0 * radius * asin ( localPsi ) ;
     } 
   else 
     { 
-      total_s = 2.0 * radius * AliHLTTPCTransform::Pi() ;
+      totalS = 2.0 * radius * AliHLTTPCTransform::Pi() ;
     } 
   
   AliHLTTPCConfMapPoint *previousHit = NULL;
@@ -491,7 +480,7 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
   //for ( startLoop() ; done() ; nextHit() ) {
   Double_t dpsi,s;
 
-  //  for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
+  //  for(hit_counter=0; hit_counter<numOfHits; hit_counter++)
   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
     {
       // AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
@@ -507,8 +496,8 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
          cHit->SetS(s);
        }
       else
-       cHit->SetS(total_s);
-      //       cHit->s = total_s ;
+       cHit->SetS(totalS);
+      //       cHit->s = totalS ;
     
       sum += cHit->GetZWeight() ;
       ss  += cHit->GetZWeight() * cHit->GetS() ;
@@ -540,7 +529,7 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
   chi2 = 0.;
   Double_t r1 ;
   
-  //for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
+  //for(hit_counter=0; hit_counter<numOfHits; hit_counter++)
   for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
     {
       //AliHLTTPCConfMapPoint *cHit = (AliHLTTPCConfMapPoint*)hits->At(hit_counter);
@@ -561,13 +550,38 @@ Int_t AliHLTTPCConfMapFit::FitLine ( )
   
   fTrack->SetTglerr(dtanl);
   fTrack->SetZ0err(dz0);
+
+    //The calculation of pT comes from "Perticle Physics Booklet": 
+  //24.11 Measurement of particle momenta in a uniform magnetic field.
+  //(dk)^2=(dk_res)^2 + (dk_ms)^2
+  //for now k_ms is 0. Need to find length of track in 3D. 
+
+  Double_t lengthXY = fTrack->GetLengthXY();
+  //Double_t lengthTot = fTrack->GetLengthTot();
+  //Double_t beta = fTrack->GetP()/TMath::Sqrt((fTrack->GetP()*fTrack->GetP())+(0.13957*0.13957));
+  //Double_t lambda = TMath::ATan(fTrack->GetTgl());
+  Double_t lengthXY2 = lengthXY*lengthXY;
+  Int_t nCluster4 = fTrack->GetNHits()+4;
+
+  Double_t Kres = 0.03/lengthXY2;
+  Kres = Kres * TMath::Sqrt(720/nCluster4);
+
+  //Double_t Kres = (0.03/(lengthXY*lengthXY))*TMath::Sqrt(720/(fTrack->GetNHits()+4));
+
+  //Double_t d = lengthTot*fTrack->GetP()*beta*TMath::Cos(lambda)*TMath::Cos(lambda);
+  //Double_t Kms = (0.016/d)*TMath::Sqrt(lengthTot/24.0);
+  Double_t Kms = 0.0;
+
+  Double_t KTot = TMath::Sqrt((Kres * Kres) + (Kms * Kms));
   
+  Double_t Pterr = (1/(0.3*AliHLTTPCTransform::GetBField()))*KTot;
+
+  fTrack->SetPterr(Pterr);
+
   return 0 ;
 } 
 
 
-// #### -B0-CHANGE-START == JMT
-
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 //    Straight Line Fit  in x-y plane
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -579,18 +593,18 @@ Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
     // with y = a' + bx' , x' = x - <x>
     // -----------------------------------------------------------------------------
 
-    Double_t S = 0.;
-    Double_t Sx = 0.;
+    Double_t s = 0.;
+    Double_t sx = 0.;
 
-    Double_t SPrime = 0.;
-    Double_t SxPrime = 0.;
-    Double_t SxxPrime = 0.;
-    Double_t SyPrime = 0.;
-    Double_t SxyPrime = 0.;
+    Double_t sPrime = 0.;
+    Double_t sxPrime = 0.;
+    Double_t sxxPrime = 0.;
+    Double_t syPrime = 0.;
+    Double_t sxyPrime = 0.;
 
     Double_t chi2 = 0.;
 
-    Int_t num_of_hits = fTrack->GetNumberOfPoints();
+    Int_t numOfHits = fTrack->GetNumberOfPoints();
 
     Int_t co=0;
     
@@ -601,28 +615,28 @@ Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
        // ** maybe not necessary, already done in ConfMapPoint
        currentHit->SetXYWeight( 1./ (Double_t)(currentHit->GetXerr()*currentHit->GetXerr() + currentHit->GetYerr()*currentHit->GetYerr()) );
        // **
-       S   += currentHit->GetXYWeight();
-       Sx  += currentHit->GetXYWeight() * currentHit->GetX();
+       s   += currentHit->GetXYWeight();
+       sx  += currentHit->GetXYWeight() * currentHit->GetX();
     }   
     
-    if(co!=num_of_hits) 
-       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
+    if(co!=numOfHits) 
+       LOG(AliHLTTPCLog::kError,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Mismatch of hits. Counter: "<<co<<" nHits: "<<numOfHits<<ENDLOG;
 
-    Double_t xav = (Double_t)Sx / S;
+    Double_t xav = (Double_t)sx / s;
     
     // Calculate weighted means
     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
        AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 
        Double_t xPrime =  currentHit->GetX() - xav;
-       SPrime   += currentHit->GetXYWeight();
-       SxPrime  += currentHit->GetXYWeight() * xPrime;
-       SxxPrime += currentHit->GetXYWeight() * xPrime * xPrime;
-       SyPrime  += currentHit->GetXYWeight() * currentHit->GetY();
-       SxyPrime += currentHit->GetXYWeight() * xPrime * currentHit->GetY();
+       sPrime   += currentHit->GetXYWeight();
+       sxPrime  += currentHit->GetXYWeight() * xPrime;
+       sxxPrime += currentHit->GetXYWeight() * xPrime * xPrime;
+       syPrime  += currentHit->GetXYWeight() * currentHit->GetY();
+       sxyPrime += currentHit->GetXYWeight() * xPrime * currentHit->GetY();
     }
 
-    Double_t det = SPrime*SxxPrime + SxPrime*SxPrime;
+    Double_t det = sPrime*sxxPrime + sxPrime*sxPrime;
 
     if (fabs(det) < 1e-20) { 
        LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineXY","TrackFit") << "Determinant == 0" << ENDLOG;  
@@ -631,11 +645,11 @@ Int_t AliHLTTPCConfMapFit::FitLineXY ( ){
        return -1 ;
     }
 
-    Double_t b   = (Double_t)(SPrime*SxyPrime - SxPrime*SyPrime) / det;        // line parameter b
-    Double_t aPrime   = (Double_t)(SxxPrime*SyPrime - SxPrime*SxyPrime) / det; // line parameter a
+    Double_t b   = (Double_t)(sPrime*sxyPrime - sxPrime*syPrime) / det;        // line parameter b
+    Double_t aPrime   = (Double_t)(sxxPrime*syPrime - sxPrime*sxyPrime) / det; // line parameter a
 
-    Double_t sigma2b = (Double_t)1. / SxxPrime;
-    //-- Double_t sigma2aprime = (Double_t)1. /SPrime;
+    Double_t sigma2b = (Double_t)1. / sxxPrime;
+    //-- Double_t sigma2aprime = (Double_t)1. /sPrime;
 
     // Get gradient angle psi of line in xy plane
     Double_t psi  = (Double_t) atan(b) ; 
@@ -689,14 +703,18 @@ Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
     Double_t S = 0.;
     Double_t Ss = 0.;
 
-    Double_t SPrime = 0.;
-    Double_t SsPrime = 0.;
-    Double_t SssPrime = 0.;
-    Double_t SzPrime = 0.;
-    Double_t SszPrime = 0.;
+    Double_t sPrime = 0.;
+    Double_t ssPrime = 0.;
+    Double_t sssPrime = 0.;
+    Double_t szPrime = 0.;
+    Double_t sszPrime = 0.;
 
     Double_t chi2 = 0.;
-    Double_t s = 0.;
+
+    // Matthias 16.10.2007
+    // what's that!!! local variables 's' and 'S'
+    // change Double_t s = 0.; -> slength
+    Double_t slength = 0.;
 
     AliHLTTPCConfMapPoint *previousHit = NULL;
   
@@ -707,15 +725,15 @@ Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
        if(currentHit != fTrack->GetFirstHit()) {
            Double_t dx = currentHit->GetX() - previousHit->GetX() ;
            Double_t dy = currentHit->GetY() - previousHit->GetY() ;
-           s = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy );
+           slength = previousHit->GetS() - (Double_t)sqrt ( dx*dx + dy*dy );
        }
        else{
            Double_t dx = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetX() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetX(); 
            Double_t dy = ((AliHLTTPCConfMapPoint *)fTrack->GetFirstHit())->GetY() - ((AliHLTTPCConfMapPoint *)fTrack->GetLastHit())->GetY();
-           s = (Double_t)sqrt ( dx*dx + dy*dy );
+           slength = (Double_t)sqrt ( dx*dx + dy*dy );
        }
 
-       currentHit->SetS(s);
+       currentHit->SetS(slength);
 
        S   += currentHit->GetZWeight();
        Ss  += currentHit->GetZWeight() * currentHit->GetS();
@@ -727,15 +745,18 @@ Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
     for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit()) {
        AliHLTTPCConfMapPoint *currentHit = (AliHLTTPCConfMapPoint*)fTrack->GetCurrentHit();
 
-       Double_t sPrime =  currentHit->GetS() - sav;
-       SPrime   += currentHit->GetZWeight();
-       SsPrime  += currentHit->GetZWeight() * sPrime;
-       SssPrime += currentHit->GetZWeight() * sPrime * sPrime;
-       SzPrime  += currentHit->GetZWeight() * currentHit->GetZ();
-       SszPrime += currentHit->GetZWeight() * sPrime * currentHit->GetZ();
+       // Matthias 20.05.2008
+       // here was a shadowed variable, sPrime is formerly defined
+       // renamed it to lsPrime ('local')
+       Double_t lsPrime =  currentHit->GetS() - sav;
+       lsPrime   += currentHit->GetZWeight();
+       ssPrime  += currentHit->GetZWeight() * lsPrime;
+       sssPrime += currentHit->GetZWeight() * lsPrime * lsPrime;
+       szPrime  += currentHit->GetZWeight() * currentHit->GetZ();
+       sszPrime += currentHit->GetZWeight() * lsPrime * currentHit->GetZ();
     }
 
-    Double_t det = SPrime*SssPrime + SsPrime*SsPrime;
+    Double_t det = sPrime*sssPrime + ssPrime*ssPrime;
 
     if (fabs(det) < 1e-20) { 
        LOG(AliHLTTPCLog::kDebug,"AliHLTTPCConfMapFit::FitLineSZ","TrackFit") << "Determinant == 0" << ENDLOG;  
@@ -744,13 +765,13 @@ Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
        return -1 ;
     }
 
-    Double_t b   = (Double_t)(SPrime*SszPrime - SsPrime*SzPrime) / det;        // line parameter b
-    Double_t aPrime   = (Double_t)(SssPrime*SzPrime - SsPrime*SszPrime) / det; // line parameter a
+    Double_t b   = (Double_t)(sPrime*sszPrime - ssPrime*szPrime) / det;        // line parameter b
+    Double_t aPrime   = (Double_t)(sssPrime*szPrime - ssPrime*sszPrime) / det; // line parameter a
 
     Double_t a = aPrime - b*sav;
 
-    Double_t sigma2b = (Double_t) 1. / SssPrime;
-    Double_t sigma2aprime = (Double_t) 1. /SPrime;
+    Double_t sigma2b = (Double_t) 1. / sssPrime;
+    Double_t sigma2aprime = (Double_t) 1. /sPrime;
 
     Double_t sigma2a = sigma2aprime + sav*sav * sigma2b*sigma2b;
 
@@ -771,4 +792,3 @@ Int_t AliHLTTPCConfMapFit::FitLineSZ ( ){
     return 0;
 }
 
-// #### -B0-CHANGE-END == JMT