<commit_before>// $Id: AliHLTTPCGMTrackParam.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
// **************************************************************************
// This file is property of and copyright by the ALICE HLT Project          *
// ALICE Experiment at CERN, All rights reserved.                           *
//                                                                          *
// Primary Authors: <NAME> <<EMAIL>> *
//                  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     *
// without fee, provided that the above copyright notice appears in all     *
// copies and that both the copyright notice and this permission notice     *
// appear in the supporting documentation. The authors make no claims       *
// about the suitability of this software for any purpose. It is            *
// provided "as is" without express or implied warranty.                    *
//                                                                          *
//***************************************************************************


#include "AliHLTTPCGMTrackParam.h"
#include "AliHLTTPCCAMath.h"
#include "AliHLTTPCGMPhysicalTrackModel.h"
#include "AliHLTTPCGMPropagator.h"
#include "AliHLTTPCGMBorderTrack.h"
#include "AliHLTTPCGMMergedTrack.h"
#include "AliHLTTPCGMPolynomialField.h"
#ifndef HLTCA_STANDALONE
#include "AliExternalTrackParam.h"
#endif
#include "AliHLTTPCCAParam.h"
#include <cmath>
#include <stdlib.h>

#define DEBUG 0
#define PRINT_TRACKS 0
#define MIRROR 1
#define DOUBLE 1
#include "TFile.h"
#include "TNtuple.h"

GPUd() void AliHLTTPCGMTrackParam::Fit
(
 const AliHLTTPCGMPolynomialField &field,
 float x[], float y[], float z[], int row[], float alpha[], const AliHLTTPCCAParam &param,
 int &N, float &Alpha, bool UseMeanPt, float maxSinPhi
 ){

  const float kRho = 1.025e-3;//0.9e-3;
  const float kRadLen = 29.532;//28.94;
  
  static int nTracks = 0;
  if (PRINT_TRACKS || DEBUG) nTracks++;

  ResetCovariance();

  AliHLTTPCGMPropagator prop;
  prop.SetMaterial( kRadLen, kRho );
  prop.SetPolynomialField( field );  
  prop.SetUseMeanMomentum( UseMeanPt );
  prop.SetContinuousTracking( param.GetContinuousTracking() );
  prop.SetMaxSinPhi( maxSinPhi );
  
  int nWays = param.GetNWays();
  int maxN = N;
  int resetT0 = CAMath::Max(10.f, CAMath::Min(40.f, 150.f / fP[4]));
  for (int iWay = 0;iWay < nWays;iWay++)
  {
    if (DEBUG) printf("Fitting track %d way %d\n", nTracks, iWay);

    const bool rejectChi2ThisRound = ( nWays == 1 || iWay == 1 );
    const bool markNonFittedClusters = rejectChi2ThisRound && !(param.HighQPtForward() < fabs(fP[4]));
    const double kDeg2Rad = 3.14159265358979323846/180.;
    const float maxSinForUpdate = CAMath::Sin(70.*kDeg2Rad);
  
    prop.SetTrack( this, Alpha);

    N = 0;
    const bool inFlyDirection = iWay & 1;
    char directionState = 2;
    bool directionChangePending = 0;
    const int wayDirection = (iWay & 1) ? -1 : 1;
    int ihit = (iWay & 1) ? (maxN - 1) : 0;
    int lastRow = 0;
    for(;ihit >= 0 && ihit<maxN;ihit += wayDirection)
    {
      if (row[ihit] < 0) continue; // hit is excluded from fit
      const int rowType = row[ihit] < 64 ? 0 : row[ihit] < 128 ? 2 : 1;
      if (DEBUG) printf("\tHit %3d Row %3d: Cluster Alpha %8.3f    , X %8.3f - Y %8.3f, Z %8.3f\n", ihit, row[ihit], alpha[ihit], x[ihit], y[ihit], z[ihit]);
      
      char newState = N > 0 ? (lastRow > row[ihit]) : 2;
      bool doubleRow = ihit + wayDirection >= 0 && ihit + wayDirection < maxN && row[ihit] == row[ihit + wayDirection];
      float xx = x[ihit];
      float yy = y[ihit];
      float zz = z[ihit];
      if (DOUBLE && doubleRow)
      {
          float count = 1.;
          do
          {
              if (alpha[ihit] != alpha[ihit + wayDirection] || fabs(y[ihit] - y[ihit + wayDirection]) > 4. || fabs(z[ihit] - z[ihit + wayDirection]) > 4.) break;
              if (DEBUG) printf("    Merging hit row %d X %f Y %f Z %f\n", row[ihit], x[ihit], y[ihit], z[ihit]);
              ihit += wayDirection;
              xx += x[ihit];
              yy += y[ihit];
              zz += z[ihit];
              count += 1.;
          } while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && row[ihit] == row[ihit + wayDirection]);
          xx /= count;
          yy /= count;
          zz /= count;
          if (DEBUG) printf("    Double row (%d hits)\n", (int) count);
      }
      
      bool changeDirection = ((directionState ^ newState) == 1);
      if (changeDirection) directionChangePending = 1;
      bool canChangeDirection = (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && (newState ^ (row[ihit] > row[ihit + wayDirection])) == 1);
      directionState = newState;
      lastRow = row[ihit];
      if (DEBUG && (changeDirection || canChangeDirection || directionChangePending)) printf("    Change direction possibly change %d can %d pending %d\n", (int) changeDirection, (int) canChangeDirection, (int) directionChangePending);
      if (DEBUG) printf("\t%17sTrack   Alpha %8.3f %s, X %8.3f - Y %8.3f, Z %8.3f   -   QPt %7.2f (%7.2f), SinPhi %5.2f (%5.2f) %28s    ---   Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f   -   YPt %8.3f SPPt %8.3f YSP %8.3f\n", "", prop.GetAlpha(), (fabs(prop.GetAlpha() - alpha[ihit]) < 0.01 ? "   " : " R!"), fX, fP[0], fP[1], fP[4], prop.GetQPt0(), fP[2], prop.GetSinPhi0(), "", sqrt(fC[0]), sqrt(fC[2]), sqrt(fC[5]), sqrt(fC[14]), fC[10], fC[12], fC[3]);
      int err = prop.PropagateToXAlpha(xx, alpha[ihit], inFlyDirection );
      if (err == -2) //Rotation failed, try to bring to new x with old alpha first, rotate, and then propagate to x, alpha
      {
          if (prop.PropagateToXAlpha(xx, prop.GetAlpha(), inFlyDirection ) == 0)
            err = prop.PropagateToXAlpha(xx, alpha[ihit], inFlyDirection );
      }
      
      if (PRINT_TRACKS)
      {
          float ac = cos(alpha[ihit]), as = sin(alpha[ihit]);
          static int init = 0;
          static TFile *file;
          static TNtuple *nt;
          if (init == 0)
          {
            file = new TFile("tracksout.root","RECREATE");
            file->cd();
            nt = new TNtuple("field","field","track:cx:cy:cz:tx:ty:tz:mx:my:mz");      
            init = 1;
          }
          nt->Fill((float) nTracks, ac * xx - as * yy, as * xx + ac * yy, zz, ac * fX - as * fP[0], as * fX + ac * fP[0], fP[1], ac * prop.Model().X() - as * prop.Model().Y(), as * prop.Model().X() + ac * prop.Model().Y(), prop.Model().Z());
          if (nTracks == 29) {
            nt->Write();
            file->Write();
            file->Close();
            exit(0);
          }
      }
    
      if (DEBUG) printf("\t%17sPropaga Alpha %8.3f    , X %8.3f - Y %8.3f, Z %8.3f   -   QPt %7.2f (%7.2f), SinPhi %5.2f (%5.2f)   ---   Res %8.3f %8.3f   ---   Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f   -   YPt %8.3f SPPt %8.3f YSP %8.3f   -   Err %d\n", "", prop.GetAlpha(), fX, fP[0], fP[1], fP[4], prop.GetQPt0(), fP[2], prop.GetSinPhi0(), fP[0] - yy, fP[1] - zz, sqrt(fC[0]), sqrt(fC[2]), sqrt(fC[5]), sqrt(fC[14]), fC[10], fC[12], fC[3], err);
      const float Bz = prop.GetBz(prop.GetAlpha(), fX, fP[0], fP[1]);
      if (MIRROR && err == 0 && (changeDirection || canChangeDirection || directionChangePending) && fP[2] * fP[4] * Bz < 0)
      {
          float sinphi = fabs(fP[2]) > 0.99 ? 0.99 : fP[2];
          float cosphi = sqrt(1 - sinphi * sinphi);
          const float mirrordY = fP[0] - 2.0 * cosphi  / (fP[4] * Bz);
          if (DEBUG) printf(" -- MiroredY: %f --> %f", fP[0], mirrordY);
          if (changeDirection || ((directionChangePending || canChangeDirection) && fabs(yy - fP[0]) > fabs(yy - mirrordY)))
          {
              if (DEBUG) printf(" - Mirroring!!!");
              prop.Mirror(inFlyDirection);
              float err2Y, err2Z;
              prop.GetErr2(err2Y, err2Z, param, zz, rowType);
              prop.Model().Y() = fP[0] = yy;
              prop.Model().Z() = fP[1] = zz;
              if (fC[0] < err2Y) fC[0] = err2Y;
              if (fC[2] < err2Z) fC[2] = err2Z;
              if (fabs(fC[5]) < 0.1) fC[5] = fC[5] > 0 ? 0.1 : -0.1;
              if (fC[9] < 1.) fC[9] = 1.;
              prop.SetTrack(this, prop.GetAlpha());

              fNDF = -3;
              directionState = 2;
              directionChangePending = 0;
              N++;
              resetT0 = CAMath::Max(10.f, CAMath::Min(40.f, 150.f / fP[4]));
              if (DEBUG) printf("\n");<commit_msg>use new GetMirroredYtrack method instead of computing manually<commit_after>// $Id: AliHLTTPCGMTrackParam.cxx 41769 2010-06-16 13:58:00Z sgorbuno $
// **************************************************************************
// This file is property of and copyright by the ALICE HLT Project          *
// ALICE Experiment at CERN, All rights reserved.                           *
//                                                                          *
// Primary Authors: <NAME> <<EMAIL>> *
//                  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     *
// without fee, provided that the above copyright notice appears in all     *
// copies and that both the copyright notice and this permission notice     *
// appear in the supporting documentation. The authors make no claims       *
// about the suitability of this software for any purpose. It is            *
// provided "as is" without express or implied warranty.                    *
//                                                                          *
//***************************************************************************


#include "AliHLTTPCGMTrackParam.h"
#include "AliHLTTPCCAMath.h"
#include "AliHLTTPCGMPhysicalTrackModel.h"
#include "AliHLTTPCGMPropagator.h"
#include "AliHLTTPCGMBorderTrack.h"
#include "AliHLTTPCGMMergedTrack.h"
#include "AliHLTTPCGMPolynomialField.h"
#ifndef HLTCA_STANDALONE
#include "AliExternalTrackParam.h"
#endif
#include "AliHLTTPCCAParam.h"
#include <cmath>
#include <stdlib.h>

#define DEBUG 0
#define PRINT_TRACKS 0
#define MIRROR 1
#define DOUBLE 1
#include "TFile.h"
#include "TNtuple.h"

GPUd() void AliHLTTPCGMTrackParam::Fit
(
 const AliHLTTPCGMPolynomialField &field,
 float x[], float y[], float z[], int row[], float alpha[], const AliHLTTPCCAParam &param,
 int &N, float &Alpha, bool UseMeanPt, float maxSinPhi
 ){

  const float kRho = 1.025e-3;//0.9e-3;
  const float kRadLen = 29.532;//28.94;
  
  static int nTracks = 0;
  if (PRINT_TRACKS || DEBUG) nTracks++;

  ResetCovariance();

  AliHLTTPCGMPropagator prop;
  prop.SetMaterial( kRadLen, kRho );
  prop.SetPolynomialField( field );  
  prop.SetUseMeanMomentum( UseMeanPt );
  prop.SetContinuousTracking( param.GetContinuousTracking() );
  prop.SetMaxSinPhi( maxSinPhi );
  
  int nWays = param.GetNWays();
  int maxN = N;
  int resetT0 = CAMath::Max(10.f, CAMath::Min(40.f, 150.f / fP[4]));
  for (int iWay = 0;iWay < nWays;iWay++)
  {
    if (DEBUG) printf("Fitting track %d way %d\n", nTracks, iWay);

    const bool rejectChi2ThisRound = ( nWays == 1 || iWay == 1 );
    const bool markNonFittedClusters = rejectChi2ThisRound && !(param.HighQPtForward() < fabs(fP[4]));
    const double kDeg2Rad = 3.14159265358979323846/180.;
    const float maxSinForUpdate = CAMath::Sin(70.*kDeg2Rad);
  
    prop.SetTrack( this, Alpha);

    N = 0;
    const bool inFlyDirection = iWay & 1;
    char directionState = 2;
    bool directionChangePending = 0;
    const int wayDirection = (iWay & 1) ? -1 : 1;
    int ihit = (iWay & 1) ? (maxN - 1) : 0;
    int lastRow = 0;
    for(;ihit >= 0 && ihit<maxN;ihit += wayDirection)
    {
      if (row[ihit] < 0) continue; // hit is excluded from fit
      const int rowType = row[ihit] < 64 ? 0 : row[ihit] < 128 ? 2 : 1;
      if (DEBUG) printf("\tHit %3d Row %3d: Cluster Alpha %8.3f    , X %8.3f - Y %8.3f, Z %8.3f\n", ihit, row[ihit], alpha[ihit], x[ihit], y[ihit], z[ihit]);
      
      char newState = N > 0 ? (lastRow > row[ihit]) : 2;
      bool doubleRow = ihit + wayDirection >= 0 && ihit + wayDirection < maxN && row[ihit] == row[ihit + wayDirection];
      float xx = x[ihit];
      float yy = y[ihit];
      float zz = z[ihit];
      if (DOUBLE && doubleRow)
      {
          float count = 1.;
          do
          {
              if (alpha[ihit] != alpha[ihit + wayDirection] || fabs(y[ihit] - y[ihit + wayDirection]) > 4. || fabs(z[ihit] - z[ihit + wayDirection]) > 4.) break;
              if (DEBUG) printf("    Merging hit row %d X %f Y %f Z %f\n", row[ihit], x[ihit], y[ihit], z[ihit]);
              ihit += wayDirection;
              xx += x[ihit];
              yy += y[ihit];
              zz += z[ihit];
              count += 1.;
          } while (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && row[ihit] == row[ihit + wayDirection]);
          xx /= count;
          yy /= count;
          zz /= count;
          if (DEBUG) printf("    Double row (%d hits)\n", (int) count);
      }
      
      bool changeDirection = ((directionState ^ newState) == 1);
      if (changeDirection) directionChangePending = 1;
      bool canChangeDirection = (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && (newState ^ (row[ihit] > row[ihit + wayDirection])) == 1);
      directionState = newState;
      lastRow = row[ihit];
      if (DEBUG && (changeDirection || canChangeDirection || directionChangePending)) printf("    Change direction possibly change %d can %d pending %d\n", (int) changeDirection, (int) canChangeDirection, (int) directionChangePending);
      if (DEBUG) printf("\t%17sTrack   Alpha %8.3f %s, X %8.3f - Y %8.3f, Z %8.3f   -   QPt %7.2f (%7.2f), SinPhi %5.2f (%5.2f) %28s    ---   Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f   -   YPt %8.3f SPPt %8.3f YSP %8.3f\n", "", prop.GetAlpha(), (fabs(prop.GetAlpha() - alpha[ihit]) < 0.01 ? "   " : " R!"), fX, fP[0], fP[1], fP[4], prop.GetQPt0(), fP[2], prop.GetSinPhi0(), "", sqrt(fC[0]), sqrt(fC[2]), sqrt(fC[5]), sqrt(fC[14]), fC[10], fC[12], fC[3]);
      int err = prop.PropagateToXAlpha(xx, alpha[ihit], inFlyDirection );
      if (err == -2) //Rotation failed, try to bring to new x with old alpha first, rotate, and then propagate to x, alpha
      {
          if (prop.PropagateToXAlpha(xx, prop.GetAlpha(), inFlyDirection ) == 0)
            err = prop.PropagateToXAlpha(xx, alpha[ihit], inFlyDirection );
      }
      
      if (PRINT_TRACKS)
      {
          float ac = cos(alpha[ihit]), as = sin(alpha[ihit]);
          static int init = 0;
          static TFile *file;
          static TNtuple *nt;
          if (init == 0)
          {
            file = new TFile("tracksout.root","RECREATE");
            file->cd();
            nt = new TNtuple("field","field","track:cx:cy:cz:tx:ty:tz:mx:my:mz");      
            init = 1;
          }
          nt->Fill((float) nTracks, ac * xx - as * yy, as * xx + ac * yy, zz, ac * fX - as * fP[0], as * fX + ac * fP[0], fP[1], ac * prop.Model().X() - as * prop.Model().Y(), as * prop.Model().X() + ac * prop.Model().Y(), prop.Model().Z());
          if (nTracks == 29) {
            nt->Write();
            file->Write();
            file->Close();
            exit(0);
          }
      }
    
      if (DEBUG) printf("\t%17sPropaga Alpha %8.3f    , X %8.3f - Y %8.3f, Z %8.3f   -   QPt %7.2f (%7.2f), SinPhi %5.2f (%5.2f)   ---   Res %8.3f %8.3f   ---   Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f   -   YPt %8.3f SPPt %8.3f YSP %8.3f   -   Err %d\n", "", prop.GetAlpha(), fX, fP[0], fP[1], fP[4], prop.GetQPt0(), fP[2], prop.GetSinPhi0(), fP[0] - yy, fP[1] - zz, sqrt(fC[0]), sqrt(fC[2]), sqrt(fC[5]), sqrt(fC[14]), fC[10], fC[12], fC[3], err);
      const float Bz = prop.GetBz(prop.GetAlpha(), fX, fP[0], fP[1]);
      if (MIRROR && err == 0 && (changeDirection || canChangeDirection || directionChangePending) && fP[2] * fP[4] * Bz < 0)
      {
          const float mirrordY = prop.GetMirroredYTrack();
          if (DEBUG) printf(" -- MiroredY: %f --> %f", fP[0], mirrordY);
          if (changeDirection || ((directionChangePending || canChangeDirection) && fabs(yy - fP[0]) > fabs(yy - mirrordY)))
          {
              if (DEBUG) printf(" - Mirroring!!!");
              prop.Mirror(inFlyDirection);
              float err2Y, err2Z;
              prop.GetErr2(err2Y, err2Z, param, zz, rowType);
              prop.Model().Y() = fP[0] = yy;
              prop.Model().Z() = fP[1] = zz;
              if (fC[0] < err2Y) fC[0] = err2Y;
              if (fC[2] < err2Z) fC[2] = err2Z;
              if (fabs(fC[5]) < 0.1) fC[5] = fC[5] > 0 ? 0.1 : -0.1;
              if (fC[9] < 1.) fC[9] = 1.;
              prop.SetTrack(this, prop.GetAlpha());

              fNDF = -3;
              directionState = 2;
              directionChangePending = 0;
              N++;
              resetT0 = CAMath::Max(10.f, CAMath::Min(40.f, 150.f / fP[4]));
              if (DEBUG) printf("\n");
{'fraction_non_alphanumeric': 0.173925584078454, 'fraction_numerical': 0.03466974329391405, 'mean_word_length': 2.2079940784603997, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '125356', 'n_tokens_mistral': 6454, 'n_tokens_neox': 6075, 'n_words': 1583}