]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Bug fixes in the fitter. Possibility of different points selection for
authormasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Dec 2009 08:20:31 +0000 (08:20 +0000)
committermasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Dec 2009 08:20:31 +0000 (08:20 +0000)
cosmics and pp events in simultaneous fit (Ruben)

ITS/AliITSAlignMille2.cxx
ITS/AliITSAlignMille2.h
ITS/AliITSTPArrayFit.cxx
ITS/AliITSTPArrayFit.h

index c7f6a9e169233bb455e75e44072b49085e1b0ab3..cded35c2349e28210f750c1ba72e6bee23f47a3a 100644 (file)
@@ -136,7 +136,6 @@ AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInf
   fCacheMatrixCurr(kMaxITSSensID+1),
   //
   fUseGlobalDelta(kFALSE),
-  fRequirePoints(kFALSE),
   fTempExcludedModule(-1),
   //
   fDefCDBpath("local://$ALICE_ROOT/OCDB"),
@@ -162,6 +161,7 @@ AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInf
   fUseLocalYErr(kFALSE),
   fBOn(kFALSE),
   fBField(0.0),
+  fDataType(kCosmics),
   fMinPntPerSens(0),
   fBug(0),
   fMilleVersion(2),
@@ -171,15 +171,21 @@ AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInf
   for (int i=3;i--;) fSigmaFactor[i] = 1.0;
   //
   // new RS
-  for (Int_t i=0; i<6; i++) {
-    fNReqLayUp[i]=0;
-    fNReqLayDown[i]=0;
-    fNReqLay[i]=0;
+  for (int i=0;i<3;i++) {
+    
   }
-  for (Int_t i=0; i<3; i++) {
-    fNReqDetUp[i]=0;
-    fNReqDetDown[i]=0;
-    fNReqDet[i]=0;
+  for (int itp=0;itp<kNDataType;itp++) {
+    fRequirePoints[itp] = kFALSE;
+    for (Int_t i=0; i<6; i++) {
+      fNReqLayUp[itp][i]=0;
+      fNReqLayDown[itp][i]=0;
+      fNReqLay[itp][i]=0;
+    }
+    for (Int_t i=0; i<3; i++) {
+      fNReqDetUp[itp][i]=0;
+      fNReqDetDown[itp][i]=0;
+      fNReqDet[itp][i]=0;
+    }
   }
   //
   if (ProcessUserInfo(userInfo)) exit(1);
@@ -630,20 +636,34 @@ Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
          int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
          int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
          int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
-         fRequirePoints = kTRUE;
-         if (recOpt == "LAYER") {
-           if (lr<0 || lr>5) {stopped = kTRUE; break;}
-           if (hb>0) fNReqLayUp[lr] = np;
-           else if (hb<0) fNReqLayDown[lr] = np;
-           else fNReqLay[lr] = np;
+         //
+         int rtp = -1; // use for run type
+         if (nrecElems>5) {
+           TString tpstr = ((TObjString*)recArr->At(5))->GetString();
+           if ( tpstr.Contains("cosmics",TString::kIgnoreCase) ) rtp = kCosmics;
+           else if ( tpstr.Contains("collision",TString::kIgnoreCase) ) rtp = kCollision;
+           else {stopped = kTRUE; break;}
          }
-         else if (recOpt == "DETECTOR") {
-           if (lr<0 || lr>2) {stopped = kTRUE; break;}
-           if (hb>0) fNReqDetUp[lr] = np;
-           else if (hb<0) fNReqDetDown[lr] = np;
-           else fNReqDet[lr] = np;
+         //
+         int tpmn= rtp<0 ? 0 : rtp;
+         int tpmx= rtp<0 ? kNDataType-1 : rtp;
+         for (int itp=tpmn;itp<=tpmx;itp++) {
+           fRequirePoints[itp]=kTRUE;
+           if (recOpt == "LAYER") {
+             if (lr<0 || lr>5) {stopped = kTRUE; break;}
+             if (hb>0) fNReqLayUp[itp][lr]=np;
+             else if (hb<0) fNReqLayDown[itp][lr]=np;
+             else fNReqLay[itp][lr]=np;
+           }
+           else if (recOpt == "DETECTOR") {
+             if (lr<0 || lr>2) {stopped = kTRUE; break;}
+             if (hb>0) fNReqDetUp[itp][lr]=np;
+             else if (hb<0) fNReqDetDown[itp][lr]=np;
+             else fNReqDet[itp][lr]=np;
+           }
+           else {stopped = kTRUE; break;}
          }
-         else {stopped = kTRUE; break;}
+         if (stopped) break;
        }
        else {stopped = kTRUE; break;}
       }
@@ -902,7 +922,7 @@ void AliITSAlignMille2::SetCurrentModule(Int_t id)
 }
 // endpepo
 //________________________________________________________________________________________________________
-void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts) 
+void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts,Int_t runtype
 {
   // set minimum number of points in specific detector or layer
   // where = LAYER or DETECTOR
@@ -910,19 +930,23 @@ void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw,
   // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
   // nreqpts = minimum number of points of that type
   ndet--;
-  if (strstr(where,"LAYER")) {
-    if (ndet<0 || ndet>5) return;
-    if (updw>0) fNReqLayUp[ndet]=nreqpts;
-    else if (updw<0) fNReqLayDown[ndet]=nreqpts;
-    else fNReqLay[ndet]=nreqpts;
-    fRequirePoints=kTRUE;
-  }
-  else if (strstr(where,"DETECTOR")) {
-    if (ndet<0 || ndet>2) return;
-    if (updw>0) fNReqDetUp[ndet]=nreqpts;
-    else if (updw<0) fNReqDetDown[ndet]=nreqpts;
-    else fNReqDet[ndet]=nreqpts;       
-    fRequirePoints=kTRUE;
+  int tpmn= runtype<0 ? 0 : runtype;
+  int tpmx= runtype<0 ? kNDataType-1 : runtype;
+  //
+  for (int itp=tpmn;itp<=tpmx;itp++) {
+    fRequirePoints[itp]=kTRUE;
+    if (strstr(where,"LAYER")) {
+      if (ndet<0 || ndet>5) return;
+      if (updw>0) fNReqLayUp[itp][ndet]=nreqpts;
+      else if (updw<0) fNReqLayDown[itp][ndet]=nreqpts;
+      else fNReqLay[itp][ndet]=nreqpts;
+    }
+    else if (strstr(where,"DETECTOR")) {
+      if (ndet<0 || ndet>2) return;
+      if (updw>0) fNReqDetUp[itp][ndet]=nreqpts;
+      else if (updw<0) fNReqDetDown[itp][ndet]=nreqpts;
+      else fNReqDet[itp][ndet]=nreqpts;        
+    }
   }
 }
 
@@ -1371,7 +1395,7 @@ AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *at
   // >> RS
   AliTrackPoint p;
   // check points in specific places
-  if (fRequirePoints) {
+  if (fRequirePoints[fDataType]) {
     Int_t nlayup[6],nlaydown[6],nlay[6];
     Int_t ndetup[3],ndetdown[3],ndet[3];
     for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
@@ -1403,14 +1427,14 @@ AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *at
     // checks minimum values
     Bool_t isok=kTRUE;
     for (Int_t j=0; j<6; j++) {
-      if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE; 
-      if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE; 
-      if (nlay[j]<fNReqLay[j]) isok=kFALSE; 
+      if (nlayup[j]<fNReqLayUp[fDataType][j]) isok=kFALSE; 
+      if (nlaydown[j]<fNReqLayDown[fDataType][j]) isok=kFALSE; 
+      if (nlay[j]<fNReqLay[fDataType][j]) isok=kFALSE; 
     }
     for (Int_t j=0; j<3; j++) {
-      if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE; 
-      if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE; 
-      if (ndet[j]<fNReqDet[j]) isok=kFALSE; 
+      if (ndetup[j]<fNReqDetUp[fDataType][j]) isok=kFALSE; 
+      if (ndetdown[j]<fNReqDetDown[fDataType][j]) isok=kFALSE; 
+      if (ndet[j]<fNReqDet[fDataType][j]) isok=kFALSE; 
     }
     if (!isok) {
       AliDebug(2,Form("Track does not meet all location point requirements!"));
@@ -1714,16 +1738,18 @@ void AliITSAlignMille2::Print(Option_t*) const
   //
   printf("Using local Y error due to the sensor thickness: %s\n",(fUseLocalYErr && fTPAFitter) ? "ON":"OFF");
   //
-  if (fRequirePoints) printf("    Required points in tracks:\n");
-  for (Int_t i=0; i<6; i++) {
-    if (fNReqLayUp[i]>0) printf("        Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
-    if (fNReqLayDown[i]>0) printf("        Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
-    if (fNReqLay[i]>0) printf("        Layer %d : %d points \n",i+1,fNReqLay[i]);
-  }
-  for (Int_t i=0; i<3; i++) {
-    if (fNReqDetUp[i]>0) printf("        Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
-    if (fNReqDetDown[i]>0) printf("        Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
-    if (fNReqDet[i]>0) printf("        Detector %d : %d points \n",i+1,fNReqDet[i]);
+  for (int itp=0;itp<kNDataType;itp++) {
+    if (fRequirePoints[itp]) printf("    Required points in %s tracks:\n",itp==kCosmics? "cosmics" : "collisions");
+    for (Int_t i=0; i<6; i++) {
+      if (fNReqLayUp[itp][i]>0) printf("        Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[itp][i]);
+      if (fNReqLayDown[itp][i]>0) printf("        Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[itp][i]);
+      if (fNReqLay[itp][i]>0) printf("        Layer %d : %d points \n",i+1,fNReqLay[itp][i]);
+    }
+    for (Int_t i=0; i<3; i++) {
+      if (fNReqDetUp[itp][i]>0) printf("        Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[itp][i]);
+      if (fNReqDetDown[itp][i]>0) printf("        Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[itp][i]);
+      if (fNReqDet[itp][i]>0) printf("        Detector %d : %d points \n",i+1,fNReqDet[itp][i]);
+    }
   }
   //  
   printf("\n    Millepede configuration parameters:\n");
@@ -2055,8 +2081,10 @@ Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track)
   // move points if prealignment is set, sort by Yglo if required
   
   fTrack=PrepareTrack(track);
-  if (!fTrack) return -1;
-
+  if (!fTrack) {
+    RemoveHelixFitConstraint();
+    return -1;
+  }
   npts = fTrack->GetNPoints();
   if (npts>kMaxPoints) {
     AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
@@ -2066,20 +2094,25 @@ Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track)
   if (fTPAFitter) {  // use dediacted fitter
     //
     fTPAFitter->AttachPoints(fTrack);
-    if (fBOn) fTPAFitter->SetBz(fBField);
+    fTPAFitter->SetBz(fBField);
+    fTPAFitter->SetTypeCosmics(IsTypeCosmics());
     if (fInitTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
     double chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
     //
     // suppress eventual constraints to not affect fit of the next track
-    fConstrCharge = 0;
-    fConstrPT = fConstrPTErr = -1;
+    RemoveHelixFitConstraint();
     //
-    if ( chi2<0 || (fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) {
+    if ( chi2<0 || (chi2>fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) { //RRR
       AliInfo("Track fit failed! skipping this track...");
+      fTrack->Print("");
+      fTPAFitter->FitHelixCrude();
+      fTPAFitter->SetFitDone();
+      fTPAFitter->Print();
       fTPAFitter->Reset();
       fTrack = NULL;
       return -5;
     }
+    fNLocal = fTPAFitter->IsFieldON() ? 5:4; // Attantion: the fitter might have decided to work in line mode
     /*
     double *pr = fTPAFitter->GetParams();
     printf("FtPar: %+.5e  %+.5e  %+.5e  %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
@@ -3415,8 +3448,9 @@ Int_t AliITSAlignMille2::ProcessUserInfo(TList* userInfo)
     //
   }  
   //
-  objStr = (TObjString*)userInfo->FindObject("BzkGauss");
-  if (objStr) {
+  TList *bzlst = (TList*)userInfo->FindObject("BzkGauss");
+  if (bzlst && bzlst->At(0)) {
+    objStr = (TObjString*)bzlst->At(0);
     SetBField( objStr->GetString().Atof() );
     AliInfo(Form("Magentic field from UserInfo: %+.2e",GetBField()));
   }
@@ -3555,6 +3589,14 @@ Int_t AliITSAlignMille2::CacheMatricesOrig()
   return 0;
 }
 
+//________________________________________________________________________________________________________
+void AliITSAlignMille2::RemoveHelixFitConstraint()
+{
+  // suppress constraint
+  fConstrCharge = 0;
+  fConstrPT = fConstrPTErr = -1;
+}
+
 //________________________________________________________________________________________________________
 void AliITSAlignMille2::ConstrainHelixFitPT(Int_t q,Double_t pt,Double_t pterr)
 {
@@ -3579,7 +3621,7 @@ void AliITSAlignMille2::ConstrainHelixFitCurv(Int_t q,Double_t crv,Double_t crve
   }
   else {
     fConstrPT    = 1./crv*fBField*kCQConv;
-    fConstrPTErr =fConstrPT/crv*crverr;
+    fConstrPTErr = crverr>1e-10 ? fConstrPT/crv*crverr : 0.;
   }
 }
 
index 5d5c9b2e297af3ff8dcc77462e610c3c574677c1..00fbe84d2d4cb4bd5511ed690ccdd1eeed73b6da 100644 (file)
@@ -40,6 +40,7 @@ class AliITSAlignMille2: public TObject
 {
  public:
   enum {kX,kY,kZ};
+  enum {kCosmics, kCollision, kNDataType};
   enum {kNLocal=5,kMaxPoints=100,
        kNParChGeom = AliITSAlignMille2Module::kMaxParGeom,
        kNParCh     = AliITSAlignMille2Module::kMaxParTot,
@@ -124,14 +125,19 @@ class AliITSAlignMille2: public TObject
   Double_t  GetMeasLoc(Int_t dim)                                 const {return fMeasLoc[dim];}
   Int_t     GetCurrentLayer()                                     const;
   void      SetBField(Double_t b=0);
+  void      SetTypeCosmics()                                            {fDataType = kCosmics;}
+  void      SetTypeCollision()                                          {fDataType = kCollision;}
+  void      SetDataType(Int_t tp=kCosmics)                              {fDataType = tp>=0&&tp< kNDataType ? tp:kCosmics;}
   void      SetUseLocalYErrors(Bool_t v=kTRUE)                          {fUseLocalYErr = v && fTPAFitter;}
   void      SetMinPointsPerSensor( Int_t n )                            {fMinPntPerSens = n>0 ? n:0;}
   Int_t     GetMinPointsPerSensor()                               const {return fMinPntPerSens;}
   void      ConstrainHelixFitPT(  Int_t q=0,Double_t pt=-1, Double_t err=-1);
   void      ConstrainHelixFitCurv(Int_t q=0,Double_t crv=-1,Double_t crverr=-1);
+  void      RemoveHelixFitConstraint();
   Double_t  GetHelixContraintCharge()                             const {return fConstrCharge;}
   Double_t  GetHelixContraintPT()                                 const {return fConstrPT;}
   Double_t  GetHelixContraintPTErr()                              const {return fConstrPTErr;}
+  Int_t     GetDataType()                                         const {return fDataType;}
   //
   TGeoHMatrix* GetSensorOrigMatrixSID(Int_t sid)                  const;
   TGeoHMatrix* GetSensorOrigMatrixVID(Int_t vid)                  const;
@@ -173,9 +179,12 @@ class AliITSAlignMille2: public TObject
   // debug stuffs
   void       FetchCluster(int ip)                                      {fTrack->GetPoint(fCluster,ip);fCluster.SetUniqueID(ip);} 
   void       SetLocalInitParams(const Double_t *par)                   {for (int i=kNLocal;i--;) fLocalInitParam[i]=par[i];}
+  Bool_t     IsTypeCosmics()                                     const {return fDataType==kCosmics;}
+  Bool_t     IsTypeCollision()                                   const {return fDataType==kCollision;}
   Double_t  *GetMeasLoc()                                        const {return (Double_t*)fMeasLoc;}
   Double_t  *GetSigmaLoc()                                       const {return (Double_t*)fSigmaLoc;}
   Double_t   GetBField()                                         const {return fBField;}
+  Bool_t     IsFieldON()                                         const {return fBOn;}
   Double_t  *GetLocalInitParam()                                 const {return (Double_t*)fLocalInitParam;}
   Double_t  *GetLocalInitParEr()                                 const {return (Double_t*)fLocalInitParEr;}
   Double_t   GetLocalDif(int par, int coor)                      const {return fDerivativeLoc[par][coor];}
@@ -199,11 +208,11 @@ class AliITSAlignMille2: public TObject
   Int_t     IsDefined(UShort_t voluid) const {return IsVIDDefined(voluid);}
   Int_t     IsContained(UShort_t voluid) const {return IsVIDContained(voluid);}
   // moved from private to public
-  void      SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts); 
+  void      SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts,Int_t runtype=-1); 
   Bool_t    InitRiemanFit();
   void      SetMinNPtsPerTrack(Int_t pts=3)  {fMinNPtsPerTrack=pts;}
   //
-  static Bool_t    IsZero(Double_t v,Double_t threshold = 1e-16)       { return TMath::Abs(v)<threshold; }
+  static Bool_t    IsZero(Double_t v,Double_t threshold = 1e-15)       { return TMath::Abs(v)<threshold; }
   static void      SetWordBit(UInt_t word,Int_t bitID)                 { word |= (1<<bitID);}
   static void      ResetWordBit(UInt_t word,Int_t bitID)               { word &= ~(1<<bitID);}
   static Bool_t    TestWordBit(UInt_t word,Int_t bitID)                { return (Bool_t)(word&(1<<bitID));}      
@@ -348,14 +357,14 @@ class AliITSAlignMille2: public TObject
   TObjArray     fCacheMatrixOrig;               // cach for original geom matrices
   TObjArray     fCacheMatrixCurr;               // cach for prealigned geom matrices
   // >> new members
-  Bool_t        fUseGlobalDelta;  // intetpret deltas as global 
-  Bool_t        fRequirePoints;   // required points in specific layers
-  Int_t         fNReqLayUp[6];    /// number of points required in layer[n] with Y>0
-  Int_t         fNReqLayDown[6];  /// number of points required in layer[n] with Y<0
-  Int_t         fNReqLay[6];      /// number of points required in layer[n] 
-  Int_t         fNReqDetUp[3];    /// number of points required in Detector[n] with Y>0
-  Int_t         fNReqDetDown[3];  /// number of points required in Detector[n] with Y<0
-  Int_t         fNReqDet[3];      /// number of points required in Detector[n]
+  Bool_t        fUseGlobalDelta;                // intetpret deltas as global 
+  Bool_t        fRequirePoints[kNDataType];     // required points in specific layers
+  Int_t         fNReqLayUp[kNDataType][6];      // number of points required in layer[n] with Y>0
+  Int_t         fNReqLayDown[kNDataType][6];    // number of points required in layer[n] with Y<0
+  Int_t         fNReqLay[kNDataType][6];        // number of points required in layer[n] 
+  Int_t         fNReqDetUp[kNDataType][3];      // number of points required in Detector[n] with Y>0
+  Int_t         fNReqDetDown[kNDataType][3];    // number of points required in Detector[n] with Y<0
+  Int_t         fNReqDet[kNDataType][3];        // number of points required in Detector[n]
   Int_t         fTempExcludedModule; /// single module temporary excluded from initial fit
   // << new members
   //
@@ -384,6 +393,7 @@ class AliITSAlignMille2: public TObject
   Bool_t        fUseLocalYErr;                  // use local Yerror due to the sensor thickness
   Bool_t        fBOn;                           // magentic field ON
   Double_t      fBField;                        // value of magnetic field
+  Int_t         fDataType;                      // is this cosmics or collision processing?
   Int_t         fMinPntPerSens;                 // min number of points per module to vary it
   Int_t         fBug;                           /// tag for temporary bug correction
   // pepo
index 5d1402c575fe4de74f8ad06b5c08db8deac3b5b9..ec6b6bf139386d9354327c6a2bce01bf2f189c4d 100644 (file)
@@ -88,7 +88,8 @@ Double_t AliITSTPArrayFit::fgRhoLITS[AliITSTPArrayFit::kMaxLrITS] = {
 AliITSTPArrayFit::AliITSTPArrayFit() :
   fPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
   fPntLast(-1),fNPBooked(0),fParAxis(-1),fCovI(0),fChi2NDF(0),
-  fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fkAxID(0),fkAxCID(0),fCurT(0),
+  fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(1.e6),
+  fkAxID(0),fkAxCID(0),fCurT(0),
   fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
 {
   // default constructor
@@ -101,8 +102,8 @@ AliITSTPArrayFit::AliITSTPArrayFit() :
 AliITSTPArrayFit::AliITSTPArrayFit(Int_t np) :
   fPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
   fPntLast(-1),fNPBooked(np),fParAxis(-1),fCovI(0),fChi2NDF(0),
-  fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fkAxID(0),fkAxCID(0),fCurT(0),
-  fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
+  fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(2.e3),
+  fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
 {
   // constructor with booking of np points
   for (int i=kMaxParam;i--;)   fParams[i] = 0;
@@ -118,7 +119,8 @@ AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) :
   TObject(src),fPoints(src.fPoints),fParSol(0),fBz(src.fBz),
   fCharge(src.fCharge),fPntFirst(src.fPntFirst),fPntLast(src.fPntLast),fNPBooked(src.fNPBooked),
   fParAxis(src.fParAxis),fCovI(0),fChi2NDF(0),fMaxIter(20),fIter(0),fEps(0),fMass(src.fMass),
-  fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
+  fSwitch2Line(src.fSwitch2Line),fMaxRforHelix(src.fMaxRforHelix),fkAxID(0),fkAxCID(0),fCurT(0),
+  fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
 {
   // copy constructor
   InitAux();
@@ -187,6 +189,7 @@ void AliITSTPArrayFit::Reset()
   fIter = 0;
   fPntFirst=fPntLast=-1; 
   SetParAxis(-1);
+  fSwitch2Line = kFALSE;
   ResetBit(kFitDoneBit|kCovInvBit);
 }
 
@@ -946,7 +949,7 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
   Bool_t eloss = IsELossON();
   //
   int np = fPntLast - fPntFirst + 1;
-  if (np<2) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
+  if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
   //
   const float *x=fPoints->GetX(),*y=fPoints->GetY(),*z=fPoints->GetZ(),*cov=fPoints->GetCov();
   //
@@ -1005,11 +1008,11 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
   Double_t yc   = -(rhs0*mi01 + rhs1*mi11 + rhs2*mi12)/2;
   Double_t rho2 =  (rhs0*mi02 + rhs1*mi12 + rhs2*mi22);
   //
-  Double_t dcen = xc*xc + yc*yc;
-  Double_t rad = dcen - rho2;
+  Double_t rad = xc*xc + yc*yc - rho2;
   rad = (rad>fgkAlmostZero) ? (TMath::Sqrt(rad)):fgkAlmostZero;
   //
   //  printf("Rad: %+e xc: %+e yc: %+e\n",rad,xc,yc);
+
   // linear circle fit --------------------------------------------------- <<<
   //
   // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
@@ -1028,14 +1031,13 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
     SetCharge( fBz<0 ? -sqb : sqb);
   }
   //
-  dcen = TMath::Sqrt(dcen);
-  fParams[kD0] = dcen-rad;
   Double_t phi = TMath::ATan2(yc,xc);
   if (sqb<0) phi += TMath::Pi();
   if      (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
   else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
   fParams[kPhi0] = phi;  
   fParams[kR0]   = sqb<0 ? -rad:rad;  
+  fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
   //
   // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
   //
@@ -1199,6 +1201,12 @@ Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr
   //
   if (!FitHelixCrude(extQ)) return -1; // get initial estimate, ignoring the errors
   //
+  if (TMath::Abs(fParams[kR0])>fMaxRforHelix && extPT<=0) {
+    fSwitch2Line = kTRUE;
+    return FitLine();
+  }
+  //
+  //
   if (!IsCovInv()) InvertPointsCovMat();    // prepare inverted errors
   if (!fParSol) fParSol = new AliParamSolver(5);
   fParSol->SetNGlobal(5);
@@ -1225,8 +1233,8 @@ Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr
       dXYZdGlo[offs + kZ] = 0;
       //
       offs = kPhi0*3;                  // dXYZ/dPhi0
-      dXYZdGlo[offs + kX] = -rrho*sn0;
-      dXYZdGlo[offs + kY] =  rrho*cs0;
+      dXYZdGlo[offs + kX] = -rrho*sn0 + fParams[kR0]*snt;
+      dXYZdGlo[offs + kY] =  rrho*cs0 - fParams[kR0]*cst;
       dXYZdGlo[offs + kZ] = 0;
       //
       offs = kR0*3;                   // dXYZ/dR0
@@ -1263,7 +1271,11 @@ Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr
     Double_t *deltaG = fParSol->GetGlobals();
     Double_t *deltaT = fParSol->GetLocals();
     for (int ipar=5;ipar--;) fParams[ipar] -= deltaG[ipar];
-    for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst];
+    for (int ip=fPntFirst;ip<=fPntLast;ip++) {
+      fCurT[ip] = CalcParPCA(ip);
+      //      printf("iter%d : delta%2d : expl: %+e | %+e | R=%+e p0=%+e\n",iter,ip,deltaT[ip-fPntFirst], fCurT[ip],fParams[kR0],fParams[kPhi0]);
+      //      fCurT[ip] -= deltaT[ip-fPntFirst];
+    }
     iter++;
     //
     fChi2NDF = CalcChi2NDF();
@@ -1603,3 +1615,13 @@ Bool_t AliITSTPArrayFit::CalcErrorMatrix()
   //
   return kFALSE;
 }
+
+//____________________________________________________
+Double_t AliITSTPArrayFit::CalcParPCA(Int_t ipnt) const
+{
+  // get parameter for the point with least weighted distance to the point
+  const double *xyz  = GetPoint(ipnt);
+  const double *covI = GetCovI(ipnt);
+  if (IsFieldON()) return GetParPCAHelix(xyz,covI);
+  else             return GetParPCALine(xyz,covI);
+}
index 5cd4385cee450ec4f9c16ef95d31292a00ad9073..8f85913f5b0c3a71282d0c490c69c707f11c0a60 100644 (file)
@@ -65,7 +65,7 @@ class AliITSTPArrayFit : public TObject
   //
   void          SetBz(Double_t bz)                                {fBz = bz;}
   Double_t      GetBz()                                     const {return fBz;}
-  Bool_t        IsFieldON()                                 const {return TMath::Abs(fBz)>1e-5;}
+  Bool_t        IsFieldON()                                 const {return TMath::Abs(fBz)>1e-5 && !IsSwitched2Line();}
   Bool_t        IsTypeCosmics()                             const {return TestBit(kCosmicsBit);}
   Bool_t        IsTypeCollision()                           const {return !IsTypeCosmics();}
   Int_t         GetCharge()                                 const {return fCharge;}
@@ -83,6 +83,7 @@ class AliITSTPArrayFit : public TObject
   Double_t      CalcChi2NDF()                               const;
   Double_t      GetChi2NDF()                                const {return fChi2NDF;}
   Double_t      GetParPCA(const double *xyz, const double *covI=0)        const;
+  Double_t      CalcParPCA(Int_t ipnt)                      const;
   Bool_t        CalcErrorMatrix();
   //
   void          GetDResDParamsLine (Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0) const;
@@ -113,7 +114,7 @@ class AliITSTPArrayFit : public TObject
   Int_t         GetNParams()                                const {return IsFieldON() ? 5:4;}
   Bool_t        InvertPointsCovMat();
   //
-  Int_t*     GetElsId() const {return fElsId;}
+  Int_t*        GetElsId() const {return fElsId;}
   Double_t*     GetElsDR() const {return fElsDR;}
   //
   Double_t*     GetCovI(Int_t ip)                           const {return fCovI + ip*6;}
@@ -125,6 +126,7 @@ class AliITSTPArrayFit : public TObject
   Double_t      GetLineOffset(Int_t axis)                   const;
   Double_t      GetLineSlope(Int_t axis)                    const;
   //
+  Bool_t        IsSwitched2Line()                           const {return fSwitch2Line;}
   Bool_t        IsELossON()                                 const {return TestBit(kELossBit)&&IsFieldON();}
   Bool_t        IsFitDone()                                 const {return TestBit(kFitDoneBit);}
   Bool_t        IsCovInv()                                  const {return TestBit(kCovInvBit);}
@@ -134,6 +136,8 @@ class AliITSTPArrayFit : public TObject
   Double_t      GetEps()                                    const {return fEps;}
   Double_t      GetMass()                                   const {return fMass;}
   //
+  void          Switch2Line(Bool_t v=kTRUE)                       {fSwitch2Line = v;}
+  void          SetMaxRforHelix(Double_t r)                       {fMaxRforHelix = r>0 ? r : 1e9;}
   void          SetCharge(Int_t q=1)                              {fCharge = q<0 ? -1:1;}
   void          SetELossON(Bool_t v=kTRUE)                        {SetBit(kELossBit,v);}
   void          SetTypeCosmics(Bool_t v=kTRUE)                    {SetBit(kCosmicsBit,v);}
@@ -187,6 +191,8 @@ class AliITSTPArrayFit : public TObject
   Int_t     fIter;                                 // real number of iterations
   Double_t  fEps;                                  // precision
   Double_t  fMass;                                 // assumed particle mass for ELoss Calculation
+  Bool_t    fSwitch2Line;                          // decided to switch to line
+  Double_t  fMaxRforHelix;                         // above this radius use straight line fit
   //
   const Int_t  *fkAxID;                            // axis IDs
   const Int_t  *fkAxCID;                           // axis combinations IDs