From cf98c13f3ecf1714b3fd385672b41dd1c28d1a42 Mon Sep 17 00:00:00 2001 From: kowal2 Date: Thu, 2 Nov 2000 07:40:20 +0000 Subject: [PATCH] New classes for handling and testing new hit structures --- TPC/AliTPCComparison.C | 131 ++++-- TPC/AliTPCTrackHits.cxx | 678 ++++++++++++++++++++++++++++++++ TPC/AliTPCTrackHits.h | 104 +++++ TPC/AliTPCTrackHitsInterfaces.h | 110 ++++++ TPC/AliTPCtest.C | 8 +- TPC/Makefile | 16 +- TPC/TestTPCTrackHits.cxx | 335 ++++++++++++++++ TPC/TestTPCTrackHits.h | 10 + 8 files changed, 1359 insertions(+), 33 deletions(-) create mode 100644 TPC/AliTPCTrackHits.cxx create mode 100644 TPC/AliTPCTrackHits.h create mode 100644 TPC/AliTPCTrackHitsInterfaces.h create mode 100644 TPC/TestTPCTrackHits.cxx create mode 100644 TPC/TestTPCTrackHits.h diff --git a/TPC/AliTPCComparison.C b/TPC/AliTPCComparison.C index 2329373d80c..74d540ac18f 100644 --- a/TPC/AliTPCComparison.C +++ b/TPC/AliTPCComparison.C @@ -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; iLoadEntry(i); + for (i=0; iLoadEntry(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; iSetAddress(&iotrack); @@ -93,29 +98,60 @@ Int_t AliTPCComparison() { TH1F *he =new TH1F("he","dE/dX for pions with 0.4Fill(ptg); - AliTPCtrack *track; - Int_t i; + AliTPCtrack *track=0; for (i=0; iUncheckedAt(i); + track=(AliTPCtrack*)tarray.UncheckedAt(i); tlab=track->GetLabel(); if (lab==TMath::Abs(tlab)) break; } if (i==nentr) { - cerr<<"Track "<GetLabel())) micount++; + } + if (micount>1) { + //cout<<"Track no. "<GetPadRowRadii(0,0); track->PropagateTo(xk); if (lab==tlab) hfound->Fill(ptg); - else { hfake->Fill(ptg); cerr<Fill(ptg); + //cerr<GetPt());track->GetPxPyPz(px,py,pz); @@ -143,11 +179,35 @@ Int_t AliTPCComparison() { } + Stat_t ng=hgood->GetEntries(); cerr<<"Good tracks "<GetEntries(); if (ng!=0) - cerr<<"Integral efficiency is about "<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; iiGet("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; iGetEntries(); + sectors_by_rows=(Int_t)TD->GetEntries(); for (i=0; iGetEvent(i)) continue; Int_t sec,row; @@ -306,9 +376,9 @@ Int_t good_tracks(GoodTrack *gt, Int_t max) { for (Int_t j=0; j1) { 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<Close(); + gBenchmark->Stop("AliTPCComparison"); + gBenchmark->Show("AliTPCComparison"); return nt; } diff --git a/TPC/AliTPCTrackHits.cxx b/TPC/AliTPCTrackHits.cxx new file mode 100644 index 00000000000..aadb4a58569 --- /dev/null +++ b/TPC/AliTPCTrackHits.cxx @@ -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 +/* + +*/ +//End_Html +// // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include "TVector3.h" +#include "TClonesArray.h" +#include "AliTPCTrackHits.h" +#include "AliTPC.h" + +#include + + + +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)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)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)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)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 ¶m = *(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"<fStackIndex<<"\n"<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 ¶m = *(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<<" "<fStackIndex<<" \t"; + cout<<" "<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+1GetSize()) + &&((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 index 00000000000..edd0c0b3527 --- /dev/null +++ b/TPC/AliTPCTrackHits.h @@ -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 index 00000000000..0623672aeb8 --- /dev/null +++ b/TPC/AliTPCTrackHitsInterfaces.h @@ -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;iDump();} +}; +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;iDump();} +}; +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;iDump();} +}; +AliClassAliHitInfo galiclass____AliHitInfo; diff --git a/TPC/AliTPCtest.C b/TPC/AliTPCtest.C index f6476b27105..7f6d197002c 100644 --- a/TPC/AliTPCtest.C +++ b/TPC/AliTPCtest.C @@ -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; diff --git a/TPC/Makefile b/TPC/Makefile index 00ef48fb143..4318578eaee 100644 --- a/TPC/Makefile +++ b/TPC/Makefile @@ -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 index 00000000000..846d108cf57 --- /dev/null +++ b/TPC/TestTPCTrackHits.cxx @@ -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<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<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;iGetEvent(i); + mybranch->GetEvent(i); + //if (arr){ + if ((arr->GetEntriesFast()>0) && size>0) { + CompareHits(arr,myhits,debug,arrd); + } + if ((i%10)==0) cout<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;iGetEntriesFast();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;iGetEntriesFast();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 "<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<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 index 00000000000..8a8fb7cfc5f --- /dev/null +++ b/TPC/TestTPCTrackHits.h @@ -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) +}; -- 2.43.5