+//#include "alles.h"
struct GoodTrack {
Int_t lab;
Int_t code;
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");
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
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);
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);
}
+
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);
hep->SetXTitle("p (Gev/c)"); hep->SetYTitle("dE/dX (Arb. Units)");
hep->Draw(); c1->cd();
+ gBenchmark->Stop("AliTPCComparison");
+ gBenchmark->Show("AliTPCComparison");
+
+
return 0;
}
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"))) {
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:
{
}
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;
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;
}
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;
}
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);
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;}
}
delete gAlice; gAlice=0;
file->Close();
+ gBenchmark->Stop("AliTPCComparison");
+ gBenchmark->Show("AliTPCComparison");
return nt;
}
--- /dev/null
+///////////////////////////////////////////////////////////////////////////////
+// //
+// 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 ¶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"<<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 ¶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<<" "<<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;
+}
+
+
+
--- /dev/null
+#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
--- /dev/null
+#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;
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;
}
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;
# 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
# 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
# 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)
--- /dev/null
+/*
+ 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);
+}
--- /dev/null
+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)
+};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff