New classes for handling and testing new hit structures
authorkowal2 <kowal2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Nov 2000 07:40:20 +0000 (07:40 +0000)
committerkowal2 <kowal2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Nov 2000 07:40:20 +0000 (07:40 +0000)
TPC/AliTPCComparison.C
TPC/AliTPCTrackHits.cxx [new file with mode: 0644]
TPC/AliTPCTrackHits.h [new file with mode: 0644]
TPC/AliTPCTrackHitsInterfaces.h [new file with mode: 0644]
TPC/AliTPCtest.C
TPC/Makefile
TPC/TestTPCTrackHits.cxx [new file with mode: 0644]
TPC/TestTPCTrackHits.h [new file with mode: 0644]

index 2329373..74d540a 100644 (file)
@@ -1,3 +1,4 @@
+//#include "alles.h"
 struct GoodTrack {
   Int_t lab;
   Int_t code;
@@ -7,6 +8,8 @@ struct GoodTrack {
 Int_t good_tracks(GoodTrack *gt, Int_t max);
 
 Int_t AliTPCComparison() {
+  Int_t i;
+   gBenchmark->Start("AliTPCComparison");
    cerr<<"Doing comparison...\n";
 
    TFile *cf=TFile::Open("AliTPCclusters.root");
@@ -20,8 +23,9 @@ Int_t AliTPCComparison() {
    ca->SetClusterType("AliTPCcluster");
    ca->ConnectTree("Segment Tree");
    Int_t nentr=Int_t(ca->GetTree()->GetEntries());
-   for (Int_t i=0; i<nentr; i++) {
-       AliSegmentID *s=ca->LoadEntry(i);
+   for (i=0; i<nentr; i++) {
+     //AliSegmentID *s=;
+       ca->LoadEntry(i);
    }
 
 // Load tracks
@@ -30,7 +34,8 @@ Int_t AliTPCComparison() {
    TObjArray tarray(2000);
    TTree *tracktree=(TTree*)tf->Get("TreeT");
    TBranch *tbranch=tracktree->GetBranch("tracks");
-   Int_t nentr=tracktree->GetEntries();
+   nentr=(Int_t)tracktree->GetEntries();
+   
    for (i=0; i<nentr; i++) {
        AliTPCtrack *iotrack=new AliTPCtrack;
        tbranch->SetAddress(&iotrack);
@@ -93,29 +98,60 @@ Int_t AliTPCComparison() {
    TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
    TH2F *hep=new TH2F("hep","dE/dX vs momentum",50,0.,2.,50,0.,200.);
 
+   Int_t mingood = ngood;  //MI change
+   Int_t * track_notfound = new Int_t[ngood];
+   Int_t   itrack_notfound =0;
+   Int_t * track_fake  = new Int_t[ngood];
+   Int_t   itrack_fake = 0;
+   Int_t * track_multifound = new Int_t[ngood];
+   Int_t * track_multifound_n = new Int_t[ngood];
+   Int_t   itrack_multifound =0;
+
    while (ngood--) {
-      Int_t lab=gt[ngood].lab,tlab;
+      Int_t lab=gt[ngood].lab,tlab=-1;
       Float_t ptg=
       TMath::Sqrt(gt[ngood].px*gt[ngood].px + gt[ngood].py*gt[ngood].py);
 
       hgood->Fill(ptg);
 
-      AliTPCtrack *track;
-      Int_t i;
+      AliTPCtrack *track=0;
       for (i=0; i<nentr; i++) {
-          track=(AliTPCtrack*)tarray->UncheckedAt(i);
+          track=(AliTPCtrack*)tarray.UncheckedAt(i);
           tlab=track->GetLabel();
           if (lab==TMath::Abs(tlab)) break;
       }
       if (i==nentr) {
-       cerr<<"Track "<<lab<<" was not found !\n";
+       //      cerr<<"Track "<<lab<<" was not found !\n";
+       track_notfound[itrack_notfound++]=lab;
         continue;
       }
+      
+      //MI change  - addition
+      Int_t micount=0;
+      Int_t mi;
+      AliTPCtrack * mitrack;
+      for (mi=0; mi<nentr; mi++) {
+       mitrack=(AliTPCtrack*)tarray.UncheckedAt(mi);          
+       if (lab==TMath::Abs(mitrack->GetLabel())) micount++;
+      }
+      if (micount>1) {
+       //cout<<"Track no. "<<lab<<" found "<<micount<<"  times\n";
+       track_multifound[itrack_multifound]=lab;
+       track_multifound_n[itrack_multifound]=micount;
+       itrack_multifound++;
+      }
+      
+
+      //
       Double_t xk=gt[ngood].x;//digp->GetPadRowRadii(0,0);
       track->PropagateTo(xk);
 
       if (lab==tlab) hfound->Fill(ptg);
-      else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
+      else { 
+       track_fake[itrack_fake++]=lab;
+       hfake->Fill(ptg); 
+       //cerr<<lab<<" fake\n";
+      }
 
       Double_t px,py,pz,pt=fabs(track->GetPt());track->GetPxPyPz(px,py,pz);
 
@@ -143,11 +179,35 @@ Int_t AliTPCComparison() {
 
    }
 
+   
    Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<<ng<<endl;
    Stat_t nf=hfound->GetEntries();
    if (ng!=0) 
-      cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
+      cerr<<"\n\n\nIntegral efficiency is about "<<nf/ng*100.<<" %\n";
+   
+   //MI change  - addition
+   cout<<"Total number of found tracks ="<<nentr<<"\n";
+   cout<<"Total number of \"good\" tracks ="<<mingood<<"\n";
+
+
+   cout<<"\nList of Not found tracks :\n";
+   //   Int_t i;
+   for ( i = 0; i< itrack_notfound; i++){
+     cout<<track_notfound[i]<<"\t";
+     if ((i%5)==4) cout<<"\n";
+   }
+   cout<<"\nList of fake  tracks :\n";
+   for ( i = 0; i< itrack_fake; i++){
+     cout<<track_fake[i]<<"\t";
+     if ((i%5)==4) cout<<"\n";
+   }
+    cout<<"\nList of multiple found tracks :\n";
+   for ( i=0; i<itrack_multifound; i++) {
+       cout<<"id.   "<<track_multifound[i]<<"     found - "<<track_multifound_n[i]<<"times\n";
+   }
+   cout<<"\n\n\n";
 
+   
    gStyle->SetOptStat(111110);
    gStyle->SetOptFit(1);
 
@@ -211,6 +271,10 @@ Int_t AliTPCComparison() {
    hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)"); 
    hep->Draw(); c1->cd();
 
+   gBenchmark->Stop("AliTPCComparison");
+   gBenchmark->Show("AliTPCComparison");
+
+
    return 0;
 }
 
@@ -218,7 +282,7 @@ Int_t AliTPCComparison() {
 Int_t good_tracks(GoodTrack *gt, Int_t max) {
    Int_t nt=0;
 
-   TFile *file=TFile::Open("rfio:galice.root");
+   TFile *file=TFile::Open("galice.root");
    if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
 
    if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
@@ -247,6 +311,13 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
    Int_t *good=new Int_t[np];
    for (Int_t ii=0; ii<np; ii++) good[ii]=0;
 
+
+   //MI change to be possible compile macro
+   //definition out of the swith statemnet
+    Int_t sectors_by_rows=0;
+    TTree *TD=0;
+    AliSimDigits da, *digits=&da;
+    Int_t *count=0;
    switch (ver) {
    case 1:
      {
@@ -280,13 +351,12 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
      }
       break;
    case 2:
-      TTree *TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
-      AliSimDigits da, *digits=&da;
+      TD=(TTree*)gDirectory->Get("TreeD_75x40_100x60");
       TD->GetBranch("Segment")->SetAddress(&digits);
-      Int_t *count = new Int_t[np];
+      count = new Int_t[np];
       Int_t i;
       for (i=0; i<np; i++) count[i]=0;
-      Int_t sectors_by_rows=(Int_t)TD->GetEntries();
+      sectors_by_rows=(Int_t)TD->GetEntries();
       for (i=0; i<sectors_by_rows; i++) {
           if (!TD->GetEvent(i)) continue;
           Int_t sec,row;
@@ -306,9 +376,9 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
           for (Int_t j=0; j<np; j++) {
               if (count[j]>1) {
                  if (sec>=digp->GetNInnerSector())
-                 if (row==nrow_up-1    ) good[j]|=0x1000;
+                  if (row==nrow_up-1    ) good[j]|=0x1000;
                  if (sec>=digp->GetNInnerSector())
-                 if (row==nrow_up-1-gap) good[j]|=0x800;
+                  if (row==nrow_up-1-gap) good[j]|=0x800;
                  good[j]++;
               }
               count[j]=0;
@@ -323,13 +393,27 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
    }
 
    TTree *TH=gAlice->TreeH();
-   TClonesArray *hits=TPC->Hits();
-   Int_t npart=TH->GetEntries();
+   //   TClonesArray *hits=TPC->Hits();
+   Int_t npart=(Int_t)TH->GetEntries();
 
    while (npart--) {
+      AliTPChit *hit0=0;
+
       TPC->ResetHits();
       TH->GetEvent(npart);
-      Int_t nhits=hits->GetEntriesFast();
+      AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
+      while (hit){
+       if (hit->fQ==0.) break;
+       hit =  (AliTPChit*) TPC->NextHit();
+      }
+      if (hit) {
+       hit0 = new AliTPChit(*hit); //Make copy of hit
+       hit = hit0;
+      }
+      else continue;
+      AliTPChit *hit1=(AliTPChit*)TPC->NextHit();      
+      if (hit1==0) continue;
+      /*      Int_t nhits=hits->GetEntriesFast();
       if (nhits==0) continue;
       AliTPChit *hit;
       Int_t j;
@@ -339,6 +423,7 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
       }
       if (j==nhits-1) continue;
       AliTPChit *hit1=(AliTPChit*)hits->UncheckedAt(j+1);
+      */
       if (hit1->fQ != 0.) continue;
       Int_t i=hit->Track();
       TParticle *p = (TParticle*)particles->UncheckedAt(i);
@@ -355,8 +440,8 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
       gt[nt].x = hit1->X()*cs + hit1->Y()*sn;
       gt[nt].y =-hit1->X()*sn + hit1->Y()*cs;
       gt[nt].z = hit1->Z();
-      nt++;
-
+      nt++;    
+       if (hit0) delete hit0;
       cerr<<i<<"                \r";
       if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
    }
@@ -365,6 +450,8 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) {
    delete gAlice; gAlice=0;
 
    file->Close();
+   gBenchmark->Stop("AliTPCComparison");
+   gBenchmark->Show("AliTPCComparison");   
    return nt;
 }
 
diff --git a/TPC/AliTPCTrackHits.cxx b/TPC/AliTPCTrackHits.cxx
new file mode 100644 (file)
index 0000000..aadb4a5
--- /dev/null
@@ -0,0 +1,678 @@
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Time Projection Chamber  track hits object                                //
+//
+//  Origin: Marian Ivanov , GSI Darmstadt
+//
+//  Class for storing simulated AliTPCHits  for given track                  //
+//             -average compression comparing to classical ClonesArray is    //
+//              around 5-7 (depending on the required hit precision)       //
+//                                                                           //
+//Begin_Html
+/*
+<img src="gif/AliTPCTrackHits.gif">
+*/
+//End_Html
+//                                                                           //
+//                                                                          //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "TVector3.h"
+#include "TClonesArray.h"    
+#include "AliTPCTrackHits.h"
+#include "AliTPC.h"
+
+#include <iostream.h>
+
+
+
+ClassImp(AliTPCTrackHits) 
+LClassImp(AliTrackHitsInfo) 
+LClassImp(AliTrackHitsParam)  
+LClassImp(AliHitInfo)
+
+Int_t AliTrackHitsInfo::fgCounter1 =0;
+Int_t AliTrackHitsInfo::fgCounter2 =0;
+Int_t AliTrackHitsParam::fgCounter1 =0;
+Int_t AliTrackHitsParam::fgCounter2 =0;
+Int_t AliHitInfo::fgCounter1 =0;
+Int_t AliHitInfo::fgCounter2 =0;
+Int_t AliTPCTrackHits::fgCounter1 =0;
+Int_t AliTPCTrackHits::fgCounter2 =0;
+const Double_t AliTPCTrackHits::fgkPrecision=1e-6;  //precision 
+const Double_t AliTPCTrackHits::fgkPrecision2=1e-20;  //precision
+
+
+/************************************************************/
+//           Interface classes                              // 
+#include "AliTPCTrackHitsInterfaces.h"
+
+
+
+
+struct AliTPCCurrentHit {
+  AliTPChit fHit;
+  UInt_t   fInfoIndex;//   - current info pointer 
+  UInt_t   fParamIndex;//  - current param pointer
+  UInt_t   fStackIndex; // - current hit stack index
+  Double_t fR;   //current Radius
+  Bool_t  fStatus; //current status    
+};   
+
+
+struct  AliTPCTempHitInfo {
+  enum    { fkStackSize = 10000};
+  AliTPCTempHitInfo();   
+  void     NewParam(Double_t r, Double_t z, Double_t fi, Int_t q);
+  void     SetHit(Double_t r, Double_t z, Double_t fi, Int_t q);
+  Double_t * GetPosition(Int_t index){return &fPositionStack[index*3];}
+  void    UpdateParam(Double_t maxdelta); //recal
+  void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
+           Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
+           Double_t fSumX4, Int_t n,
+             Double_t &a, Double_t &b, Double_t &c);
+  void  Fit(AliTrackHitsParam * param);
+  Double_t fSumDr;    //
+  Double_t fSumDr2;   //
+  Double_t fSumDr3;   //  
+  Double_t fSumDr4;   //
+  Double_t fSumDFi;  //
+  Double_t fSumDFiDr; //  
+  Double_t fSumDFiDr2;//
+  Double_t fSumDZ;     //
+  Double_t fSumDZDr;  //
+  Double_t fSumDZDr2;  //
+  Double_t fOldR;     //previos r
+  Double_t fPositionStack[3*fkStackSize];  //position stack 
+  UInt_t   fQStack[fkStackSize];           //Q stack
+  UInt_t fStackIndex;   //current stack index 
+  UInt_t fInfoIndex;    //current track info index
+  UInt_t fParamIndex;   //current track parameters index
+  AliTrackHitsInfo  * fInfo; //current track info
+  AliTrackHitsParam * fParam; //current track param
+};
+
+
+AliTPCTempHitInfo::AliTPCTempHitInfo()
+{
+  //
+  //set to default value
+  fSumDr=fSumDr2=fSumDr3=fSumDr4=
+    fSumDFi=fSumDFiDr=fSumDFiDr2=
+    fSumDZ=fSumDZDr=fSumDZDr2=0;  
+  fStackIndex = 0;
+  fInfoIndex  = 0;
+  fParamIndex = 0;
+}
+
+
+void AliTPCTempHitInfo::NewParam(Double_t r, Double_t z, Double_t fi, Int_t q)
+{
+  //
+  //reset stack and sum parameters
+  //store line initial point
+  fSumDr=fSumDr2=fSumDr3=fSumDr4=
+    fSumDFi=fSumDFiDr=fSumDFiDr2=
+    fSumDZ=fSumDZDr=fSumDZDr2=0;  
+  fStackIndex=0;
+  fParam->fR = r;
+  fOldR = r;
+  fParam->fZ = z;
+  fParam->fFi = fi;
+  fParam->fAn = 0.;
+  fParam->fAd = 0.;
+  fParam->fTheta =0.;
+  fParam->fThetaD =0.;
+  SetHit(r,z,fi,q);
+}
+
+void AliTPCTempHitInfo::SetHit(Double_t r, Double_t z, Double_t fi, Int_t q)
+{
+  //
+  //add hit to the stack
+  //recalculate new estimete of line parameters
+  Double_t *f = GetPosition(fStackIndex);  
+  f[0] = r;
+  f[1] = z;
+  f[2] = fi;
+  fQStack[fStackIndex]=q;
+  if (fStackIndex==0) return;
+  Double_t dr  = (r-fParam->fR);
+  if (TMath::Abs(dr)<AliTPCTrackHits::fgkPrecision) dr =AliTPCTrackHits::fgkPrecision;
+  Double_t dfi = fi-fParam->fFi;
+  Double_t dz  = z -fParam->fZ; 
+  Double_t dr2 =dr*dr;
+  Double_t dr3 =dr2*dr;
+  Double_t dr4 =dr3*dr;
+  fSumDr +=dr;
+  fSumDr2+=dr2;
+  fSumDr3+=dr3;
+  fSumDr4+=dr4;
+  fSumDFi +=dfi;
+  fSumDFiDr+=dfi*dr;
+  fSumDFiDr2+=dfi*dr2;
+  fSumDZ +=dz;
+  fSumDZDr+=dz*dr;
+  fSumDZDr2+=dz*dr2;
+  
+  //update fit parameters
+  //
+  Double_t det = fSumDr2*fSumDr4-fSumDr3*fSumDr3;
+  if (TMath::Abs(det)<AliTPCTrackHits::fgkPrecision2) return;
+  if ( ( fStackIndex>1 )  ){
+    fParam->fAn = (fSumDr4*fSumDFiDr-fSumDr3*fSumDFiDr2)/det;
+    fParam->fAd = (fSumDr2*fSumDFiDr2-fSumDr3*fSumDFiDr)/det;
+  }
+  else
+    fParam->fAn = fSumDFiDr/fSumDr2;
+  if ( ( fStackIndex>1 )  ){
+    fParam->fTheta = (fSumDr4*fSumDZDr-fSumDr3*fSumDZDr2)/det;
+    fParam->fThetaD= (fSumDr2*fSumDZDr2-fSumDr3*fSumDZDr)/det;
+  }
+  else
+    fParam->fTheta = fSumDZDr/fSumDr2; 
+}
+
+
+void   AliTPCTempHitInfo::UpdateParam(Double_t maxdelta)
+{
+  //recalc parameters not fixing origin point
+  if (fStackIndex>5){ 
+    Double_t a,b,c;
+    a=b=c=0;
+    Fit2(fSumDFi, fSumDFiDr, fSumDFiDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
+        fStackIndex, a,b,c);
+    if (TMath::Abs(a)<maxdelta){
+      fParam->fFi +=a/fParam->fR;    
+      fParam->fAn = b;    
+      fParam->fAd = c;                  
+    }
+    Fit2(fSumDZ, fSumDZDr, fSumDZDr2, fSumDr,fSumDr2,fSumDr3,fSumDr4,
+        fStackIndex, a,b,c) ;   
+    if (TMath::Abs(a)<maxdelta){
+      fParam->fZ +=a;    
+      fParam->fTheta = b;    
+      fParam->fThetaD = c;   
+    }                         
+  }
+      
+}
+void   AliTPCTempHitInfo::Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
+           Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
+           Double_t fSumX4, Int_t n,
+           Double_t &a, Double_t &b, Double_t &c)
+{
+  //fit of second order
+  Double_t det = 
+    n* (fSumX2*fSumX4-fSumX3*fSumX3) -
+    fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
+    fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
+    
+  if (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) {    
+    a = 
+      (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
+       fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
+       fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
+    b=
+      (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
+      fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
+      fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
+    c=
+      (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
+       fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
+       fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;  
+  }
+}
+
+void   AliTPCTempHitInfo::Fit(AliTrackHitsParam * param)
+{
+  // fit fixing first and the last point 
+  //result stored in new param
+  Double_t dx2  = (GetPosition(fStackIndex))[0]-fParam->fR;
+  Double_t det = fSumDr4+dx2*fSumDr2-2*dx2*fSumDr3;
+  if ( (TMath::Abs(det)> AliTPCTrackHits::fgkPrecision) &&
+       ((TMath::Abs(dx2)> AliTPCTrackHits::fgkPrecision))){
+    Double_t dfi2 = (GetPosition(fStackIndex))[1]-fParam->fFi;
+    param->fAd = (fSumDFiDr2+dfi2*fSumDr-dx2*fSumDFiDr-dfi2*fSumDr3/dx2)/det;
+    param->fAn  = (dfi2-param->fAd*dx2*dx2)/dx2;
+    
+    Double_t dz2 = (GetPosition(fStackIndex))[1]-fParam->fZ;
+    param->fTheta = (fSumDZDr2+dz2*fSumDr-dx2*fSumDZDr-dz2*fSumDr3/dx2)/det;
+    param->fTheta  = (dz2-param->fAd*dx2*dx2)/dx2;
+  }
+  
+}
+
+
+
+
+AliTPCTrackHits::AliTPCTrackHits()
+{
+  //
+  //default constructor
+  //
+  const Float_t kHitPrecision=0.002; //default precision for hit position in cm
+  const Float_t kStep =0.003;  //30 mum step 
+  const UShort_t kMaxDistance =100;  //maximum distance 100  
+
+  fPrecision=kHitPrecision; //precision in cm
+  fStep = kStep; //step size
+  fMaxDistance = kMaxDistance; //maximum distance
+  fTempInfo =0;
+  fTrackHitsInfo = new AliObjectArray("AliTrackHitsInfo"); 
+  fTrackHitsParam = new AliObjectArray("AliTrackHitsParam");
+  fHitsPosAndQ = new TArrayOfArray_vStack("AliHitInfo");
+  fCurrentHit = new AliTPCCurrentHit;
+  fgCounter1++;
+  fgCounter2++;
+
+} 
+
+AliTPCTrackHits::~AliTPCTrackHits()
+{
+  //
+  //default destructor
+  //
+  if (fTrackHitsInfo) delete fTrackHitsInfo;
+  if (fTrackHitsParam) delete fTrackHitsParam;
+  if (fHitsPosAndQ) delete fHitsPosAndQ;
+  if (fCurrentHit) delete fCurrentHit;
+  if (fTempInfo) delete fTempInfo;
+  fgCounter1--;
+}
+
+void AliTPCTrackHits::Clear()
+{
+  //
+  //clear object 
+  //  fTrackHitsInfo->Clear();
+  //fTrackHitsParam->Clear();
+  
+  fTrackHitsInfo->Resize(0);
+  fTrackHitsParam->Resize(0);
+  fHitsPosAndQ->Clear();
+
+  if (fTempInfo){
+    delete fTempInfo; 
+    fTempInfo =0;
+  } 
+}
+
+
+void AliTPCTrackHits::AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x, 
+             Double_t y, Double_t z,Int_t q)
+{
+  Double_t r = TMath::Sqrt(x*x+y*y);
+  Double_t fi = TMath::ACos(x/r);
+  if (y<0) fi*=-1.;
+    AddHit(volumeID,trackID,r,z,fi,q);
+}
+
+void AliTPCTrackHits::AddHit(Int_t volumeID, Int_t trackID, 
+                            Double_t r, Double_t z, Double_t fi, Int_t q)
+{
+  //
+  Bool_t diff=kFALSE;
+  if (!fTempInfo) { //initialsation of track
+    fTempInfo = new AliTPCTempHitInfo;
+    //
+    if (fTrackHitsInfo->GetCapacity()<10) fTrackHitsInfo->Reserve(10);
+    fTrackHitsInfo->Resize(1);
+    fTempInfo->fInfoIndex =0;
+    if (fTrackHitsParam->GetCapacity()<100) fTrackHitsParam->Reserve(100);    
+    fTrackHitsParam->Resize(1);
+    //
+    fTempInfo->fInfo = 
+      (AliTrackHitsInfo*) (fTrackHitsInfo->At(0));
+    fTempInfo->fInfo->fVolumeID = volumeID;
+    fTempInfo->fInfo->fTrackID = trackID;
+    fTempInfo->fInfo->fHitParamIndex =0;     
+    fTempInfo->fInfoIndex = 0;
+    //
+    fTempInfo->fParam = 
+      (AliTrackHitsParam*) (fTrackHitsParam->At(0));
+    fTempInfo->fParamIndex = 0;
+    fTempInfo->NewParam(r,z,fi,q);
+    return;
+  }
+  
+  Int_t size = fHitsPosAndQ->GetSize();
+  if (size>(Int_t)fTempInfo->fParamIndex) {
+    fTempInfo->fParamIndex++;
+    if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
+      fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
+    fTempInfo->fParam = 
+      (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
+    fTempInfo->NewParam(r,z,fi,q);
+    return;
+  }
+  
+
+  // if new volume or new trackID  
+  if ( (volumeID!=fTempInfo->fInfo->fVolumeID) || 
+       (trackID!=fTempInfo->fInfo->fTrackID)){
+    diff=kTRUE;
+
+    FlushHitStack(kTRUE);
+      
+    fTempInfo->fInfoIndex++;
+    if (fTempInfo->fInfoIndex+1>fTrackHitsInfo->GetSize()) 
+      fTrackHitsInfo->Resize(fTempInfo->fInfoIndex+1);
+    fTempInfo->fInfo = 
+      (AliTrackHitsInfo*) (fTrackHitsInfo->At(fTempInfo->fInfoIndex));      
+    fTempInfo->fInfo->fVolumeID = volumeID;
+    fTempInfo->fInfo->fTrackID = trackID;
+    fTempInfo->fInfo->fHitParamIndex =fTempInfo->fParamIndex+1;  
+    // FlushHitStack(kTRUE);
+
+    fTempInfo->fParamIndex++;
+    if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
+      fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
+    fTempInfo->fParam = 
+      (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
+    fTempInfo->NewParam(r,z,fi,q);
+    return;
+  }
+     
+  //calculate current fit precission to next point
+  AliTrackHitsParam &param = *(fTempInfo->fParam);
+  Double_t dd=0;
+  Double_t dl=0;
+  Double_t ratio=0;
+  Double_t dr,dz,dfi,ddz,ddfi;
+  Double_t drhit,ddl;
+  dr=dz=dfi=ddz=ddfi=0;
+  drhit = r-fTempInfo->fOldR;
+  { 
+    //Double_t dfi2 = param.fAn+2*param.fAd*(r-param.fR); 
+    Double_t dfi2 = param.fAn;
+    dfi2*=dfi2*fTempInfo->fOldR*fTempInfo->fOldR;
+    //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(r-param.fR);
+    Double_t ddz2 =  param.fTheta;
+    ddz2*=ddz2;
+    ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
+  }
+  dl = fStep * Short_t(TMath::Nint(drhit*ratio/fStep));
+  ddl = dl - drhit*ratio; 
+  fTempInfo->fOldR += dl/ratio; 
+
+  if (fTempInfo->fStackIndex>2){     
+    dr = r-param.fR;        
+    dz =  z-param.fZ;  
+    dfi = fi-param.fFi;
+    ddz = dr*param.fTheta+dr*dr*param.fThetaD-dz;
+    ddfi= dr*param.fAn+dr*dr*param.fAd-dfi;    
+    dd  = TMath::Sqrt(ddz*ddz+r*r*ddfi*ddfi+ddl*ddl); 
+    //
+  }        
+  //safety factor 1.25
+  if ( ( (dd*1.25>fPrecision) ) ||  
+       (fTempInfo->fStackIndex+4>fTempInfo->fkStackSize) || 
+       (TMath::Abs(dl/fStep)>fMaxDistance)  ) 
+    diff=kTRUE;
+  else{
+    fTempInfo->fStackIndex++;   
+    fTempInfo->SetHit(r,z,fi,q);
+    return;
+  }  
+  //if parameter changed 
+  if (FlushHitStack(kFALSE)){   //if full buffer flushed
+    fTempInfo->fParamIndex++;
+    if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
+      fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
+    fTempInfo->fParam = 
+      (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));  
+    fTempInfo->NewParam(r,z,fi,q);
+  }
+  else{
+    fTempInfo->fStackIndex++;
+    fTempInfo->SetHit(r,z,fi,q);              
+  }
+}   
+
+Bool_t AliTPCTrackHits::FlushHitStack(Bool_t force)
+{
+  //
+  //write fHitsPosAndQ information from the stack to te arrays
+  if (!fTempInfo) return kFALSE; 
+  Int_t size = fHitsPosAndQ->GetSize();
+
+  if ( (size>0)&&(size!=(Int_t)fTempInfo->fParamIndex)) return kFALSE;  
+
+  if (fHitsPosAndQ->Push(fTempInfo->fStackIndex+1)!=fTempInfo->fParamIndex){
+    cout<<"internal error - contact MI\n";
+    return kFALSE;
+  }
+  AliHitInfo * info;
+  AliTrackHitsParam & param = *(fTempInfo->fParam);
+  //recalculate track parameter not fixing first point
+  fTempInfo->UpdateParam(fStep/4.);
+  //fTempInfo->Fit(fTempInfo->fParam);  //- fixing the first and the last point
+
+  Double_t oldr = param.fR; 
+  //cout<<"C3"<<fTempInfo->fStackIndex<<"\n"<<flush;
+  UInt_t i;
+  Double_t dd;
+  for (i=0; i <= fTempInfo->fStackIndex; i++){
+    Double_t * position = fTempInfo->GetPosition(i);
+    Double_t   dr = position[0]-oldr;
+    Double_t   ratio; 
+    { 
+      //Double_t dfi2 = param.fAn+2*param.fAd*(position[0]-param.fR);
+      Double_t dfi2 = param.fAn;
+      dfi2*=dfi2*oldr*oldr;
+      //Double_t ddz2 =  param.fTheta+2*param.fThetaD*(position[0]-param.fR);
+      Double_t ddz2 =  param.fTheta;
+      ddz2*=ddz2;
+      ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
+    }
+
+    Double_t dl = fStep*(Short_t)TMath::Nint(dr*ratio/fStep);  
+    dr = dl/ratio; 
+    oldr+=dr;
+    //calculate precission
+    AliTrackHitsParam &param = *(fTempInfo->fParam);    
+    //real deltas
+    Double_t dr1=  position[0]-param.fR;
+    Double_t dz =  position[1]-param.fZ;
+
+    Double_t dfi = position[2]-param.fFi;
+    //extrapolated deltas
+    Double_t dr2 = oldr-param.fR; 
+    Double_t ddr = dr2-dr1;
+    Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
+    Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
+    dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
+
+    if ( (dd>fPrecision) ){ 
+      if (i==0){
+       param.fAn = 0;
+       param.fAd = 0;
+       param.fTheta =0;
+        param.fThetaD =0;
+       Double_t ddz = dr2*param.fTheta+dr2*dr2*param.fThetaD-dz;
+       Double_t ddfi= dr2*param.fAn+dr2*dr2*param.fAd-dfi;    
+       dl = 0;
+       dd = TMath::Sqrt(ddz*ddz+oldr*oldr*ddfi*ddfi+ddr*ddr); 
+      }
+      else
+       break;
+    }
+
+    info = (AliHitInfo*)(fHitsPosAndQ->At(fTempInfo->fParamIndex,i));
+    info->fHitDistance = Short_t(TMath::Nint(dl/fStep));
+    info->fCharge = Short_t(fTempInfo->fQStack[i]);
+    /*
+    cout<<"C2";
+    cout<<" "<<fTempInfo->fStackIndex<<" \t";
+    cout<<" "<<i<<" \t";
+    cout<<position[0]<<"\t";
+    cout<<position[1]<<"\t"; 
+    cout<<position[2]<<"\t";
+    cout<<param.fAn<<"\t";
+    cout<<param.fTheta<<"\t";
+    cout<<dr1<<"\t"<<ddr<<"\t"<<ddz<<"\t"<<ddfi<<"\t"<<dd<<"\n"<<flush;    
+    */
+  }    
+  
+  if (i<=fTempInfo->fStackIndex){ //if previous iteration not succesfull 
+    fHitsPosAndQ->Resize(fTempInfo->fParamIndex,i);
+    //
+    fTempInfo->fParamIndex++;
+    if (fTempInfo->fParamIndex+1>fTrackHitsParam->GetSize()) 
+      fTrackHitsParam->Resize(fTempInfo->fParamIndex+1);    
+    fTempInfo->fParam = 
+       (AliTrackHitsParam*) (fTrackHitsParam->At(fTempInfo->fParamIndex));    
+    Double_t * p = fTempInfo->GetPosition(i);
+    UInt_t index2 = fTempInfo->fStackIndex;
+    fTempInfo->NewParam(p[0],p[1],p[2],fTempInfo->fQStack[i]);        
+    if (i+1<=index2) FlushHitStack2(i+1,index2);
+
+    if (force) return      FlushHitStack(kTRUE);      
+    return kFALSE;
+  }  
+  return kTRUE;
+} 
+
+void AliTPCTrackHits::FlushHitStack2(Int_t index1, Int_t index2)
+{
+  //
+  // second iteration flush stack
+  // call only for hits where first iteration were not succesfully interpolated  
+  Double_t * positionstack = new Double_t[3*(index2-index1+1)];
+  UInt_t   * qstack        = new UInt_t[index2-index1+1];
+  memcpy(positionstack, &fTempInfo->fPositionStack[3*index1],
+        (3*(index2-index1+1))*sizeof(Double_t));
+  memcpy(qstack, &fTempInfo->fQStack[index1],(index2-index1+1)*sizeof(UInt_t));
+  Double_t *p = positionstack;
+  for (Int_t j=0; j<=index2-index1;j++){ 
+    fTempInfo->fStackIndex++;
+    fTempInfo->SetHit(p[3*j+0],p[3*j+1],p[3*j+2],qstack[j]);
+  }  
+  delete []positionstack;
+  delete []qstack;
+}
+
+
+   
+
+
+
+
+  
+
+Bool_t AliTPCTrackHits::First()
+{
+  //
+  //set Current hit for the first hit
+  //
+  AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(0);
+  AliTrackHitsParam *param = (AliTrackHitsParam *)fTrackHitsParam->At(0);
+  AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(0,0);
+
+  if (!(info) || !(param) || !(hinfo) ) {
+    fCurrentHit->fStatus = kFALSE;
+    return kFALSE;
+  }
+
+  fCurrentHit->fInfoIndex  = 0;
+  fCurrentHit->fParamIndex = 0;
+  fCurrentHit->fStackIndex = 0;
+
+  fCurrentHit->fHit.fSector = info->fVolumeID;
+  fCurrentHit->fHit.SetTrack(info->fTrackID);
+  fCurrentHit->fHit.SetX(param->fR*TMath::Cos(param->fFi));
+  fCurrentHit->fHit.SetY(param->fR*TMath::Sin(param->fFi));
+  fCurrentHit->fHit.SetZ(param->fZ); 
+  fCurrentHit->fHit.fQ = hinfo->fCharge;
+   
+  fCurrentHit->fR = param->fR;
+  
+  return fCurrentHit->fStatus = kTRUE;
+}
+
+Bool_t AliTPCTrackHits::Next()
+{
+  //
+  //  
+  if (!(fCurrentHit->fStatus)) 
+    return kFALSE;
+
+  fCurrentHit->fStackIndex++;
+  AliHitInfo * hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,
+                                                     fCurrentHit->fStackIndex);
+  AliTrackHitsInfo *info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
+  AliTrackHitsParam *param =  (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
+
+  if (!hinfo) {
+    hinfo = (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex+1, 0);
+    if (!hinfo) 
+      return fCurrentHit->fStatus = kFALSE;
+    if (hinfo){ 
+      fCurrentHit->fParamIndex++;
+      fCurrentHit->fStackIndex = 0;
+      param = (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex);
+      if (!param) 
+       return fCurrentHit->fStatus = kFALSE;     
+      fCurrentHit->fR = param->fR;
+
+      if ((fCurrentHit->fInfoIndex+1<fTrackHitsInfo->GetSize())
+       &&((info+1)->fHitParamIndex<=fCurrentHit->fParamIndex)){
+       fCurrentHit->fInfoIndex++;
+       info = (AliTrackHitsInfo *)fTrackHitsInfo->At(fCurrentHit->fInfoIndex);
+       if (!info) 
+         return fCurrentHit->fStatus = kFALSE;
+       fCurrentHit->fHit.fSector = info->fVolumeID;
+       fCurrentHit->fHit.SetTrack(info->fTrackID);
+      }
+    }  
+  } 
+  Double_t ratio;
+  { 
+    //    Double_t dfi2 = param->fAn+2*param->fAd*(fCurrentHit->fR-param->fR);
+    Double_t dfi2 = param->fAn;
+    dfi2*=dfi2*fCurrentHit->fR*fCurrentHit->fR;
+    //    Double_t ddz2 = param->fTheta+2*param->fThetaD*(fCurrentHit->fR-param->fR);
+    Double_t ddz2 =  param->fTheta;
+    ddz2*=ddz2;
+    ratio = TMath::Sqrt(1.+ dfi2+ ddz2);  
+  }
+
+  fCurrentHit->fHit.fQ = hinfo->fCharge;
+  fCurrentHit->fR += fStep*hinfo->fHitDistance/ratio;
+  Double_t dR = fCurrentHit->fR - param->fR;
+  //Double_t dR =0;
+  Double_t fi = param->fFi + (param->fAn*dR+param->fAd*dR*dR);
+  Double_t z  = param->fZ + (param->fTheta*dR+param->fThetaD*dR*dR);
+  
+  fCurrentHit->fHit.SetX(fCurrentHit->fR*TMath::Cos(fi));
+  fCurrentHit->fHit.SetY(fCurrentHit->fR*TMath::Sin(fi));
+  fCurrentHit->fHit.SetZ(z);  
+  return kTRUE;
+}
+  
+AliTPChit * AliTPCTrackHits::GetHit()
+{
+  //
+   return (fCurrentHit->fStatus)? &fCurrentHit->fHit:0;
+  //return &fCurrentHit->fHit;
+
+}  
+
+AliTrackHitsParam * AliTPCTrackHits::GetParam()
+{
+  //
+  return (fCurrentHit->fStatus)? (AliTrackHitsParam *)fTrackHitsParam->At(fCurrentHit->fParamIndex) :0;
+} 
+
+AliHitInfo * AliTPCTrackHits::GetHitInfo()
+{
+  //
+  return (fCurrentHit->fStatus)? 
+    (AliHitInfo *)fHitsPosAndQ->At(fCurrentHit->fParamIndex,fCurrentHit->fStackIndex) :0;
+} 
+
+
+
diff --git a/TPC/AliTPCTrackHits.h b/TPC/AliTPCTrackHits.h
new file mode 100644 (file)
index 0000000..edd0c0b
--- /dev/null
@@ -0,0 +1,104 @@
+#ifndef ALITPCTRACKHITS_H
+#define ALITPCTRACKHITS_H
+////////////////////////////////////////////////
+//  Manager class for TPC   clusters                   //
+////////////////////////////////////////////////
+
+#include "AliCTypes.h"
+#include "AliSegmentID.h"
+#include "AliArrayS.h"
+#include "AliTPC.h"
+#include "TVector3.h"
+#include "AliObjectArray.h"
+#include "TArrayOfArray.h"
+
+class TClonesArray;
+class AliArrayS;
+class AliTPChit;
+class AliTPCTempHitInfo;
+class AliTPCCurrentHit;
+
+
+class AliTrackHitsInfo {
+public:
+  AliTrackHitsInfo(){fgCounter1++;fgCounter2++;}
+  ~AliTrackHitsInfo(){fgCounter1--;}
+  Int_t   fTrackID;  //track ID
+  Int_t   fVolumeID;   //volume ID
+  UInt_t   fHitParamIndex; //corresponding index  
+  static Int_t fgCounter1;
+  static Int_t fgCounter2;  
+  LClassDef(AliTrackHitsInfo,1)  
+};
+
+
+class AliTrackHitsParam {
+public:
+  AliTrackHitsParam(){fgCounter1++;fgCounter2++;}
+  ~AliTrackHitsParam(){fgCounter1--;}
+  Float_t fR;  //radius
+  Float_t fZ;  //z position
+  Float_t fFi; //radial angle
+  Float_t fAn; //angle with  the radial vector
+  Float_t fAd; //derivation of angle
+  Float_t fTheta; //theta angle
+  Float_t fThetaD; //theta angle derivation
+  static Int_t fgCounter1;
+  static Int_t fgCounter2;  
+  LClassDef(AliTrackHitsParam,1)  
+};
+
+
+class AliHitInfo {
+public:
+  AliHitInfo(){fgCounter1++;fgCounter2++;}
+  ~AliHitInfo(){fgCounter1--;}
+  Short_t fHitDistance; //distance to previous hit
+  Short_t fCharge; //deponed charge
+  static Int_t fgCounter1;
+  static Int_t fgCounter2;  
+  LClassDef(AliHitInfo,1)
+};
+
+
+
+class AliTPCTrackHits : public TObject{
+public:
+  AliTPCTrackHits(); 
+  ~AliTPCTrackHits();
+  void Clear();
+  void AddHitKartez(Int_t volumeID, Int_t trackID, Double_t x, 
+                   Double_t y, Double_t z,Int_t q);
+  void AddHit(Int_t volumeID, Int_t trackID, Double_t r, 
+             Double_t z, Double_t fi,Int_t q);
+  Bool_t First(); //set current hit to first hit 
+  Bool_t Next();  //set current hit to next
+  AliTPChit * GetHit();
+  AliTrackHitsParam * GetParam();
+  AliHitInfo * GetHitInfo();
+  Int_t  GetEntriesFast() { return fHitsPosAndQ ? fHitsPosAndQ->ArraySize():0;}
+  void SetHitPrecision(Double_t prec) {fPrecision=prec;}
+  void SetStepPrecision(Double_t prec) {fStep=prec;}
+  void SetMaxDistance(UInt_t distance) {fMaxDistance = distance;}
+  Bool_t  FlushHitStack(Bool_t force=kTRUE);    //
+public:
+  void FlushHitStack2(Int_t index1, Int_t index2);   //
+  AliObjectArray * fTrackHitsInfo;  //quick information about track
+  AliObjectArray * fTrackHitsParam;  //hit information  
+  TArrayOfArray_vStack * fHitsPosAndQ;  //position information
+
+  Double_t fPrecision;  // required precision
+  Double_t fStep;       //unit step size
+  UInt_t fMaxDistance;   //maximal distance between two connected hits 
+  AliTPCTempHitInfo * fTempInfo; //!information about track
+  AliTPCCurrentHit  * fCurrentHit; //!information about current hit 
+  static const Double_t fgkPrecision;  //precision 
+  static const Double_t fgkPrecision2;  //precision
+  static Int_t fgCounter1;
+  static Int_t fgCounter2;  
+  ClassDef(AliTPCTrackHits,1) 
+};
+
+
+#endif //ALITPCTRACKHITS_H
diff --git a/TPC/AliTPCTrackHitsInterfaces.h b/TPC/AliTPCTrackHitsInterfaces.h
new file mode 100644 (file)
index 0000000..0623672
--- /dev/null
@@ -0,0 +1,110 @@
+#include "AliClassInfo.h"
+#include "AliTPCTrackHits.h"
+
+/************************************************/
+/* Automaticaly generated interface for class     
+                 AliTrackHitsInfo                                
+**************************************************/
+
+
+class AliClassAliTrackHitsInfo : public AliClassInfo {
+public:
+       AliClassAliTrackHitsInfo(){
+         SetName("AliTrackHitsInfo");
+         SetTitle("Interface for AliTrackHitsInfo class ");
+         fgList.Add(this);
+         fSize = sizeof(AliTrackHitsInfo);
+       }
+       const char * GetClassName(){ return "AliTrackHitsInfo";}
+       virtual TClass* GetClass(){return AliTrackHitsInfo::Class();}
+       void CTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliTrackHitsInfo * last = &((AliTrackHitsInfo*)pointer)[size];
+         AliTrackHitsInfo * p = (AliTrackHitsInfo*)pointer;
+         while (p!=last) new (p++)AliTrackHitsInfo;
+       }
+       void DTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliTrackHitsInfo * last = &((AliTrackHitsInfo*)pointer)[size];
+         AliTrackHitsInfo * p = (AliTrackHitsInfo*)pointer;
+         while (p!=last) (p++)->~AliTrackHitsInfo();
+       }
+       void StreamBuffer(TBuffer &b,const void * pointer, UInt_t size=1)
+       {
+         for (UInt_t i=0;i<size;i++) ((AliTrackHitsInfo*)pointer)[i].Streamer(b);
+       }
+         void ObjectDump(void *p) {((AliTrackHitsInfo*)p)->Dump();}
+};
+AliClassAliTrackHitsInfo galiclass____AliTrackHitsInfo; 
+
+/************************************************/
+/* Automaticaly generated interface for class     
+                 AliTrackHitsParam                                
+**************************************************/
+
+
+class AliClassAliTrackHitsParam : public AliClassInfo {
+public:
+       AliClassAliTrackHitsParam(){
+         SetName("AliTrackHitsParam");
+         SetTitle("Interface for AliTrackHitsParam class ");
+         fgList.Add(this);
+         fSize = sizeof(AliTrackHitsParam);
+       }
+       const char * GetClassName(){ return "AliTrackHitsParam";}
+       virtual TClass* GetClass(){return AliTrackHitsParam::Class();}
+       void CTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliTrackHitsParam * last = &((AliTrackHitsParam*)pointer)[size];
+         AliTrackHitsParam * p = (AliTrackHitsParam*)pointer;
+         while (p!=last) new (p++)AliTrackHitsParam;
+       }
+       void DTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliTrackHitsParam * last = &((AliTrackHitsParam*)pointer)[size];
+         AliTrackHitsParam * p = (AliTrackHitsParam*)pointer;
+         while (p!=last) (p++)->~AliTrackHitsParam();
+       }
+       void StreamBuffer(TBuffer &b,const void * pointer, UInt_t size=1)
+       {
+         for (UInt_t i=0;i<size;i++) ((AliTrackHitsParam*)pointer)[i].Streamer(b);
+       }
+         void ObjectDump(void *p) {((AliTrackHitsParam*)p)->Dump();}
+};
+AliClassAliTrackHitsParam galiclass____AliTrackHitsParam; 
+
+/************************************************/
+/* Automaticaly generated interface for class     
+                 AliHitInfo                                
+**************************************************/
+
+
+class AliClassAliHitInfo : public AliClassInfo {
+public:
+       AliClassAliHitInfo(){
+         SetName("AliHitInfo");
+         SetTitle("Interface for AliHitInfo class ");
+         fgList.Add(this);
+         fSize = sizeof(AliHitInfo);
+       }
+       const char * GetClassName(){ return "AliHitInfo";}
+       virtual TClass* GetClass(){return AliHitInfo::Class();}
+       void CTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliHitInfo * last = &((AliHitInfo*)pointer)[size];
+         AliHitInfo * p = (AliHitInfo*)pointer;
+         while (p!=last) new (p++)AliHitInfo;
+       }
+       void DTORBuffer(void * pointer, UInt_t size=1)
+       {
+         AliHitInfo * last = &((AliHitInfo*)pointer)[size];
+         AliHitInfo * p = (AliHitInfo*)pointer;
+         while (p!=last) (p++)->~AliHitInfo();
+       }
+       void StreamBuffer(TBuffer &b,const void * pointer, UInt_t size=1)
+       {
+         for (UInt_t i=0;i<size;i++) ((AliHitInfo*)pointer)[i].Streamer(b);
+       }
+         void ObjectDump(void *p) {((AliHitInfo*)p)->Dump();}
+};
+AliClassAliHitInfo galiclass____AliHitInfo; 
index f6476b2..7f6d197 100644 (file)
@@ -17,8 +17,8 @@ Int_t AliTPCtest() {
      gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCHits2Digits.C");
      if (rc=AliTPCHits2Digits()) return rc;
 
-     gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCDisplayDigits.C");
-     if (rc=AliTPCDisplayDigits(1,1)) return rc;
+     //     gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCDisplayDigits.C");
+     //     if (rc=AliTPCDisplayDigits(1,1)) return rc;
    }
 
 
@@ -26,8 +26,8 @@ Int_t AliTPCtest() {
    gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCFindClusters.C");
    if (rc=AliTPCFindClusters()) return rc;
 
-   gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCDisplayClusters.C");
-   if (rc=AliTPCDisplayClusters()) return rc;
+   //  gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCDisplayClusters.C");
+   // if (rc=AliTPCDisplayClusters()) return rc;
 
    gROOT->LoadMacro("$(ALICE_ROOT)/TPC/AliTPCFindTracks.C");
    if (rc=AliTPCFindTracks()) return rc;
index 00ef48f..4318578 100644 (file)
@@ -9,17 +9,17 @@ PACKAGE = TPC
 
 # C++ sources
 
-SRCS          = AliDigits.cxx   AliSegmentArray.cxx AliTPCClusterFinder.cxx  AliTPC.cxx  AliTPCv0.cxx    AliTPCv1.cxx      AliTPCv2.cxx  \
+SRCS          = AliTPCClusterFinder.cxx  AliTPC.cxx  AliTPCv0.cxx    AliTPCv1.cxx      AliTPCv2.cxx  \
                AliTPCv3.cxx AliDetectorParam.cxx AliTPCParam.cxx \
                AliTPCParamSR.cxx AliTPCParamCR.cxx  AliTPCPRF2D.cxx \
-               AliTPCRF1D.cxx AliArrayI.cxx AliArrayS.cxx \
-               AliSegmentID.cxx  \
+               AliTPCRF1D.cxx \
                AliSimDigits.cxx AliDigitsArray.cxx  AliTPCDigitsArray.cxx \
                AliComplexCluster.cxx AliClusters.cxx AliClustersArray.cxx \
                 AliTPCClustersRow.cxx \
                AliTPCClustersArray.cxx AliH2F.cxx \
                 AliTPCcluster.cxx AliTPCclusterer.cxx \
-                AliTPCtrack.cxx AliTPCtracker.cxx
+                AliTPCtrack.cxx AliTPCtracker.cxx \
+               AliTPCTrackHits.cxx
 
 
 # C++ Headers
@@ -46,8 +46,8 @@ OBJS          = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
 
 # C++ compilation flags
 
-CXXFLAGS      = $(CXXOPTS) -g -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include
-#CXXFLAGS      = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/
+CXXFLAGS      = $(CXXOPTS) -g -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include \
+               -I$(ALICE_ROOT)/CONTAINERS
  
 # FORTRAN compilation flags
 
@@ -58,10 +58,12 @@ FFLAGS      = $(FOPT)
 # Target
 
 SLIBRARY       = $(LIBDIR)/libTPC.$(SL)
+ALIBRARY      = $(LIBDIR)/static/libTPC.a
 
-default:       $(SLIBRARY)
+default:       $(SLIBRARY) 
 
 $(LIBDIR)/libTPC.$(SL):  $(OBJS)
+$(LIBDIR)/static/libTPC.a:  $(OBJS)
 
 $(DICT):               $(HDRS)
 
diff --git a/TPC/TestTPCTrackHits.cxx b/TPC/TestTPCTrackHits.cxx
new file mode 100644 (file)
index 0000000..846d108
--- /dev/null
@@ -0,0 +1,335 @@
+/*
+  Author : MI
+  macro to compare TClonesArray hits with interpolated hits
+  ConvertHits1 read 
+*/
+
+#include "alles.h"
+#include "AliObjectArray.h"
+#include "AliTPCTrackHits.h"
+#include "AliArrayBranch.h"
+#include "TestTPCTrackHits.h"
+ClassImp(AliTPChitD)
+void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits,  Bool_t debug, TClonesArray *arrd=0);
+AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits);
+
+void ConvertHits(const char * benchmark, Bool_t debug)
+{      
+  //
+  //convert - not parametrised hits stored in Clones array
+  //to
+
+
+  TFile f("galice.root","update");  
+  TClonesArray *arr = new TClonesArray("AliTPChit",100);
+  TTree * treeCl =(TTree*)f.Get("TreeH0");
+  TBranch *branch = treeCl->GetBranch("TPC"); 
+  AliTPCTrackHits *  myhits = new AliTPCTrackHits;
+  branch->SetAddress(&arr);
+  //particle  
+  TTree * treeP = (TTree*)f.Get("TreeK0");
+  TBranch *branchP = treeP->GetBranch("Particles"); 
+  TClonesArray *arrp = new TClonesArray("TParticle",100);
+  branchP->SetAddress(&arrp);
+  branchP->GetEvent(0);
+  //
+  TFile f2("treeh.root","recreate");
+  f2.SetCompressionLevel(10);
+  TTree * treeh2 = new TTree("TreeTPCH","TreeTPCH");
+  //treeh2->AliBranch("TPC","AliTPCTrackHits", &myhits,4096);
+  treeh2->GetListOfBranches()->Add(new AliObjectBranch("TPC","AliTPCTrackHits", 
+        &myhits,treeh2,4096,1));
+
+  TFile f3("treehcl.root","recreate");
+  f3.SetCompressionLevel(2);
+  TTree * treeh3 = new TTree("TreeTPCH","TreeTPCH");
+  treeh3->Branch("TPC", &arr,64000,2);
+
+  gBenchmark->Start(benchmark); 
+  Int_t trsize = (Int_t)treeCl->GetEntries();
+  for (Int_t i=0;i<300;i++){    
+    arr->Clear(); 
+    Int_t size = branch->GetEvent(i);     
+    //if (size>0){
+    if ((size>0)&& arr->GetEntriesFast()>0) {    
+      f.cd();
+      myhits->Clear();
+      myhits =MakeTrack(arr,arrp,myhits);            
+      CompareHits(arr,myhits,debug);
+      f2.cd();
+      treeh2->Fill();
+      f3.cd();
+      treeh3->Fill();      
+      if ((i%10)==0) cout<<i<<"\n"; 
+    }    
+  }
+  f3.cd();
+  treeh3->Write();
+  f2.cd();
+  treeh2->Write();
+  f2.Close();
+  f.Close();
+  delete myhits;
+  gBenchmark->Stop(benchmark);
+  gBenchmark->Show(benchmark); 
+}  
+
+
+
+void CompareHits(const char * benchmark, Bool_t debug)
+{        
+  //compare persistent hits
+  TFile f2("treeh.root");
+  TTree * treeh2 = (TTree*)f2.Get("TreeTPCH");
+  AliTPCTrackHits *  myhits = new AliTPCTrackHits ; 
+  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC"); 
+  mybranch->SetAddress(&myhits);
+
+    
+  TClonesArray *arr = new TClonesArray("AliTPChit",100);
+  TFile f3("treehcl.root");
+  TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
+  TBranch *branch = treeh3->GetBranch("TPC"); 
+  branch->SetAddress(&arr);
+   
+  cout<<"Lets go!\n"; 
+  gBenchmark->Start(benchmark);
+  Int_t trsize = (Int_t)treeh3->GetEntries();
+  for (Int_t i=0;i<300;i++){    
+    Int_t size = branch->GetEvent(i);     
+    mybranch->GetEvent(i);
+    //if (arr){
+    if ((arr->GetEntriesFast()>0) && size>0) {
+      CompareHits(arr,myhits,debug);
+      if ((i%10)==0) cout<<i<<"\n"; 
+    }    
+  }
+  gBenchmark->Stop(benchmark);
+  gBenchmark->Show(benchmark); 
+}  
+
+
+void CompareHitsG(const char * benchmark, Bool_t debug)
+{        
+  //compare persistent hits
+  TFile f2("galice.root");
+  TTree * treeh2 = (TTree*)f2.Get("TreeH0");
+  AliTPCTrackHits *  myhits = new AliTPCTrackHits ; 
+  AliObjectBranch *mybranch = (AliObjectBranch*)treeh2->GetBranch("TPC2"); 
+  mybranch->SetAddress(&myhits);
+
+    
+  TClonesArray *arr = new TClonesArray("AliTPChit",100);
+  //TFile f3("treehcl.root");
+  //TTree * treeh3 = (TTree*)f3.Get("TreeTPCH");
+  TBranch *branch = treeh2->GetBranch("TPC"); 
+  branch->SetAddress(&arr);
+
+  TFile f3("treehdelta.root","recreate");
+  f3.SetCompressionLevel(2);
+  TTree * treeh3 = new TTree("DelataH","DeltaH");
+  TClonesArray *arrd = new TClonesArray("AliTPChitD",100);
+  treeh3->Branch("TPC", &arrd,64000,2);
+   
+  cout<<"Lets go!\n"; 
+  gBenchmark->Start(benchmark);
+  Int_t trsize = treeh2->GetEntries();
+  for (Int_t i=0;i<trsize;i++){    
+    Int_t size = branch->GetEvent(i);     
+    mybranch->GetEvent(i);
+    //if (arr){
+    if ((arr->GetEntriesFast()>0) && size>0) {
+      CompareHits(arr,myhits,debug,arrd);      
+    }  
+    if ((i%10)==0) cout<<i<<"\n";
+    treeh3->Fill();         
+  }
+  treeh3->Write();
+  f3.Close();
+  gBenchmark->Stop(benchmark);
+  gBenchmark->Show(benchmark); 
+}  
+
+
+   
+AliTPCTrackHits * MakeTrack(TClonesArray * arr, TClonesArray * arrp, AliTPCTrackHits *myhits)  
+{
+  //
+  //make track wit hits  according 
+  AliTPChit * hit;
+  //  AliTPCTrackHits * myhits= new AliTPCTrackHits;
+  myhits->SetHitPrecision(0.002);
+  myhits->SetStepPrecision(0.003);  
+  myhits->SetMaxDistance(100);  
+  myhits->FlushHitStack(kTRUE);
+  for (Int_t i=0;i<arr->GetEntriesFast();i++){
+    hit = (AliTPChit*)arr->At(i);    
+    if (hit){
+      TParticle *p = (TParticle*)arrp->At(hit->GetTrack());
+      Float_t momentum = TMath::Sqrt(p->Px()*p->Px()+p->Py()*p->Py());
+      Float_t ran= 100.*gRandom->Rndm();
+      if (ran<1.)  myhits->FlushHitStack(kTRUE);
+      if (momentum<0.01) {//if small momentum - not precise
+       myhits->SetHitPrecision(0.05);
+       myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
+                            hit->Z(), hit->fQ+1000);
+      }
+      else {
+       myhits->SetHitPrecision(0.002);
+       myhits->AddHitKartez(hit->fSector,hit->GetTrack(), hit->X(), hit->Y(),
+                            hit->Z(), hit->fQ);
+      }
+    }
+  }
+  myhits->FlushHitStack();
+  return myhits;  
+}
+
+
+
+void CompareHits(TClonesArray * arr, AliTPCTrackHits * myhits, Bool_t debug, TClonesArray *arrd)
+{     
+  //
+  // if debug option  kTRUE
+  // compare hits and write result  to the stdoutput
+  AliTPChit * hit, *hit2;
+  if (arrd) arrd->Clear();
+
+  for (Int_t i=0;i<arr->GetEntriesFast();i++){
+    hit = (AliTPChit*)arr->At(i);  
+    if (hit){
+      if (i==0) myhits->First();
+      else myhits->Next();
+      hit2 = myhits->GetHit();   
+      if (!hit2) {
+       hit2=0;
+       if (hit) cout<<"Error _ hits "<<i<<"didn't find\n";
+      }
+
+      if (hit&&hit2){
+       // if (hit){
+       
+       AliTrackHitsParam *param= myhits->GetParam();
+       AliHitInfo * info = myhits->GetHitInfo();
+       //      
+       if (arrd) {
+         TClonesArray &larrd = *arrd;
+         AliTPChitD * h = new(larrd[i]) AliTPChitD;
+         h->SetX(hit->X());
+         h->SetY(hit->Y());
+         h->SetZ(hit->Z());
+         h->SetTrack(hit->GetTrack());
+         h->fQ = hit->fQ;
+         h->fSector = hit->fSector;
+         AliTPChit * hitd = h->GetDelta();
+         hitd->SetX(hit->X()-hit2->X());
+         hitd->SetY(hit->Y()-hit2->Y());
+         hitd->SetZ(hit->Z()-hit2->Z());
+         hitd->SetTrack(hit->GetTrack()-hit2->GetTrack());
+         hitd->fQ = hit->fQ-hit2->fQ;
+         hitd->fSector = hit->fSector-hit2->fSector;
+       }
+
+       if (debug){
+         Float_t dd = 
+           TMath::Sqrt(
+                       (hit->X()-hit2->X())*(hit->X()-hit2->X())+
+                       (hit->Y()-hit2->Y())*(hit->Y()-hit2->Y())+
+                       (hit->Z()-hit2->Z())*(hit->Z()-hit2->Z())); 
+         printf("C1\t%d\t%d\t%d\t%d\t",
+                hit->fSector,hit2->fSector,hit->GetTrack(),hit2->GetTrack());
+         printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
+                hit->X(),hit2->X(), hit->Y(),hit2->Y(),hit->Z(),hit2->Z());
+         printf("%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t%3.6f\t",
+                dd, param->fR,param->fZ,param->fFi,param->fAn,
+                param->fAd,param->fTheta);
+         printf("%d\t%d\t%d\n",
+                (Int_t)info->fHitDistance, (Int_t)hit->fQ, (Int_t)hit2->fQ);
+       }
+      }        
+    }
+  }
+} 
+
+
+
+
+//Filter for paw visualisation of results
+//
+//sed filter  
+ //sed '/*/d' dd.txt | sed '/C/d'|  cat -n > out.txt
+ // sed -n '/C2/p' dd.txt | sed 's/C1//' | cat -n >out2.txt 
+ // sed -n '/C3/p' dd.txt | sed 's/C2//' | cat -n >out3.txt 
+
+//filter   
+//myfilter C1 dd.txt >out1.txt
+//myfilter C2 dd.txt >out2.txt
+//myfilter C3 dd.txt >out3.txt
+// sed -n 1,50000p
+// sed -n 1,50000p dd.txt | myfilter C1  | sed -n 1,50000p >out1.txt
+
+//myfilter C1 dd.txt | sed -n 1,50000p >out1.txt
+//myfilter C2 dd.txt | sed -n 1,50000p >out2.txt
+//myfilter C3 dd.txt | sed -n 1,50000p >out3.txt 
+/*
+paw visualisation
+ Nt/Create 1 'Test of Match' 21 ! ! event v1 v2 t1 t2  x1 x2 y1  y2 z1 z2 dd r z fi an ad theta dist q1 q2
+  Nt/Read 1 out1.txt 
+ nt 
+ Nt/Create 2 'Test of Match2' 13 ! ! event stacksize i r  z fi alfa theta dr1 ddr ddz ddrfi dd
+  Nt/Read 2 out2.txt
+ Nt/Create 3 'Test of Match2' 2 ! ! event stacksize 
+  Nt/Read 3 out3.txt     
+*/ 
+
+void   Fit2(Double_t fSumY, Double_t fSumYX, Double_t fSumYX2,
+           Double_t fSumX,  Double_t fSumX2, Double_t fSumX3, 
+           Double_t fSumX4, Int_t n,
+           Double_t &a, Double_t &b, Double_t &c)
+{
+  //
+  //recalc parameters not fixing origin point
+  Double_t det = 
+    n* (fSumX2*fSumX4-fSumX3*fSumX3) -
+    fSumX*      (fSumX*fSumX4-fSumX3*fSumX2)+
+    fSumX2*     (fSumX*fSumX3-fSumX2*fSumX2);
+    
+  if (TMath::Abs(det)>0){    
+    a = 
+      (fSumY * (fSumX2*fSumX4-fSumX3*fSumX3)-
+       fSumX *(fSumYX*fSumX4-fSumYX2*fSumX3)+
+       fSumX2*(fSumYX*fSumX3-fSumYX2*fSumX2))/det; 
+    b=
+      (n*(fSumYX*fSumX4-fSumX3*fSumYX2)-
+      fSumY*(fSumX*fSumX4-fSumX3*fSumX2)+
+      fSumX2*(fSumX*fSumYX2-fSumYX*fSumX2))/det;
+    c=
+      (n*(fSumX2*fSumYX2-fSumYX*fSumX3)-
+       fSumX*(fSumX*fSumYX2-fSumYX*fSumX2)+
+       fSumY*(fSumX*fSumX3-fSumX2*fSumX2))/det;  
+    cout<<a<<"\t"<<b<<"\t"<<c<<"\n";
+  }
+}
+
+void TestFit(Float_t a, Float_t b, Float_t c, Int_t n)
+{
+  Double_t fSumY,fSumYX,fSumYX2,fSumX, fSumX2,fSumX3, fSumX4;
+  fSumY = fSumYX = fSumYX2 = fSumX =  fSumX2 = fSumX3 =  fSumX4 =0;
+  Double_t a2,b2,c2;
+  for (Int_t i=0;i<n; i++){
+    Float_t  x = gRandom->Rndm();
+    Float_t  y = a+b*x+c*x*x;
+    fSumY+=y;
+    fSumYX+=y*x;
+    fSumYX2+=y*x*x;
+    fSumX += x;
+    fSumX2 += x*x;
+    fSumX3 += x*x*x;
+    fSumX4 += x*x*x*x;
+  }
+  Fit2(fSumY,fSumYX,fSumYX2, fSumX,fSumX2,fSumX3,fSumX4,n,a2,b2,c2);  
+}
diff --git a/TPC/TestTPCTrackHits.h b/TPC/TestTPCTrackHits.h
new file mode 100644 (file)
index 0000000..8a8fb7c
--- /dev/null
@@ -0,0 +1,10 @@
+void ConvertHits(const char * benchmark="0", Bool_t debug=kFALSE);
+void CompareHits(const char * benchmark="1", Bool_t debug=kFALSE);
+
+class AliTPChitD : public AliTPChit {
+public:
+  AliTPChit * GetDelta() {return &fDelta;}
+private:
+  AliTPChit fDelta;     //delta of hit information 
+  ClassDef(AliTPChitD,1) 
+};