/*
$Log$
+Revision 1.2 2001/02/05 14:43:13 hristov
+Compare() declared const
+
Revision 1.1 2000/10/05 16:17:27 kowal2
New class replacing AliCluster
}
ClassImp(AliDigitCluster)
-ClassImp(AliDifCluster)
+ClassImp(AliTPCExactPoint)
+ClassImp(AliTPCClusterPoint)
+ClassImp(AliTPCTrackerPoint)
+ClassImp(AliTPCTrackPoint)
+ClassImp(AliTPCTrackPointRef)
+
+
/* $Id$ */
#include "TObject.h"
+#include "TMath.h"
class AliComplexCluster : public TObject {
public:
public:
AliComplexCluster() {
fTracks[0]=fTracks[1]=fTracks[2]=0;
- fX=fY=fQ=fSigmaX2=fSigmaY2=0.;
+ fX=fY=fQ=fSigmaX2=fSigmaY2=fSigmaXY=fArea=fMax=0.;
}
virtual ~AliComplexCluster() {;}
Bool_t IsSortable() const;
ClassDef(AliDigitCluster,1) // Tclusters
};
-class AliDifCluster : public AliDigitCluster {
-public:
- Float_t fDx; //delta x
- Float_t fDy; //delta y
- Float_t fAngleX;//x angle
- Float_t fAngleY;//y angle
- Float_t fOrigQ; //original charge
- Int_t fGenerTrack; //track number of nearest track
- ClassDef(AliDifCluster,1) // // Cluster manager
+
+
+class AliTPCTrackerPoint {
+ private:
+ Short_t fTX; // x position of the cluster in cm - 10 mum prec
+ Short_t fTZ; // current prolongation in Z in cm - 10 mum prec.
+ Short_t fTY; // current prolongation in Y in cm - 10 mum prec.
+ Char_t fTAngleZ; // angle
+ Char_t fTAngleY; // angle
+ public:
+ AliTPCTrackerPoint(){fTX=0; fTY=0; fTZ=0; fTAngleZ=0; fTAngleY=0;}
+ inline Float_t GetX() {return (fTX*0.01);}
+ inline Float_t GetZ() {return (fTZ*0.01);}
+ inline Float_t GetY() {return (fTY*0.01);}
+ inline Float_t GetAngleZ() {return (Float_t(fTAngleZ)*0.02);}
+ inline Float_t GetAngleY() {return (Float_t(fTAngleY)*0.02);}
+ //
+ void SetX(Float_t x){ fTX = Short_t(TMath::Nint(x*100.));}
+ void SetY(Float_t y){ fTY = Short_t(TMath::Nint(y*100.));}
+ void SetZ(Float_t z){ fTZ = Short_t(TMath::Nint(z*100.));}
+ void SetAngleZ(Float_t anglez) {fTAngleZ = Char_t(TMath::Nint(anglez*50.));}
+ void SetAngleY(Float_t angley) {fTAngleY = Char_t(TMath::Nint(angley*50.));}
+ //
+ public:
+ ClassDef(AliTPCTrackerPoint,1)
};
+class AliTPCClusterPoint {
+ private:
+ Short_t fCZ; // current cluster position Z in cm - 100 mum precision
+ Short_t fCY; // current cluster position Y in cm - 100 mum precision
+ UChar_t fErrZ; // z error estimate - in mm - 50 mum precision
+ UChar_t fErrY; // y error estimate - in mm - 50 mum precision
+ UChar_t fSigmaZ; // shape Z - normalised shape - normaliziation 1 - precision 2 percent
+ UChar_t fSigmaY; // shape Y - normalised shape - normaliziation 1 - precision 2 percent
+ UShort_t fQ; // total charge in cluster
+ UShort_t fMax; // charge at maximum
+ Char_t fCType; // type of the cluster
+ public:
+ AliTPCClusterPoint(){fCZ=fCY=fSigmaZ=fSigmaY=fErrZ=fErrY=fQ=fMax=fCType=0;}
+ inline Float_t GetZ() {return (fCZ*0.01);}
+ inline Float_t GetY() {return (fCY*0.01);}
+ inline Float_t GetErrZ() {return (fErrZ*0.005);}
+ inline Float_t GetErrY() {return (fErrY*0.005);}
+ inline Float_t GetSigmaZ() {return (fSigmaZ*0.02);}
+ inline Float_t GetSigmaY() {return (fSigmaY*0.02);}
+ inline Int_t GetType() {return fCType;}
+ inline Int_t GetMax() {return fMax;}
+ inline Float_t GetQ() {return fQ;}
+
+ //
+ void SetY(Float_t y){ fCY = Short_t(TMath::Nint(y*100.));}
+ void SetZ(Float_t z){ fCZ = Short_t(TMath::Nint(z*100.));}
+ void SetErrZ(Float_t errz) {fErrZ = UChar_t(TMath::Nint(errz*200.));}
+ void SetErrY(Float_t erry) {fErrY = UChar_t(TMath::Nint(erry*200.));}
+ void SetSigmaZ(Float_t sigmaz) {fSigmaZ = UChar_t(TMath::Nint(sigmaz*50.));}
+ void SetSigmaY(Float_t sigmay) {fSigmaY = UChar_t(TMath::Nint(sigmay*50.));}
+ void SetQ(Float_t q) {fQ = UShort_t(q);}
+ void SetMax(Float_t max) {fMax = UShort_t(max);}
+ void SetType(Char_t type) {fCType = type;}
+
+ //
+ public:
+ ClassDef(AliTPCClusterPoint,1)
+};
+
+
+class AliTPCExactPoint : public TObject{
+ public:
+ AliTPCExactPoint(){fEZ=fEY=fEAngleZ=fEAngleY=fEAmp=fEPrim=fTrackID=0;}
+ Float_t fEZ; // current "exact" position according simulation
+ Float_t fEY; // current "exact" position according simulation
+ Float_t fEAngleZ; // angle Z
+ Float_t fEAngleY; // angle Y
+ Float_t fEAmp; // total charge deposited in row
+ Float_t fEPrim; // primary charge deposited in row
+ Int_t fTrackID; // id of the track
+ ClassDef(AliTPCExactPoint,1)
+};
+
+
+class AliTPCTrackPoint: public TObject{
+ public:
+ AliTPCTrackPoint(){ fIsShared = kFALSE;}
+ AliTPCClusterPoint & GetCPoint(){return fCPoint;}
+ AliTPCTrackerPoint & GetTPoint(){return fTPoint;}
+ public:
+ AliTPCClusterPoint fCPoint;
+ AliTPCTrackerPoint fTPoint;
+ Bool_t fIsShared; // indicate sharing of the point between several tracks
+ ClassDef(AliTPCTrackPoint,1)
+};
+
+
+class AliTPCTrackPointRef: public AliTPCTrackPoint{
+ public:
+ AliTPCExactPoint & GetExactPoint(){return fEPoint;}
+ AliTPCExactPoint & GetNearestPoint(){return fNPoint;}
+ public:
+ AliTPCExactPoint fEPoint; //exact point belonging to track
+ AliTPCExactPoint fNPoint; //nearest point
+ ClassDef(AliTPCTrackPointRef,1)
+};
+
+
#endif //ALICOMPLEXCLUSTER_H
/*
$Log$
+Revision 1.4 2000/10/05 16:04:21 kowal2
+Forward declarations
+
Revision 1.3 2000/06/30 12:07:49 kowal2
Updated from the TPC-PreRelease branch
#include "TClonesArray.h"
#include "AliTPC.h"
#include "TRandom.h"
-#include "AliTPCClusterFinder.h"
ClassImp(AliH2F)
TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
{
- //find peaks and write it in form of AliTPCcluster to array
-
- //firstly we need to create object for cluster finding
- //and fill it with contents of histogram
- AliTPCClusterFinder cfinder;
- cfinder.SetThreshold(threshold);
- cfinder.SetNoise(noise);
- cfinder.GetHisto(this);
- return cfinder.FindPeaks3();
+
+ //
+ // not implemented
+ //
+ return 0;
}
void AliH2F::ClearSpectrum()
+/****************************************************************************
+ * Origin: M.Ivanov, CERN, Marian.Ivanov@cern.ch *
+ ****************************************************************************/
+
+
+
#ifndef __CINT__
#include "alles.h"
#include "AliComplexCluster.h"
--- /dev/null
+#ifndef __CINT__
+#include "alles.h"
+#include "AliComplexCluster.h"
+//#include "AliTPCclusterM.h"
+#include "AliTPCtrackerMI.h"
+#include "AliTPCclusterMI.h"
+#endif
+
+
+Int_t AliTPCCompareTracks(Int_t eventn, Bool_t all = kFALSE) {
+ cerr<<"Comparing tracks...\n";
+ //CONNECT FILES
+ TFile *file=TFile::Open("galice.root");
+ if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 1;}
+ //
+ TFile *ftracks=TFile::Open("AliTPCtracks.root","update");
+ if (!ftracks->IsOpen()){cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
+ //
+ AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60_150x60");
+ if (!param) {cerr<<"TPC parameters have not been found !\n"; return 2;}
+ //
+ TFile *fout=TFile::Open("AliTPCTracksDif.root","new");
+ if (!fout->IsOpen()) fout=TFile::Open("AliTPCClustersDif.root","recreate");
+ //
+ //connect exact clusters
+ file->cd();
+ AliTPCClustersArray *ce=new AliTPCClustersArray;
+ ce->Setup(param);
+ ce->SetClusterType("AliComplexCluster");
+ char cname[100];
+ sprintf(cname,"TreeCExact_%d",eventn);
+ ce->ConnectTree(cname);
+ //
+ //connect reconstructed tracks
+ ftracks->cd();
+ sprintf(cname,"Seeds");
+ TTree * treetracks = (TTree*)ftracks->Get(cname);
+ TBranch * branchtracks = treetracks->GetBranch("seeds");
+ //
+ //load seeds to the memory
+ Int_t trackmap[500000][4]; // map of tracks corresponding to given track number
+ memset(trackmap,0,sizeof(Int_t)*4*500000);
+ Int_t ntracks = treetracks->GetEntries();
+ TObjArray * trackarr= new TObjArray(ntracks);
+ Int_t nproces = TMath::Min(ntracks,4000);
+
+ // for (Int_t itrack =0; itrack<ntracks; itrack++){
+ for (Int_t itrack =0; itrack<nproces; itrack++){
+ AliTPCseed * seed = new AliTPCseed;
+ //
+ seed->fPoints = new TClonesArray("AliTPCTrackPoint",200);
+ branchtracks->SetAddress(&seed);
+ branchtracks->GetEntry(itrack);
+
+ //if (seed->fRemoval>0 && (itrack%4) ) continue;
+ trackarr->AddLast(seed);
+
+ //crete array with exact position information
+ seed->fEPoints = new TClonesArray("AliTPCTrackPointRef",1);
+ seed->fEPoints->ExpandCreateFast(200);
+
+ Int_t label = TMath::Abs(seed->GetLabel());
+ Int_t i;
+ if (label>500000) {
+ //too big track label
+ }else{
+ for (i=0;i<4;i++) {
+ if ( trackmap[label][i]==0)
+ break;
+ }
+ if(i<4) trackmap[label][i]=itrack;
+ }
+ }
+
+ //add information about exact positions
+ Int_t nrows=Int_t(ce->GetTree()->GetEntries());
+
+ for (Int_t n=0; n<nrows; n++) {
+ AliSegmentID *se=ce->LoadEntry(n);
+ Int_t sec,row;
+ param->AdjustSectorRow(se->GetID(),sec,row);
+ //
+ AliTPCClustersRow* clrowe = ce->GetRow(sec,row);
+ //
+ Int_t ncl=clrowe->GetArray()->GetEntriesFast();
+ const Int_t kNIS=param->GetNInnerSector(), kNOS=param->GetNOuterSector();
+ Int_t npads, sign;
+ if (sec < kNIS) {
+ npads = param->GetNPadsLow(row);
+ sign = (sec < kNIS/2) ? 1 : -1;
+ } else {
+ npads = param->GetNPadsUp(row);
+ sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
+ }
+
+ Int_t trrow = row;
+ if (sec>=param->GetNInnerSector()) trrow += param->GetNRowLow();
+
+ while (ncl--) {
+ AliComplexCluster *cl=(AliComplexCluster*)(clrowe->GetArray()->At(ncl));
+ Double_t x=param->GetPadRowRadii(sec,row);
+ Double_t y=cl->fY;
+ y = ( y + 0.5 - 0.5*npads) *param->GetPadPitchWidth(sec);
+ Double_t z=cl->fX*param->GetZWidth();
+ z = sign*(param->GetZLength() - z);
+ Float_t cs, sn, tmp;
+ param->AdjustCosSin(sec,cs,sn);
+ for (Int_t i=0;i<3;i++){
+ if (cl->fTracks[0]<500000) if (trackmap[cl->fTracks[0]][i]) {
+ AliTPCseed * seed = (AliTPCseed*)trackarr->At(trackmap[cl->fTracks[0]][i]);
+ TClonesArray * clarr = seed->fPoints;
+ if (!clarr){
+ //printf("Seed %d without info",trackmap[cl->fTracks[0]][i]);
+ continue;
+ }
+ AliTPCTrackPoint * trcluster = (AliTPCTrackPoint*)(seed->fPoints->At(trrow));
+ AliTPCTrackPointRef * ecluster = (AliTPCTrackPointRef*)(seed->fEPoints->At(trrow));
+ //
+ ecluster->GetCPoint() = trcluster->GetCPoint();
+ ecluster->GetTPoint() = trcluster->GetTPoint();
+ //
+ AliTPCExactPoint & epoint = ecluster->GetExactPoint();
+ /*
+ trcluster->fEZ = z;
+ trcluster->fEY = y;
+ trcluster->fAngleY = cl->fSigmaY2*param->GetPadPitchLength(sec);
+ trcluster->fAngleZ = cl->fSigmaX2*param->GetPadPitchLength(sec);
+ trcluster->fEAmp = cl->fQ;
+ trcluster->fEPrim = cl->fMax;
+ */
+ epoint.fEZ = z;
+ epoint.fEY = y;
+ epoint.fEAngleY = cl->fSigmaY2*param->GetPadPitchLength(sec);
+ epoint.fEAngleZ = cl->fSigmaX2*param->GetPadPitchLength(sec);
+ epoint.fEAmp = cl->fQ;
+ epoint.fEPrim = cl->fMax;
+ }
+ }
+ }
+ // cr->ClearRow(sec,row);
+ ce->ClearRow(sec,row);
+ }
+
+ // make new tree - with tracks - exact position
+ fout->cd();
+ TTree * treenew = new TTree("Seedref","Seedref");
+ AliTPCseed * ioseed = (AliTPCseed*)trackarr->At(0);
+ TBranch * br = treenew->Branch("seeds","AliTPCseed",&ioseed,32000,99);
+
+ for (Int_t itrack =0; itrack<ntracks; itrack++){
+ ioseed = (AliTPCseed*) trackarr->At(itrack);
+ br->SetAddress(&ioseed);
+ treenew->Fill();
+ }
+
+ treenew->Write();
+ printf("1\n");
+ delete ce;
+ printf("2\n");
+ //delete cr;
+ // carray->GetTree()->Write();
+ printf("3\n");
+ // delete carray;
+ printf("4\n");
+ // cf->Close();
+ printf("5\n");
+ fout->Close();
+ printf("6\n");
+ ftracks->Close();
+ printf("7\n");
+ return 0;
+}
#ifndef __CINT__
#include "alles.h"
#include "AliTPCtracker.h"
+ #include <AliMagF.h>
#endif
struct GoodTrackTPC {
Int_t AliTPCComparison(Int_t event=0) {
cerr<<"Doing comparison...\n";
-
+ AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
const Int_t MAX=20000;
Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event);
break;
case 2:
{
+ TFile fff("digits.root");
+
char dname[100]; sprintf(dname,"TreeD_75x40_100x60_150x60_%d",event);
TTree *TD=(TTree*)gDirectory->Get(dname);
AliSimDigits da, *digits=&da;
cerr<<"gAlice have not been found on galice.root !\n";
exit(5);
}
-
+
+ TFile *fdigit = TFile::Open("digits.root");
+ file->cd();
AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
Int_t ver = TPC->IsVersion();
case 2:
{
sprintf(treeName,"TreeD_75x40_100x60_150x60_%d",event);
- TD=(TTree*)gDirectory->Get(treeName);
+ // TD=(TTree*)gDirectory->Get(treeName);
+ TD=(TTree*)fdigit->Get(treeName);
TD->GetBranch("Segment")->SetAddress(&digits);
count = new Int_t[np];
Int_t i;
Int_t sec,row;
digp->AdjustSectorRow(digits->GetID(),sec,row);
digits->First();
+ digits->ExpandTrackBuffer();
do { //Many thanks to J.Chudoba who noticed this
Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
- Short_t dig = digits->GetDigit(it,ip);
- Int_t idx0=digits->GetTrackID(it,ip,0);
- Int_t idx1=digits->GetTrackID(it,ip,1);
- Int_t idx2=digits->GetTrackID(it,ip,2);
+ // Short_t dig = digits->GetDigit(it,ip);
+ Short_t dig = digits->CurrentDigit();
+
+ Int_t idx0=digits->GetTrackIDFast(it,ip,0)-2;
+ Int_t idx1=digits->GetTrackIDFast(it,ip,1)-2;
+ Int_t idx2=digits->GetTrackIDFast(it,ip,2)-2;
if (idx0>=0 && dig>=zero && idx0<np ) count[idx0]+=1; //background events - higher track ID's
if (idx1>=0 && dig>=zero && idx1<np ) count[idx1]+=1;
if (idx2>=0 && dig>=zero && idx2<np ) count[idx2]+=1;
for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
part = gAlice->Particle(iprim);
+ Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+ if (ptot<0.01) continue;
char *xxx=strstr(part->GetName(),"XXX");
if(xxx)continue;
if(part->T()>timecut)continue;
- Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+ //Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
Bool_t prmch = kTRUE; // candidate primary track
TFile *in=TFile::Open("rfio:galice.root");
if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
+ TFile *ind=TFile::Open("digits.root");
+ if (!ind->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
+
if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
cerr<<"gAlice have not been found on galice.root !\n";
return 3;
{
// delete gAlice; gAlice=0;
AliTPCv2 tpc;
- tpc.SetParam(dig); timer.Start(); cwd->cd();
+ tpc.SetParam(dig); timer.Start(); ind->cd();
for (Int_t i=0;i<n;i++){
printf("Processing event %d\n",i);
tpc.Digits2Clusters(out,i);
--- /dev/null
+/****************************************************************************
+ * Origin: M.Ivanov marian.ivanov@cern.ch *
+ ****************************************************************************/
+
+/*
+
+ macro to create array of clusters from TPC digits
+ input files - galice.root
+ digits.root - file with digits - usualy use link to galice.root
+ - in splitted mode - neccesary to create link to proper file
+
+ output file - AliTPCclusters.root
+ - to be used by AliTPCTrackFinderMI.C
+
+ Warning - if cluster file AliTPCclusters.root already exist - macro exit and don't produce anything
+
+
+*/
+
+
+#ifndef __CINT__
+#include <iostream.h>
+#include "AliRun.h"
+#include "AliTPCv1.h"
+#include "AliTPCv2.h"
+#include "AliTPCParam.h"
+#include "AliTPCclustererMI.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "TTree.h"
+#endif
+
+Int_t AliTPCFindClustersMI(Int_t n=1) {
+ TFile *out=TFile::Open("AliTPCclusters.root","new");
+ if (!out->IsOpen()) {cerr<<"Delete old AliTPCclusters.root !\n"; return 1;}
+ TFile *in=TFile::Open("rfio:galice.root");
+ if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
+
+ TFile *ind=TFile::Open("rfio:digits.root");
+ if (!ind->IsOpen()) {cerr<<"Can't open digits file !\n"; return 2;}
+
+
+ if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+ cerr<<"gAlice have not been found on galice.root !\n";
+ return 3;
+ }
+
+ TDirectory *cwd = gDirectory;
+
+ AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
+ Int_t ver = TPC->IsVersion();
+ cerr<<"TPC version "<<ver<<" has been found !\n";
+
+ AliTPCParam *dig=(AliTPCParam *)in->Get("75x40_100x60_150x60");
+ if (!dig) {cerr<<"TPC parameters have not been found !\n"; return 4;}
+
+ TStopwatch timer;
+
+ switch (ver) {
+ case 1:
+ cerr<<"Making clusters...\n";
+ {
+ AliTPCv1 &tpc=*((AliTPCv1*)TPC);
+ tpc.SetParam(dig); timer.Start(); cwd->cd();
+ for(Int_t i=0;i<n;i++){
+ printf("Processing event %d\n",i);
+ gAlice->GetEvent(i);
+ tpc.Hits2Clusters(out,i);
+ }
+ }
+ break;
+ case 2:
+ cerr<<"Looking for clusters...\n";
+ {
+ // delete gAlice; gAlice=0;
+ AliTPCv2 tpc;
+ tpc.SetParam(dig); timer.Start(); cwd->cd();
+
+
+ for (Int_t i=0;i<n;i++){
+ AliTPCclustererMI clusterer;
+ char dname[100];
+ char cname[100];
+ sprintf(dname,"TreeD_75x40_100x60_150x60_%d",i);
+ sprintf(cname,"TreeC_TPC_%d",i);
+ TTree * input = (TTree*)ind->Get(dname);
+ out->cd();
+ TTree * output = new TTree(cname,cname);
+
+ printf("Processing event %d\n",i);
+ clusterer.SetInput(input);
+ clusterer.SetOutput(output);
+ clusterer.Digits2Clusters(dig, i);
+ //tpc.Digits2Clusters(out,i);
+ // AliTPCclusterer::Digits2Clusters(dig, out, i);
+ }
+ }
+ break;
+ default:
+ cerr<<"Invalid TPC version !\n";
+ return 5;
+ }
+
+ timer.Stop(); timer.Print();
+
+ delete gAlice; gAlice=0;
+
+ out->Close();
+
+ in->Close();
+
+ return 0;
+}
Int_t AliTPCFindTracks(Int_t eventn=1) {
cerr<<"Looking for tracks...\n";
+ TFile f("galice.root");
+ gAlice = (AliRun*)f.Get("gAlice");
TFile *out=TFile::Open("AliTPCtracks.root","new");
if (!out->IsOpen()) {cerr<<"Delete old AliTPCtracks.root !\n"; return 1;}
--- /dev/null
+/****************************************************************************
+ * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
+ ****************************************************************************/
+
+#ifndef __CINT__
+#include <iostream.h>
+#include "AliTPCParam.h"
+#include "AliTPCtrackerMI.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#endif
+
+Int_t AliTPCFindTracks(Int_t eventn=1) {
+ cerr<<"Looking for tracks...\n";
+ TFile f("galice.root");
+ gAlice = (AliRun*)f.Get("gAlice");
+ AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
+ TFile *out=TFile::Open("AliTPCtracks.root","new");
+ if (!out->IsOpen()) {cerr<<"Delete old AliTPCtracks.root !\n"; return 1;}
+
+ TFile *in=TFile::Open("AliTPCclusters.root");
+ if (!in->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 2;}
+
+ AliTPCParam *par=(AliTPCParam*)in->Get("75x40_100x60_150x60");
+ if (!par) {cerr<<"Can't get TPC parameters !\n"; return 3;}
+
+ TStopwatch timer;
+
+ Int_t rc=0;
+ for (Int_t i=0;i<eventn;i++){
+ printf("Processing event %d\n",i);
+ AliTPCtrackerMI *tracker = new AliTPCtrackerMI(par,i);
+ //delete tracker;
+ //tracker = new AliTPCtrackerMI(par,i);
+ //Double_t xyz[]={0.,0.,0.}; tracker->SetVertex(xyz); //primary vertex
+ rc=tracker->Clusters2Tracks(0,out);
+ delete tracker;
+ }
+ timer.Stop(); timer.Print();
+
+ delete par; //Thanks to Mariana Bondila
+
+ in->Close();
+ out->Close();
+
+ return rc;
+}
+/****************************************************************************
+ * Origin: M.Ivanov, CERN, Marian.Ivanov@cern.ch *
+ ****************************************************************************/
+
Int_t AliTPCHits2ExactClusters(Int_t nevent=1)
{
+ //
+ // loop over hits and find itersection of the trak with pad rows
+ // so called exact clusters stored in the structure similar to hits
+ // it's used for comparison purpose - to have benchmark of cluster
+ // finder and tracker
- // new version by J.Belikov
// Connect the Root Galice file containing Geometry, Kine and Hits
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+#include "AliTPCclusterMI.h"
+
+ClassImp(AliTPCclusterMI)
+ClassImp(AliTPCclusterLMI)
+
+Bool_t AliTPCclusterMI::IsSortable() const
+{
+ //
+ //
+ return kTRUE;
+
+}
+
+Int_t AliTPCclusterMI::Compare(const TObject* obj) const
+{
+ //
+ // compare according y
+ AliTPCclusterMI * o2 = (AliTPCclusterMI*)obj;
+ return (o2->GetY()>fY)? -1:1;
+}
--- /dev/null
+#ifndef ALITPCCLUSTERMI_H
+#define ALITPCCLUSTERMI_H
+
+//-------------------------------------------------------
+// TPC Cluster Class
+//
+// Origin: Marian Ivanov
+//-------------------------------------------------------
+
+#include "AliCluster.h"
+#include "TMath.h"
+
+//_____________________________________________________________________________
+class AliTPCclusterMI : public AliCluster {
+public:
+ AliTPCclusterMI():AliCluster(){fQ=0;}
+ AliTPCclusterMI(Int_t *lab, Float_t *hit) : AliCluster(lab,hit) {fQ = (Short_t)hit[4];}
+ virtual Bool_t IsSortable() const;
+ virtual Int_t Compare(const TObject* obj) const;
+ void Use() {fQ=-fQ;}
+ void SetQ(Float_t q) {fQ=(Short_t)q;}
+ void SetType(Char_t type) {fType=type;}
+ void SetMax(Short_t max) {fMax=max;}
+ Int_t IsUsed() const {return (fQ<0) ? 1 : 0;}
+
+ Float_t GetQ() const {return TMath::Abs(fQ);}
+ Float_t GetMax() const {return fMax;}
+ Char_t GetType()const {return fType;}
+private:
+ Short_t fQ ; //Q of cluster (in ADC counts)
+ Char_t fType; //type of the cluster 0 means golden
+ Short_t fMax; //maximal amplitude in cluster
+
+ ClassDef(AliTPCclusterMI,1) // Time Projection Chamber clusters
+};
+
+class AliTPCclusterLMI {
+ private:
+ Short_t fCZ; // current cluster position Z in cm - 100 mum precision
+ Short_t fCY; // current cluster position Y in cm - 100 mum precision
+ UChar_t fSigmaZ; // shape Z - normalised shape - normaliziation 1 - precision 2 percent
+ UChar_t fSigmaY; // shape Y - normalised shape - normaliziation 1 - precision 2 percent
+ UShort_t fQ; // total charge in cluster
+ UShort_t fMax; // charge at maximum
+ Char_t fCType; // type of the cluster
+ Int_t fLabel[3]; // track indexes
+ public:
+ AliTPCclusterLMI(){fCZ=fCY=fSigmaZ=fSigmaY=fQ=fMax=fCType=0;}
+ inline Float_t GetZ() {return (fCZ*0.01);}
+ inline Float_t GetY() {return (fCY*0.01);}
+ inline Float_t GetSigmaZ() {return (fSigmaZ*0.02);}
+ inline Float_t GetSigmaY() {return (fSigmaY*0.02);}
+ inline Int_t GetType() {return fCType;}
+ inline Int_t GetMax() {return fMax;}
+ inline Float_t GetQ() {return fQ;}
+ inline Int_t GelLabel(Int_t i){return fLabel[i];}
+ //
+ void SetY(Float_t y){ fCY = Short_t(TMath::Nint(y*100.));}
+ void SetZ(Float_t z){ fCZ = Short_t(TMath::Nint(z*100.));}
+ void SetSigmaZ(Float_t sigmaz) {fSigmaZ = UChar_t(TMath::Nint(sigmaz*50.));}
+ void SetSigmaY(Float_t sigmay) {fSigmaY = UChar_t(TMath::Nint(sigmay*50.));}
+ void SetQ(Float_t q) {fQ = UShort_t(q);}
+ void SetMax(Float_t max) {fMax = UShort_t(max);}
+ void SetType(Char_t type) {fCType = type;}
+ void SetLabels(Int_t labels[3]){fLabel[0] = labels[0];fLabel[1] = labels[1];fLabel[2] = labels[2];}
+ //
+ public:
+ ClassDef(AliTPCclusterLMI,1)
+};
+
+
+
+#endif
+
+
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+
+//-------------------------------------------------------
+// Implementation of the TPC clusterer
+//
+// Origin: Marian Ivanov
+//-------------------------------------------------------
+
+#include "AliTPCclustererMI.h"
+#include "AliTPCclusterMI.h"
+#include <TObjArray.h>
+#include <TFile.h>
+#include "AliTPCClustersArray.h"
+#include "AliTPCClustersRow.h"
+#include "AliDigits.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include <iostream.h>
+#include <TTree.h>
+
+ClassImp(AliTPCclustererMI)
+
+
+
+AliTPCclustererMI::AliTPCclustererMI()
+{
+ fInput =0;
+ fOutput=0;
+}
+void AliTPCclustererMI::SetInput(TTree * tree)
+{
+ //
+ // set input tree with digits
+ //
+ fInput = tree;
+ if (!fInput->GetBranch("Segment")){
+ cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
+ fInput=0;
+ return;
+ }
+}
+
+void AliTPCclustererMI::SetOutput(TTree * tree)
+{
+ //
+ //
+ fOutput= tree;
+ AliTPCClustersRow clrow;
+ AliTPCClustersRow *pclrow=&clrow;
+ clrow.SetClass("AliTPCclusterMI");
+ clrow.SetArray(1); // to make Clones array
+ fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
+}
+
+
+Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){
+ // sigma y2 = in digits - we don't know the angle
+ Float_t z = iz*fParam->GetZWidth();
+ Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
+ (fPadWidth*fPadWidth);
+ Float_t sres = 0.25;
+ Float_t res = sd2+sres;
+ return res;
+}
+
+
+Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){
+ //sigma z2 = in digits - angle estimated supposing vertex constraint
+ Float_t z = iz*fZWidth;
+ Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
+ Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
+ angular*=angular;
+ angular/=12.;
+ Float_t sres = fParam->GetZSigma()/fZWidth;
+ sres *=sres;
+ Float_t res = angular +sd2+sres;
+ return res;
+}
+
+void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t m,
+AliTPCclusterMI &c)
+{
+ Int_t i0=k/max; //central pad
+ Int_t j0=k%max; //central time bin
+
+ // set pointers to data
+ //Int_t dummy[5] ={0,0,0,0,0};
+ Int_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2
+ Int_t * resmatrix[5];
+ for (Int_t di=-2;di<=2;di++){
+ matrix[di+2] = &bins[k+di*max];
+ resmatrix[di+2] = &fResBins[k+di*max];
+ }
+ //build matrix with virtual charge
+ Float_t sigmay2= GetSigmaY2(j0);
+ Float_t sigmaz2= GetSigmaZ2(j0);
+
+ Float_t vmatrix[5][5];
+ vmatrix[2][2] = matrix[2][0];
+ c.SetType(0);
+ c.SetMax(Short_t(vmatrix[2][2])); // write maximal amplitude
+ for (Int_t di =-1;di <=1;di++)
+ for (Int_t dj =-1;dj <=1;dj++){
+ Float_t amp = matrix[di+2][dj];
+ if ( (amp<2) && (fLoop<2)){
+ // if under threshold - calculate virtual charge
+ Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
+ amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
+ if (amp>2) amp = 2;
+ vmatrix[2+di][2+dj]=amp;
+ vmatrix[2+2*di][2+2*dj]=0;
+ if ( (di*dj)!=0){
+ //DIAGONAL ELEMENTS
+ vmatrix[2+2*di][2+dj] =0;
+ vmatrix[2+di][2+2*dj] =0;
+ }
+ continue;
+ }
+ if (amp<4){
+ //if small amplitude - below 2 x threshold - don't consider other one
+ vmatrix[2+di][2+dj]=amp;
+ vmatrix[2+2*di][2+2*dj]=0; // don't take to the account next bin
+ if ( (di*dj)!=0){
+ //DIAGONAL ELEMENTS
+ vmatrix[2+2*di][2+dj] =0;
+ vmatrix[2+di][2+2*dj] =0;
+ }
+ continue;
+ }
+ //if bigger then take everything
+ vmatrix[2+di][2+dj]=amp;
+ vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
+ if ( (di*dj)!=0){
+ //DIAGONAL ELEMENTS
+ vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
+ vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
+ }
+ }
+
+
+
+ Float_t sumw=0;
+ Float_t sumiw=0;
+ Float_t sumi2w=0;
+ Float_t sumjw=0;
+ Float_t sumj2w=0;
+ //
+ for (Int_t i=-2;i<=2;i++)
+ for (Int_t j=-2;j<=2;j++){
+ Float_t amp = vmatrix[i+2][j+2];
+
+ sumw += amp;
+ sumiw += i*amp;
+ sumi2w += i*i*amp;
+ sumjw += j*amp;
+ sumj2w += j*j*amp;
+ }
+ //
+ Float_t meani = sumiw/sumw;
+ Float_t mi2 = sumi2w/sumw-meani*meani;
+ Float_t meanj = sumjw/sumw;
+ Float_t mj2 = sumj2w/sumw-meanj*meanj;
+ //
+ Float_t ry = mi2/sigmay2;
+ Float_t rz = mj2/sigmaz2;
+
+ //
+ if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
+ if ( (ry <1.2) && (rz<1.2) ) {
+ //if cluster looks like expected
+ //+1.2 deviation from expected sigma accepted
+ // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
+
+ meani +=i0;
+ meanj +=j0;
+ //set cluster parameters
+ c.SetQ(sumw);
+ c.SetY(meani*fPadWidth);
+ c.SetZ(meanj*fZWidth);
+ c.SetSigmaY2(mi2);
+ c.SetSigmaZ2(mj2);
+ AddCluster(c);
+ //remove cluster data from data
+ for (Int_t di=-2;di<=2;di++)
+ for (Int_t dj=-2;dj<=2;dj++){
+ resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]);
+ if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
+ }
+ resmatrix[2][0] =0;
+ return;
+ }
+ //
+ //unfolding when neccessary
+ //
+
+ Int_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3
+ Int_t dummy[7]={0,0,0,0,0,0};
+ for (Int_t di=-3;di<=3;di++){
+ matrix2[di+3] = &bins[k+di*max];
+ if ((k+di*max)<3) matrix2[di+3] = &dummy[3];
+ if ((k+di*max)>fMaxBin-3) matrix2[di+3] = &dummy[3];
+ }
+ Float_t vmatrix2[5][5];
+ Float_t sumu;
+ Float_t overlap;
+ UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
+ //
+ // c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
+ meani +=i0;
+ meanj +=j0;
+ //set cluster parameters
+ c.SetQ(sumu);
+ c.SetY(meani*fPadWidth);
+ c.SetZ(meanj*fZWidth);
+ c.SetSigmaY2(mi2);
+ c.SetSigmaZ2(mj2);
+ c.SetType(Char_t(overlap)+1);
+ AddCluster(c);
+
+ //unfolding 2
+ meani-=i0;
+ meanj-=j0;
+ if (gDebug>4)
+ printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
+}
+
+
+
+void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
+ Float_t & sumu, Float_t & overlap )
+{
+ //
+ //unfold cluster from input matrix
+ //data corresponding to cluster writen in recmatrix
+ //output meani and meanj
+
+ //take separatelly y and z
+
+ Float_t sum3i[7] = {0,0,0,0,0,0,0};
+ Float_t sum3j[7] = {0,0,0,0,0,0,0};
+
+ for (Int_t k =0;k<7;k++)
+ for (Int_t l = -1; l<=1;l++){
+ sum3i[k]+=matrix2[k][l];
+ sum3j[k]+=matrix2[l+3][k-3];
+ }
+ Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
+ //
+ //unfold y
+ Float_t sum3wi = 0; //charge minus overlap
+ Float_t sum3wio = 0; //full charge
+ Float_t sum3iw = 0; //sum for mean value
+ for (Int_t dk=-1;dk<=1;dk++){
+ sum3wio+=sum3i[dk+3];
+ if (dk==0){
+ sum3wi+=sum3i[dk+3];
+ }
+ else{
+ Float_t ratio =1;
+ if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
+ sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
+ Float_t xm2 = sum3i[-dk+3];
+ Float_t xm1 = sum3i[+3];
+ Float_t x1 = sum3i[2*dk+3];
+ Float_t x2 = sum3i[3*dk+3];
+ Float_t w11 = TMath::Max(4.*xm1-xm2,0.000001);
+ Float_t w12 = TMath::Max(4 *x1 -x2,0.);
+ ratio = w11/(w11+w12);
+ for (Int_t dl=-1;dl<=1;dl++)
+ mratio[dk+1][dl+1] *= ratio;
+ }
+ Float_t amp = sum3i[dk+3]*ratio;
+ sum3wi+=amp;
+ sum3iw+= dk*amp;
+ }
+ }
+ meani = sum3iw/sum3wi;
+ Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
+
+
+
+ //unfold z
+ Float_t sum3wj = 0; //charge minus overlap
+ Float_t sum3wjo = 0; //full charge
+ Float_t sum3jw = 0; //sum for mean value
+ for (Int_t dk=-1;dk<=1;dk++){
+ sum3wjo+=sum3j[dk+3];
+ if (dk==0){
+ sum3wj+=sum3j[dk+3];
+ }
+ else{
+ Float_t ratio =1;
+ if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
+ (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
+ Float_t xm2 = sum3j[-dk+3];
+ Float_t xm1 = sum3j[+3];
+ Float_t x1 = sum3j[2*dk+3];
+ Float_t x2 = sum3j[3*dk+3];
+ Float_t w11 = TMath::Max(4.*xm1-xm2,0.000001);
+ Float_t w12 = TMath::Max(4 *x1 -x2,0.);
+ ratio = w11/(w11+w12);
+ for (Int_t dl=-1;dl<=1;dl++)
+ mratio[dl+1][dk+1] *= ratio;
+ }
+ Float_t amp = sum3j[dk+3]*ratio;
+ sum3wj+=amp;
+ sum3jw+= dk*amp;
+ }
+ }
+ meanj = sum3jw/sum3wj;
+ Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
+ overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
+ sumu = (sum3wj+sum3wi)/2.;
+
+ if (overlap ==3) {
+ //if not overlap detected remove everything
+ for (Int_t di =-2; di<=2;di++)
+ for (Int_t dj =-2; dj<=2;dj++){
+ recmatrix[di+2][dj+2] = matrix2[3+di][dj];
+ }
+ }
+ else{
+ for (Int_t di =-1; di<=1;di++)
+ for (Int_t dj =-1; dj<=1;dj++){
+ Float_t ratio =1;
+ if (mratio[di+1][dj+1]==1){
+ recmatrix[di+2][dj+2] = matrix2[3+di][dj];
+ if (TMath::Abs(di)+TMath::Abs(dj)>1){
+ recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
+ recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
+ }
+ recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
+ }
+ else
+ {
+ //if we have overlap in direction
+ recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
+ if (TMath::Abs(di)+TMath::Abs(dj)>1){
+ ratio = TMath::Min(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1),1.);
+ recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
+ //
+ ratio = TMath::Min(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1),1.);
+ recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
+ }
+ else{
+ ratio = recmatrix[di+2][dj+2]/matrix2[3][0];
+ recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
+ }
+ }
+ }
+ }
+ if (gDebug>4)
+ printf("%f\n", recmatrix[2][2]);
+
+}
+
+Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
+{
+ //
+ // estimate max
+ Float_t sumteor= 0;
+ Float_t sumamp = 0;
+
+ for (Int_t di = -1;di<=1;di++)
+ for (Int_t dj = -1;dj<=1;dj++){
+ if (vmatrix[2+di][2+dj]>2){
+ Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
+ sumteor += teor*vmatrix[2+di][2+dj];
+ sumamp += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
+ }
+ }
+ Float_t max = sumamp/sumteor;
+ return max;
+}
+
+void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
+ //
+ // transform cluster to the global coordinata
+ // add the cluster to the array
+ //
+ Float_t meani = c.GetY()/fPadWidth;
+ Float_t meanj = c.GetZ()/fZWidth;
+
+ Int_t ki = TMath::Nint(meani-3);
+ if (ki<0) ki=0;
+ if (ki>=fMaxPad) ki = fMaxPad-1;
+ Int_t kj = TMath::Nint(meanj-3);
+ if (kj<0) kj=0;
+ if (kj>=fMaxTime-3) kj=fMaxTime-4;
+ // ki and kj shifted to "real" coordinata
+ c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
+ c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
+ c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
+
+
+ Float_t s2 = c.GetSigmaY2();
+ Float_t w=fParam->GetPadPitchWidth(fSector);
+
+ c.SetSigmaY2(s2*w*w);
+ s2 = c.GetSigmaZ2();
+ w=fZWidth;
+ c.SetSigmaZ2(s2*w*w);
+ c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
+ c.SetZ(fZWidth*(meanj-3));
+ c.SetZ(c.GetZ() - 3.*fParam->GetZSigma()); // PASA delay
+ c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
+
+ if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
+ //c.SetSigmaY2(c.GetSigmaY2()*25.);
+ //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
+ c.SetType(-(c.GetType()+3)); //edge clusters
+ }
+ if (fLoop==2) c.SetType(100);
+
+ TClonesArray * arr = fRowCl->GetArray();
+ AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+
+ fNcluster++;
+}
+
+
+//_____________________________________________________________________________
+void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn)
+{
+ //-----------------------------------------------------------------
+ // This is a simple cluster finder.
+ //-----------------------------------------------------------------
+ TDirectory *savedir=gDirectory;
+
+ if (fInput==0){
+ cerr<<"AliTPC::Digits2Clusters(): input tree not initialised !\n";
+ return;
+ }
+
+ if (fOutput==0) {
+ cerr<<"AliTPC::Digits2Clusters(): output tree not initialised !\n";
+ return;
+ }
+
+ AliSimDigits digarr, *dummy=&digarr;
+ fRowDig = dummy;
+ fInput->GetBranch("Segment")->SetAddress(&dummy);
+ Stat_t nentries = fInput->GetEntries();
+
+ fMaxTime=par->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after
+
+ fParam = par;
+ ((AliTPCParam*)par)->Write(par->GetTitle());
+
+ Int_t nclusters = 0;
+
+ for (Int_t n=0; n<nentries; n++) {
+ fInput->GetEvent(n);
+ Int_t row;
+ if (!par->AdjustSectorRow(digarr.GetID(),fSector,row)) {
+ cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
+ continue;
+ }
+
+ AliTPCClustersRow *clrow= new AliTPCClustersRow();
+ fRowCl = clrow;
+ clrow->SetClass("AliTPCclusterMI");
+ clrow->SetArray(1);
+
+ clrow->SetID(digarr.GetID());
+ fOutput->GetBranch("Segment")->SetAddress(&clrow);
+ fRx=par->GetPadRowRadii(fSector,row);
+
+
+ const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
+ fZWidth = fParam->GetZWidth();
+ if (fSector < kNIS) {
+ fMaxPad = par->GetNPadsLow(row);
+ fSign = (fSector < kNIS/2) ? 1 : -1;
+ fPadLength = par->GetPadPitchLength(fSector,row);
+ fPadWidth = par->GetPadPitchWidth();
+ } else {
+ fMaxPad = par->GetNPadsUp(row);
+ fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
+ fPadLength = par->GetPadPitchLength(fSector,row);
+ fPadWidth = par->GetPadPitchWidth();
+ }
+
+
+ fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after
+ fBins =new Int_t[fMaxBin];
+ fResBins =new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop
+ memset(fBins,0,sizeof(Int_t)*fMaxBin);
+
+ if (digarr.First()) //MI change
+ do {
+ Short_t dig=digarr.CurrentDigit();
+ if (dig<=par->GetZeroSup()) continue;
+ Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
+ fBins[i*fMaxTime+j]=dig;
+ } while (digarr.Next());
+ digarr.ExpandTrackBuffer();
+
+ //add virtual charge at the edge
+ for (Int_t i=0; i<fMaxTime; i++){
+ Float_t amp1 = fBins[i+3*fMaxTime];
+ Float_t amp0 =0;
+ if (amp1>0){
+ Float_t amp2 = fBins[i+4*fMaxTime];
+ if (amp2==0) amp2=0.5;
+ Float_t sigma2 = GetSigmaY2(i);
+ amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
+ if (gDebug>4) printf("\n%f\n",amp0);
+ }
+ fBins[i+2*fMaxTime] = Int_t(amp0);
+ amp0 = 0;
+ amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
+ if (amp1>0){
+ Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
+ if (amp2==0) amp2=0.5;
+ Float_t sigma2 = GetSigmaY2(i);
+ amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
+ if (gDebug>4) printf("\n%f\n",amp0);
+ }
+ fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0);
+ }
+
+ memcpy(fResBins,fBins, fMaxBin*2);
+ //
+ fNcluster=0;
+ //first loop - for "gold cluster"
+ fLoop=1;
+ Int_t *b=&fBins[-1]+2*fMaxTime;
+ Int_t crtime = (fParam->GetZLength()-1.05*fRx)/fZWidth-5;
+
+ for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
+ b++;
+ if (*b<8) continue; //threshold form maxima
+ if (i%fMaxTime<crtime) {
+ Int_t delta = -(i%fMaxTime)+crtime;
+ b+=delta;
+ i+=delta;
+ continue;
+ }
+
+ if (!IsMaximum(*b,fMaxTime,b)) continue;
+ AliTPCclusterMI c;
+ Int_t dummy;
+ MakeCluster(i, fMaxTime, fBins, dummy,c);
+ //}
+ }
+ //memcpy(fBins,fResBins, fMaxBin*2);
+ //second loop - for rest cluster
+ /*
+ fLoop=2;
+ b=&fResBins[-1]+2*fMaxTime;
+ for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
+ b++;
+ if (*b<25) continue; // bigger threshold for maxima
+ if (!IsMaximum(*b,fMaxTime,b)) continue;
+ AliTPCclusterMI c;
+ Int_t dummy;
+ MakeCluster(i, fMaxTime, fResBins, dummy,c);
+ //}
+ }
+ */
+
+ fOutput->Fill();
+ delete clrow;
+ nclusters+=fNcluster;
+ delete[] fBins;
+ delete[] fResBins;
+ }
+ cerr<<"Number of found clusters : "<<nclusters<<" \n";
+
+ fOutput->Write();
+ savedir->cd();
+}
--- /dev/null
+#ifndef ALITPCCLUSTERERMI_H
+#define ALITPCCLUSTERERMI_H
+
+//-------------------------------------------------------
+// TPC clusterer
+//
+// Origin: Marian Ivanov
+//-------------------------------------------------------
+#include <Rtypes.h>
+#include <TObject.h>
+#define kMAXCLUSTER 2500
+
+class TFile;
+class AliTPCParam;
+class AliTPCclusterMI;
+class AliTPCClustersRow;
+class AliSimDigits;
+class TTree;
+
+class AliTPCclustererMI : public TObject{
+public:
+ AliTPCclustererMI();
+ virtual void Digits2Clusters(const AliTPCParam *par, Int_t eventn=1);
+ virtual void SetInput(TTree * tree); // set input tree with digits
+ virtual void SetOutput(TTree * tree); //set output tree with
+private:
+ Bool_t IsMaximum(Int_t k, Int_t max, const Int_t *bins);
+ void MakeCluster2(Int_t k,Int_t max,Int_t *bins,UInt_t m,
+ AliTPCclusterMI &c);
+ void MakeCluster(Int_t k,Int_t max,Int_t *bins,UInt_t m,
+ AliTPCclusterMI &c);
+ Float_t GetSigmaY2(Int_t iz);
+ Float_t GetSigmaZ2(Int_t iz);
+ Float_t FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz);
+ void AddCluster(AliTPCclusterMI &c); // add the cluster to the array
+ void UnfoldCluster(Int_t * matrix[7], Float_t recmatrix[5][5],
+ Float_t & meani, Float_t & meanj, Float_t & sum, Float_t &overlap );
+
+
+
+ Int_t * fBins; //!digits array
+ Int_t * fResBins; //!digits array with res. after 1 finder
+ Int_t fLoop; //loop - cf in 2 loops
+ Int_t fMaxBin;
+ Int_t fMaxTime;
+ Int_t fMaxPad;
+ Int_t fSector; //!current sector
+ Float_t fSign; //!current sign
+ Float_t fRx; // current radius
+ Float_t fPadWidth; // the width of the pad
+ Float_t fPadLength; // the width of the pad
+ Float_t fZWidth; //the z bin width
+
+ TTree * fInput; //!input tree with digits - object not owner
+ TTree * fOutput; //!output tree with digits - object not owner
+ AliTPCClustersRow * fRowCl; //! current cluster row
+ AliSimDigits * fRowDig; //! current digits row
+ const AliTPCParam * fParam; //! tpc parameters
+ Int_t fNcluster; // number of clusters - for given row
+ ClassDef(AliTPCclustererMI,1) // Time Projection Chamber digits
+};
+
+inline Bool_t AliTPCclustererMI::IsMaximum(Int_t q,Int_t max,const Int_t *bins){
+ //is this a local maximum ?
+ if (bins[-max] >= q) return kFALSE;
+ if (bins[-1 ] >= q) return kFALSE;
+ if (bins[+max] > q) return kFALSE;
+ if (bins[+1 ] > q) return kFALSE;
+ if (bins[-max-1] >= q) return kFALSE;
+ if (bins[+max-1] >= q) return kFALSE;
+ if (bins[+max+1] > q) return kFALSE;
+ if (bins[-max+1] >= q) return kFALSE;
+ return kTRUE;
+}
+
+
+
+//-----------------------------------------------------------------
+
+#endif
+
+
Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i);
void ResetCovariance();
-private:
+protected:
Double_t fX; // X-coordinate of this track (reference plane)
Double_t fAlpha; // Rotation angle the local (TPC sector)
// coordinate system and the global ALICE one.
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+
+
+
+
+/*
+ AliTPC parallel tracker -
+ How to use? -
+ run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
+ run AliTPCFindTracksMI.C macro - to find tracks
+ tracks are written to AliTPCtracks.root file
+ for comparison also seeds are written to the same file - to special branch
+*/
+
+//-------------------------------------------------------
+// Implementation of the TPC tracker
+//
+// Origin: Marian Ivanov Marian.Ivanov@cern.ch
+//
+//-------------------------------------------------------
+
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <iostream.h>
+
+#include "AliTPCtrackerMI.h"
+#include "AliTPCclusterMI.h"
+#include "AliTPCParam.h"
+#include "AliTPCClustersRow.h"
+#include "AliComplexCluster.h"
+#include "TStopwatch.h"
+
+
+ClassImp(AliTPCseed)
+
+AliTPCclusterTracks::AliTPCclusterTracks(){
+ // class for storing overlaping info
+ fTrackIndex[0]=-1;
+ fTrackIndex[1]=-1;
+ fTrackIndex[2]=-1;
+ fDistance[0]=1000;
+ fDistance[1]=1000;
+ fDistance[2]=1000;
+}
+
+
+
+
+
+Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
+
+ Int_t sec=(i&0xff000000)>>24;
+ Int_t row = (i&0x00ff0000)>>16;
+ track->fRow=(i&0x00ff0000)>>16;
+ track->fSector = sec;
+ // Int_t index = i&0xFFFF;
+ if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
+ track->fClusterIndex[track->fRow] = i;
+ track->fFirstPoint = row;
+ //
+
+ AliTPCTrackPoint *trpoint =track->GetTrackPoint(track->fRow);
+ Float_t angle2 = track->GetSnp()*track->GetSnp();
+ angle2 = TMath::Sqrt(angle2/(1-angle2));
+ //
+ //SET NEW Track Point
+ //
+ if (c!=0){
+ //if we have a cluster
+ trpoint->GetCPoint().SetY(c->GetY());
+ trpoint->GetCPoint().SetZ(c->GetZ());
+ //
+ trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
+ trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
+ //
+ trpoint->GetCPoint().SetType(c->GetType());
+ trpoint->GetCPoint().SetQ(c->GetQ());
+ trpoint->GetCPoint().SetMax(c->GetMax());
+ //
+ trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
+ trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
+ //
+ }
+ trpoint->GetTPoint().SetX(track->GetX());
+ trpoint->GetTPoint().SetY(track->GetY());
+ trpoint->GetTPoint().SetZ(track->GetZ());
+ //
+ trpoint->GetTPoint().SetAngleY(angle2);
+ trpoint->GetTPoint().SetAngleZ(track->GetTgl());
+
+
+
+ if (chi2>10){
+ // printf("suspicious chi2 %f\n",chi2);
+ }
+ // if (track->fIsSeeding){
+ track->fErrorY2 *= 1.2;
+ track->fErrorY2 += 0.0064;
+ track->fErrorZ2 *= 1.2;
+ track->fErrorY2 += 0.005;
+
+ //}
+
+ return track->Update(c,chi2,i);
+
+}
+//_____________________________________________________________________________
+AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn):
+AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
+{
+ //---------------------------------------------------------------------
+ // The main TPC tracker constructor
+ //---------------------------------------------------------------------
+ fInnerSec=new AliTPCSector[fkNIS];
+ fOuterSec=new AliTPCSector[fkNOS];
+
+ Int_t i;
+ for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
+ for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
+
+ fN=0; fSectors=0;
+
+ fClustersArray.Setup(par);
+ fClustersArray.SetClusterType("AliTPCclusterMI");
+
+ char cname[100];
+ if (eventn==-1) {
+ sprintf(cname,"TreeC_TPC");
+ }
+ else {
+ sprintf(cname,"TreeC_TPC_%d",eventn);
+ }
+
+ fClustersArray.ConnectTree(cname);
+
+ fEventN = eventn;
+ fSeeds=0;
+ fNtracks = 0;
+ fParam = par;
+}
+
+//_____________________________________________________________________________
+AliTPCtrackerMI::~AliTPCtrackerMI() {
+ //------------------------------------------------------------------
+ // TPC tracker destructor
+ //------------------------------------------------------------------
+ delete[] fInnerSec;
+ delete[] fOuterSec;
+ if (fSeeds) {
+ fSeeds->Delete();
+ delete fSeeds;
+ }
+}
+
+
+Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
+ //
+ //
+ Float_t snoise2;
+ Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
+
+ //cluster "quality"
+ Float_t rsigmay = 1;
+ Int_t ctype = 0;
+
+ //standard if we don't have cluster - take MIP
+ const Float_t chmip = 50.;
+ Float_t amp = chmip/0.3;
+ Float_t nel;
+ Float_t nprim;
+ if (cl){
+ amp = cl->GetQ();
+ rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
+ ctype = cl->GetType();
+ }
+
+
+ Float_t landau=2 ; //landau fluctuation part
+ Float_t gg=2; // gg fluctuation part
+ Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
+
+
+ if (fSectors==fInnerSec){
+ snoise2 = 0.0004;
+ nel = 0.268*amp;
+ nprim = 0.155*amp;
+ gg = (2+0.0002*amp)/nel;
+ landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+ if (landau>1) landau=1;
+ }
+ else {
+ snoise2 = 0.0004;
+ nel = 0.3*amp;
+ nprim = 0.133*amp;
+ gg = (2+0.0002*amp)/nel;
+ landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+ if (landau>1) landau=1;
+ }
+
+
+ Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
+ Float_t angle2 = seed->GetSnp()*seed->GetSnp();
+ angle2 = angle2/(1-angle2);
+ Float_t angular = landau*angle2*padlength*padlength/12.;
+ Float_t res = sdiff + angular;
+
+
+ if ((ctype==0) && (fSectors ==fOuterSec))
+ res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
+
+ if ((ctype==0) && (fSectors ==fInnerSec))
+ res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
+
+
+ if ((ctype>0))
+ res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
+
+ if (ctype<0)
+ res*=2.4; // overestimate error 2 times
+
+ res+= snoise2;
+
+ if (res<2*snoise2)
+ res = 2*snoise2;
+
+ seed->SetErrorY2(res);
+ return res;
+}
+
+
+
+
+Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
+ //
+ //
+ Float_t snoise2;
+ Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
+ //signal quality
+ Float_t rsigmaz=1;
+ Int_t ctype =0;
+
+ const Float_t chmip = 50.;
+ Float_t amp = chmip/0.3;
+ Float_t nel;
+ Float_t nprim;
+ if (cl){
+ amp = cl->GetQ();
+ rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
+ ctype = cl->GetType();
+ }
+
+ //
+ Float_t landau=2 ; //landau fluctuation part
+ Float_t gg=2; // gg fluctuation part
+ Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
+
+ if (fSectors==fInnerSec){
+ snoise2 = 0.0004;
+ nel = 0.268*amp;
+ nprim = 0.155*amp;
+ gg = (2+0.0002*amp)/nel;
+ landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+ if (landau>1) landau=1;
+ }
+ else {
+ snoise2 = 0.0004;
+ nel = 0.3*amp;
+ nprim = 0.133*amp;
+ gg = (2+0.0002*amp)/nel;
+ landau = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+ if (landau>1) landau=1;
+ }
+ Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
+
+ Float_t angle = seed->GetTgl();
+ Float_t angular = landau*angle*angle*padlength*padlength/12.;
+ Float_t res = sdiff + angular;
+
+ if ((ctype==0) && (fSectors ==fOuterSec))
+ res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
+
+ if ((ctype==0) && (fSectors ==fInnerSec))
+ res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
+ if ((ctype>0))
+ res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064; //0.31+0.147*ctype;
+ if (ctype<0)
+ res*=1.3;
+ if ((ctype<0) &&<70)
+ res*=1.3;
+
+ res += snoise2;
+ if (res<2*snoise2)
+ res = 2*snoise2;
+
+ seed->SetErrorZ2(res);
+ return res;
+}
+
+
+
+Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
+{
+ //-----------------------------------------------------------------
+ // This function find proloncation of a track to a reference plane x=xk.
+ // doesn't change internal state of the track
+ //-----------------------------------------------------------------
+
+ Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
+ // Double_t y1=fP0, z1=fP1;
+ Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
+ Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
+
+ y = fP0;
+ z = fP1;
+ y += dx*(c1+c2)/(r1+r2);
+ z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
+ return 0;
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
+{
+ //-----------------------------------------------------------------
+ // This function calculates a predicted chi2 increment.
+ //-----------------------------------------------------------------
+ //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+ Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+ r00+=fC00; r01+=fC10; r11+=fC11;
+
+ Double_t det=r00*r11 - r01*r01;
+ if (TMath::Abs(det) < 1.e-10) {
+ Int_t n=GetNumberOfClusters();
+ if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+ return 1e10;
+ }
+ Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+
+ Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+
+ return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
+}
+
+
+//_________________________________________________________________________________________
+
+
+Int_t AliTPCseed::Compare(const TObject *o) const {
+ //-----------------------------------------------------------------
+ // This function compares tracks according to the sector - for given sector according z
+ //-----------------------------------------------------------------
+ AliTPCseed *t=(AliTPCseed*)o;
+ if (t->fSector>fSector) return -1;
+ if (t->fSector<fSector) return 1;
+
+ Double_t z2 = t->GetZ();
+ Double_t z1 = GetZ();
+ if (z2>z1) return 1;
+ if (z2<z1) return -1;
+ return 0;
+}
+
+
+
+
+
+
+//_____________________________________________________________________________
+Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
+ //-----------------------------------------------------------------
+ // This function associates a cluster with this track.
+ //-----------------------------------------------------------------
+ // Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+ //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
+ Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+
+ r00+=fC00; r01+=fC10; r11+=fC11;
+ Double_t det=r00*r11 - r01*r01;
+ Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+ Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
+ Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
+ Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
+ Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
+ Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
+
+ Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+ Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
+ if (TMath::Abs(cur*fX-eta) >= 0.99999) {
+ Int_t n=GetNumberOfClusters();
+ if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
+ return 0;
+ }
+
+ fP0 += k00*dy + k01*dz;
+ fP1 += k10*dy + k11*dz;
+ fP2 = eta;
+ fP3 += k30*dy + k31*dz;
+ fP4 = cur;
+
+ Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
+ Double_t c12=fC21, c13=fC31, c14=fC41;
+
+ fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
+ fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
+ fC40-=k00*c04+k01*c14;
+
+ fC11-=k10*c01+k11*fC11;
+ fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
+ fC41-=k10*c04+k11*c14;
+
+ fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
+ fC42-=k20*c04+k21*c14;
+
+ fC33-=k30*c03+k31*c13;
+ fC43-=k40*c03+k41*c13;
+
+ fC44-=k40*c04+k41*c14;
+
+ Int_t n=GetNumberOfClusters();
+ fIndex[n]=index;
+ SetNumberOfClusters(n+1);
+ SetChi2(GetChi2()+chisq);
+
+ return 1;
+}
+
+
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
+ Double_t x2,Double_t y2,
+ Double_t x3,Double_t y3)
+{
+ //-----------------------------------------------------------------
+ // Initial approximation of the track curvature
+ //-----------------------------------------------------------------
+ Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+ Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+ (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+ Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+ (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+ Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+ return -xr*yr/sqrt(xr*xr+yr*yr);
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
+ Double_t x2,Double_t y2,
+ Double_t x3,Double_t y3)
+{
+ //-----------------------------------------------------------------
+ // Initial approximation of the track curvature times center of curvature
+ //-----------------------------------------------------------------
+ Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+ Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+ (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+ Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+ (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+ Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+ return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
+}
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
+ Double_t x2,Double_t y2,
+ Double_t z1,Double_t z2)
+{
+ //-----------------------------------------------------------------
+ // Initial approximation of the tangent of the track dip angle
+ //-----------------------------------------------------------------
+ return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
+}
+
+
+void AliTPCtrackerMI::LoadClusters()
+{
+ //
+ // load clusters to the memory
+ Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+ for (Int_t i=0; i<j; i++) {
+ fClustersArray.LoadEntry(i);
+ }
+}
+
+void AliTPCtrackerMI::UnloadClusters()
+{
+ //
+ // load clusters to the memory
+ Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+ for (Int_t i=0; i<j; i++) {
+ fClustersArray.ClearSegment(i);
+ }
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::LoadOuterSectors() {
+ //-----------------------------------------------------------------
+ // This function fills outer TPC sectors with clusters.
+ //-----------------------------------------------------------------
+ UInt_t index;
+ Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+ for (Int_t i=0; i<j; i++) {
+ // AliSegmentID *s=fClustersArray.LoadEntry(i);
+ AliSegmentID *s= fClustersArray.At(i);
+ Int_t sec,row;
+ AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
+ par->AdjustSectorRow(s->GetID(),sec,row);
+ if (sec<fkNIS*2) continue;
+ AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
+ Int_t ncl=clrow->GetArray()->GetEntriesFast();
+ while (ncl--) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
+ index=(((sec<<8)+row)<<16)+ncl;
+ fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
+ }
+ }
+ fN=fkNOS;
+ fSectors=fOuterSec;
+}
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::LoadInnerSectors() {
+ //-----------------------------------------------------------------
+ // This function fills inner TPC sectors with clusters.
+ //-----------------------------------------------------------------
+ UInt_t index;
+ Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+ for (Int_t i=0; i<j; i++) {
+ // AliSegmentID *s=fClustersArray.LoadEntry(i);
+ AliSegmentID *s=fClustersArray.At(i);
+ Int_t sec,row;
+ AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
+ par->AdjustSectorRow(s->GetID(),sec,row);
+ if (sec>=fkNIS*2) continue;
+ AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
+ Int_t ncl=clrow->GetArray()->GetEntriesFast();
+ while (ncl--) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
+ index=(((sec<<8)+row)<<16)+ncl;
+ fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
+ }
+ }
+
+ fN=fkNIS;
+ fSectors=fInnerSec;
+}
+
+Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
+ //-----------------------------------------------------------------
+ // This function tries to find a track prolongation to next pad row
+ //-----------------------------------------------------------------
+ // Double_t xt=t.GetX();
+ // Int_t row = fSectors->GetRowNumber(xt)-1;
+ // if (row < nr) return 1; // don't prolongate if not information until now -
+ //
+ Double_t x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
+ if (!t.PropagateTo(x)) return 0;
+ // update current
+ t.fCurrentSigmaY = GetSigmaY(&t);
+ t.fCurrentSigmaZ = GetSigmaZ(&t);
+ //
+ AliTPCclusterMI *cl=0;
+ UInt_t index=0;
+ const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+ Double_t sy2=ErrY2(&t)*2;
+ Double_t sz2=ErrZ2(&t)*2;
+
+
+ Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
+ Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
+ Double_t y=t.GetY(), z=t.GetZ();
+
+ if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+ t.fInDead = kTRUE;
+ return 0;
+ }
+ else
+ {
+ t.fNFoundable++;
+ }
+ //calculate
+ Float_t maxdistance = roady*roady + roadz*roadz;
+ if (krow) {
+ for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
+ if (c->GetZ() > z+roadz) break;
+ if ( (c->GetY()-y) > roady ) continue;
+ Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
+ if (maxdistance>distance) {
+ maxdistance = distance;
+ cl=c;
+ index=krow.GetIndex(i);
+ }
+ }
+ }
+ if (cl) {
+ //Double_t sy2=
+ ErrY2(&t,cl);
+ //Double_t sz2=
+ ErrZ2(&t,cl);
+ Double_t chi2= t.GetPredictedChi2(cl);
+ UpdateTrack(&t,cl,chi2,index);
+ } else {
+ if (y > ymax) {
+ t.fRelativeSector= (t.fRelativeSector+1) % fN;
+ if (!t.Rotate(fSectors->GetAlpha()))
+ return 0;
+ } else if (y <-ymax) {
+ t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
+ if (!t.Rotate(-fSectors->GetAlpha()))
+ return 0;
+ }
+ }
+ return 1;
+}
+
+
+Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex, Int_t nr) {
+ //-----------------------------------------------------------------
+ // This function tries to find a track prolongation to next pad row
+ //-----------------------------------------------------------------
+ t.fCurrentCluster = 0;
+ t.fCurrentClusterIndex1 = 0;
+ t.fCurrentClusterIndex2 = 0;
+
+ Double_t xt=t.GetX();
+ Int_t row = fSectors->GetRowNumber(xt)-1;
+ if (row < nr) return 1; // don't prolongate if not information until now -
+ Double_t x=fSectors->GetX(nr);
+ if (!t.PropagateTo(x)) return 0;
+ // update current
+ t.fCurrentSigmaY = GetSigmaY(&t);
+ t.fCurrentSigmaZ = GetSigmaZ(&t);
+
+ AliTPCclusterMI *cl=0;
+ UInt_t index=0;
+ AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+ //
+ Double_t y=t.GetY(), z=t.GetZ();
+ Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
+ Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
+ //
+
+ Float_t maxdistance = 1000000;
+ if (krow) {
+ for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
+ if (c->GetZ() > z+roadz) break;
+ if (TMath::Abs(c->GetY()-y)>roady) continue;
+
+ //krow.UpdateClusterTrack(i,trindex,&t);
+
+ Float_t dy2 = (c->GetY()- t.GetY());
+ dy2*=dy2;
+ Float_t dz2 = (c->GetZ()- t.GetZ());
+ dz2*=dz2;
+ //
+ Float_t distance = dy2+dz2;
+ //
+ if (distance > maxdistance) continue;
+ maxdistance = distance;
+ cl=c;
+ index=i;
+ }
+ }
+ t.fCurrentCluster = cl;
+ t.fCurrentClusterIndex1 = krow.GetIndex(index);
+ t.fCurrentClusterIndex2 = index;
+ return 1;
+}
+
+
+Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
+ //-----------------------------------------------------------------
+ // This function tries to find a track prolongation to next pad row
+ //-----------------------------------------------------------------
+ AliTPCseed & t = *((AliTPCseed*)(fSeeds->At(trindex)));
+ AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+ // Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
+ Double_t y=t.GetY();
+ Double_t ymax=fSectors->GetMaxY(nr);
+
+ if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+ t.fInDead = kTRUE;
+ return 0;
+ }
+ else
+ {
+ t.fNFoundable++;
+ }
+
+ if (t.fCurrentCluster) {
+ // Float_t l=fSectors->GetPadPitchLength();
+ // AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
+
+ Double_t sy2=ErrY2(&t,t.fCurrentCluster);
+ Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
+
+
+ Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
+ Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
+
+ Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
+ Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
+
+ Double_t rdistance = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
+
+
+ // printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
+ if (rdistance>4) return 0;
+
+ if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
+ return 0; //suspisiouce - will be changed
+
+ if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
+ // strict cut on overlaped cluster
+ return 0; //suspisiouce - will be changed
+
+ if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
+ && t.fCurrentCluster->GetType()<0)
+ return 0;
+
+ // t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
+ UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
+
+ } else {
+ if (y > ymax) {
+ t.fRelativeSector= (t.fRelativeSector+1) % fN;
+ if (!t.Rotate(fSectors->GetAlpha()))
+ return 0;
+ } else if (y <-ymax) {
+ t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
+ if (!t.Rotate(-fSectors->GetAlpha()))
+ return 0;
+ }
+ }
+ return 1;
+}
+
+
+
+
+/*
+Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
+{
+ //-----------------------------------------------------------------
+ // fast prolongation mathod -
+ // don't update track only after step clusters
+ //-----------------------------------------------------------------
+ Double_t xt=t.GetX();
+ //
+ Double_t alpha=t.GetAlpha();
+ alpha =- fSectors->GetAlphaShift();
+ if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+ if (alpha < 0. ) alpha += 2.*TMath::Pi();
+ t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+ Int_t row0 = fSectors->GetRowNumber(xt);
+ Double_t x = fSectors->GetX(row0);
+ Double_t ymax = fSectors->GetMaxY(row0);
+ //
+ Double_t sy2=ErrY2(&t)*2;
+ Double_t sz2=ErrZ2(&t)*2;
+ Double_t roady =3.*sqrt(t.GetSigmaY2() + sy2);
+ Double_t roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
+ Float_t maxdistance = roady*roady + roadz*roadz;
+ t.fCurrentSigmaY = GetSigmaY(&t);
+ t.fCurrentSigmaZ = GetSigmaZ(&t);
+ //
+ Int_t nclusters = 0;
+ Double_t y;
+ Double_t z;
+ Double_t yy[200]; //track prolongation
+ Double_t zz[200];
+ Double_t cy[200]; // founded cluster position
+ Double_t cz[200];
+ Double_t sy[200]; // founded cluster error
+ Double_t sz[200];
+ Bool_t hitted[200]; // indication of cluster presence
+ //
+
+ //
+ for (Int_t drow = step; drow>=0; drow--) {
+ Int_t row = row0-drow;
+ if (row<0) break;
+ Double_t x = fSectors->GetX(row);
+ Double_t ymax = fSectors->GetMaxY(row);
+ t.GetProlongation(x,y,z);
+ yy[drow] =y;
+ zz[drow] =z;
+ const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
+ if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+ t.fInDead = kTRUE;
+ break;
+ }
+ else
+ {
+ t.fNFoundable++;
+ }
+
+ //find nearest cluster
+ AliTPCclusterMI *cl= 0;
+ if (krow) {
+ for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+ AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
+ if (c->GetZ() > z+roadz) break;
+ if ( (c->GetY()-y) > roady ) continue;
+ Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
+ if (maxdistance>distance) {
+ maxdistance = distance;
+ cl=c;
+ // index=krow.GetIndex(i);
+ }
+ }
+ } // end of seearch
+ //update cluster information
+ if (cl){
+ cy[drow] = cl->GetY();
+ cz[drow] = cl->GetZ();
+ sy[drow] = ErrY2(&t,cl);
+ sz[drow] = ErrZ2(&t,cl);
+ hitted[drow] = kTRUE;
+ nclusters++;
+ }
+ else
+ hitted[drow] = kFALSE;
+ }
+ //if we have information - update track
+ if (nclusters>0){
+ Float_t sumyw0 = 0;
+ Float_t sumzw0 = 0;
+ Float_t sumyw = 0;
+ Float_t sumzw = 0;
+ Float_t sumrow = 0;
+ for (Int_t i=0;i<step;i++)
+ if (hitted[i]){
+ sumyw0+= 1/sy[i];
+ sumzw0+= 1/sz[i];
+ //
+ sumyw+= (cy[i]-yy[i])/sy[i];
+ sumzw+= (cz[i]-zz[i])/sz[i];
+ sumrow+= i;
+ }
+ Float_t dy = sumyw/sumyw0;
+ Float_t dz = sumzw/sumzw0;
+ Float_t mrow = sumrow/nclusters+row0;
+ Float_t x = fSectors->GetX(mrow);
+ t.PropagateTo(x);
+ AliTPCclusterMI cvirtual;
+ cvirtual.SetZ(dz+t.GetZ());
+ cvirtual.SetY(dy+t.GetY());
+ t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
+ t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
+ Float_t chi2 = t.GetPredictedChi2(&cvirtual);
+ t.Update(&cvirtual,chi2,0);
+ Int_t ncl = t.GetNumberOfClusters();
+ ncl = ncl-1+nclusters;
+ t.SetN(ncl);
+ }
+ return 1;
+}
+*/
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
+ //-----------------------------------------------------------------
+ // This function tries to find a track prolongation.
+ //-----------------------------------------------------------------
+ Double_t xt=t.GetX();
+ //
+ Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+ if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+ if (alpha < 0. ) alpha += 2.*TMath::Pi();
+ t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+
+ for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
+
+ if (FollowToNext(t,nr)==0) {
+ }
+ }
+ return 1;
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
+ //-----------------------------------------------------------------
+ // This function tries to find a track prolongation.
+ //-----------------------------------------------------------------
+ Double_t xt=t.GetX();
+ //
+ Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+ if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+ if (alpha < 0. ) alpha += 2.*TMath::Pi();
+ t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+
+ for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
+ if (t.GetSnp()>0.2)
+ FollowToNext(t,nr);
+ }
+ return 1;
+}
+
+
+
+
+
+Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
+{
+ //
+ //
+ sum1=0;
+ sum2=0;
+ Int_t sum=0;
+ if (s1->fSector!=s2->fSector) return 0;
+ //
+ Float_t dz2 =(s1->GetZ() - s2->GetZ());
+ dz2*=dz2;
+ Float_t dy2 =(s1->GetY() - s2->GetY());
+ dy2*=dy2;
+ Float_t distance = TMath::Sqrt(dz2+dy2);
+ if (distance>5.) return 0; // if there are far away - not overlap - to reduce combinatorics
+
+ for (Int_t i=0;i<160;i++){
+ if (s1->fClusterIndex[i]>0) sum1++;
+ if (s2->fClusterIndex[i]>0) sum2++;
+ if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
+ sum++;
+ }
+ }
+
+ Float_t summin = TMath::Min(sum1,sum2);
+ Float_t ratio = sum/Float_t(summin);
+ return ratio;
+}
+
+void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
+{
+ //
+ //
+ if (s1->fSector!=s2->fSector) return;
+ //
+ Float_t dz2 =(s1->GetZ() - s2->GetZ());
+ dz2*=dz2;
+ Float_t dy2 =(s1->GetY() - s2->GetY());
+ dy2*=dy2;
+ Float_t distance = TMath::Sqrt(dz2+dy2);
+ if (distance>15.) return ; // if there are far away - not overlap - to reduce combinatorics
+ //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
+ // TClonesArray &pointarray1 = *(s1->fPoints);
+ //TClonesArray &pointarray2 = *(s2->fPoints);
+ //
+ for (Int_t i=0;i<160;i++){
+ if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
+ // AliTPCTrackPoint *p1 = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
+ //AliTPCTrackPoint *p2 = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i));
+ AliTPCTrackPoint *p1 = s1->GetTrackPoint(i);
+ AliTPCTrackPoint *p2 = s2->GetTrackPoint(i);;
+ p1->fIsShared = kTRUE;
+ p2->fIsShared = kTRUE;
+ }
+ }
+}
+
+
+
+
+void AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared=kFALSE){
+
+
+
+ //
+ // remove overlap - used removal factor - removal index stored in the track
+ arr->Sort(); // sorting according z
+ arr->Expand(arr->GetEntries());
+ Int_t nseed=arr->GetEntriesFast();
+ // printf("seeds \t%p \t%d\n",arr, nseed);
+ // arr->Expand(arr->GetEntries()); //remove 0 pointers
+ nseed = arr->GetEntriesFast();
+ Int_t removed = 0;
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
+ if (!pt) {
+ continue;
+ }
+ if (!(pt->IsActive())) continue;
+ for (Int_t j=i+1; j<nseed; j++){
+ AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
+ if ((pt2) && pt2->IsActive())
+ if (pt->fSector == pt2->fSector)
+ if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
+ Int_t sum1,sum2;
+ Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
+ if (ratio>factor){
+ // if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex); // arr->RemoveAt(j);
+ Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
+ Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
+ removed++;
+ if (TMath::Abs(ratio3)>0.025){ // if much more points
+ if (sum1>sum2) pt2->Desactivate(removalindex);
+ else {
+ pt->Desactivate(removalindex); // arr->RemoveAt(i);
+ break;
+ }
+ }
+ else{ //decide on mean chi2
+ if (ratio2<1)
+ pt2->Desactivate(removalindex);
+ else {
+ pt->Desactivate(removalindex); // arr->RemoveAt(i);
+ break;
+ }
+ }
+
+ } // if suspicious ratio
+ }
+ else
+ break;
+ }
+ }
+ // printf("removed\t%d\n",removed);
+ Int_t good =0;
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
+ if (!pt) break;
+ if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
+ //desactivate tracks with small number of points
+ // printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
+ pt->Desactivate(10); //desactivate - small muber of points
+ }
+ if (!(pt->IsActive())) continue;
+ good++;
+ }
+
+
+ if (shared)
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
+ if (!pt) continue;
+ if (!(pt->IsActive())) continue;
+ for (Int_t j=i+1; j<nseed; j++){
+ AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
+ if ((pt2) && pt2->IsActive()) {
+ if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
+ SignShared(pt,pt2);
+ }
+ }
+ }
+ fNtracks = good;
+ printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
+
+
+}
+void AliTPCtrackerMI::MakeSeedsAll()
+{
+ if (fSeeds == 0) fSeeds = new TObjArray;
+ TObjArray * arr;
+ for (Int_t sec=0;sec<fkNOS;sec+=3){
+ arr = MakeSeedsSectors(sec,sec+3);
+ Int_t nseed = arr->GetEntriesFast();
+ for (Int_t i=0;i<nseed;i++)
+ fSeeds->AddLast(arr->RemoveAt(i));
+ }
+ // fSeeds = MakeSeedsSectors(0,fkNOS);
+}
+
+TObjArray * AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
+{
+ //
+ // loop over all sectors and make seed
+ //find track seeds
+ TStopwatch timer;
+ Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
+ Int_t nrows=nlow+nup;
+ Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
+ // if (fSeeds==0) fSeeds = new TObjArray;
+ TObjArray * arr = new TObjArray;
+
+ for (Int_t sec=sec1; sec<sec2;sec++){
+ MakeSeeds(arr, sec, nup-1, nup-1-gap);
+ MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
+ }
+ Int_t nseed=arr->GetEntriesFast();
+ gap=Int_t(0.3*nrows);
+ // continue seeds
+ Int_t i;
+ for (i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
+ if (!pt) continue;
+ if (FollowProlongation(t,nup-gap)) {
+ pt->fIsSeeding =kFALSE;
+ continue;
+ }
+ delete arr->RemoveAt(i);
+ }
+ //
+ //remove seeds which overlaps
+ RemoveOverlap(arr,0.6,1);
+ //delete seeds - which were sign
+ nseed=arr->GetEntriesFast();
+ for (i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
+ //, &t=*pt;
+ if (!pt) continue;
+ if ((pt->IsActive()) && pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
+ //FollowBackProlongation(t,nup);
+ continue;
+ }
+ delete arr->RemoveAt(i);
+ }
+ return arr;
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
+ //-----------------------------------------------------------------
+ // This function creates track seeds.
+ //-----------------------------------------------------------------
+ // if (fSeeds==0) fSeeds=new TObjArray(15000);
+
+ Double_t x[5], c[15];
+
+ Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
+ Double_t cs=cos(alpha), sn=sin(alpha);
+
+ Double_t x1 =fOuterSec->GetX(i1);
+ Double_t xx2=fOuterSec->GetX(i2);
+ //
+ // for (Int_t ns=0; ns<fkNOS; ns++)
+ Int_t ns =sec;
+ {
+ Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
+ Int_t nm=fOuterSec[ns][i2];
+ Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
+ const AliTPCRow& kr1=fOuterSec[ns][i1];
+ AliTPCRow& kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
+ AliTPCRow& kr22 = fOuterSec[(ns)%fkNOS][i2];
+ AliTPCRow& kr23 = fOuterSec[(ns+1)%fkNOS][i2];
+
+ for (Int_t is=0; is < kr1; is++) {
+ Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
+ Double_t x3=GetX(), y3=GetY(), z3=GetZ();
+
+ Float_t anglez = (z1-z3)/(x1-x3);
+ Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
+
+ for (Int_t js=0; js < nl+nm+nu; js++) {
+ const AliTPCclusterMI *kcl;
+ Double_t x2, y2, z2;
+ if (js<nl) {
+ if (js==0) {
+ js = kr21.Find(extraz-15.);
+ if (js>=nl) continue;
+ }
+ kcl=kr21[js];
+ z2=kcl->GetZ();
+ if ((extraz-z2)>10) continue;
+ if ((extraz-z2)<-10) {
+ js = nl;
+ continue;
+ }
+ y2=kcl->GetY();
+ x2= xx2*cs+y2*sn;
+ y2=-xx2*sn+y2*cs;
+ } else
+ if (js<nl+nm) {
+ if (js==nl) {
+ js = nl+kr22.Find(extraz-15.);
+ if (js>=nl+nm) continue;
+ }
+ kcl=kr22[js-nl];
+ z2=kcl->GetZ();
+ if ((extraz-z2)>10) continue;
+ if ((extraz-z2)<-10) {
+ js = nl+nm;
+ continue;
+ }
+ x2=xx2; y2=kcl->GetY();
+ } else {
+ //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
+ if (js==nl+nm) {
+ js = nl+nm+kr23.Find(extraz-15.);
+ if (js>=nl+nm+nu) break;
+ }
+ kcl=kr23[js-nl-nm];
+ z2=kcl->GetZ();
+ if ((extraz-z2)>10) continue;
+ if ((extraz-z2)<-10) {
+ break;
+ }
+ y2=kcl->GetY();
+ x2=xx2*cs-y2*sn;
+ y2=xx2*sn+y2*cs;
+ }
+
+ Double_t zz=z1 - anglez*(x1-x2);
+ if (TMath::Abs(zz-z2)>10.) continue;
+
+ Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
+ if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
+
+ x[0]=y1;
+ x[1]=z1;
+ x[4]=f1(x1,y1,x2,y2,x3,y3);
+ if (TMath::Abs(x[4]) >= 0.0066) continue;
+ x[2]=f2(x1,y1,x2,y2,x3,y3);
+ //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
+ x[3]=f3(x1,y1,x2,y2,z1,z2);
+ if (TMath::Abs(x[3]) > 1.2) continue;
+ Double_t a=asin(x[2]);
+ Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
+ if (TMath::Abs(zv-z3)>10.) continue;
+
+ Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
+ Double_t sy2=kcl->GetSigmaY2()*2, sz2=kcl->GetSigmaZ2()*4;
+ //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
+ Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
+ //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
+
+ Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
+ Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
+ Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
+ Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
+ Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
+ Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
+ Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
+ Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
+ Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
+ Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
+
+ c[0]=sy1;
+ c[1]=0.; c[2]=sz1;
+ c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
+ c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
+ c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
+ c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
+ c[13]=f30*sy1*f40+f32*sy2*f42;
+ c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
+
+ UInt_t index=kr1.GetIndex(is);
+ AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
+ track->fIsSeeding = kTRUE;
+ Int_t rc=FollowProlongation(*track, i2);
+ //FollowProlongationFast(*track, 5);
+ //FollowProlongationFast(*track, 5);
+ //FollowProlongationFast(*track, 5);
+ //FollowProlongationFast(*track, 5);
+ //Int_t rc = 1;
+
+ track->fLastPoint = i1; // first cluster in track position
+ if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
+ else arr->AddLast(track);
+ }
+ }
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
+ //-----------------------------------------------------------------
+ // This function reades track seeds.
+ //-----------------------------------------------------------------
+ TDirectory *savedir=gDirectory;
+
+ TFile *in=(TFile*)inp;
+ if (!in->IsOpen()) {
+ cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
+ return 1;
+ }
+
+ in->cd();
+ TTree *seedTree=(TTree*)in->Get("Seeds");
+ if (!seedTree) {
+ cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
+ cerr<<"can't get a tree with track seeds !\n";
+ return 2;
+ }
+ AliTPCtrack *seed=new AliTPCtrack;
+ seedTree->SetBranchAddress("tracks",&seed);
+
+ if (fSeeds==0) fSeeds=new TObjArray(15000);
+
+ Int_t n=(Int_t)seedTree->GetEntries();
+ for (Int_t i=0; i<n; i++) {
+ seedTree->GetEvent(i);
+ fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
+ }
+
+ delete seed;
+ delete seedTree;
+ savedir->cd();
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
+ //-----------------------------------------------------------------
+ // This is a track finder.
+ //-----------------------------------------------------------------
+ TDirectory *savedir=gDirectory;
+
+ if (inp) {
+ TFile *in=(TFile*)inp;
+ if (!in->IsOpen()) {
+ cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
+ return 1;
+ }
+ }
+
+ if (!out->IsOpen()) {
+ cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
+ return 2;
+ }
+
+ out->cd();
+
+ char tname[100];
+ sprintf(tname,"TreeT_TPC_%d",fEventN);
+ TTree tracktree(tname,"Tree with TPC tracks");
+ TTree seedtree("Seeds","Seeds");
+ AliTPCtrack *iotrack=0;
+ AliTPCseed *ioseed=0;
+ tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
+ TStopwatch timer;
+
+ printf("Loading clusters \n");
+ LoadClusters();
+ printf("Time for loading clusters: \t");timer.Print();timer.Start();
+
+ printf("Loading outer sectors\n");
+ LoadOuterSectors();
+ printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
+
+ printf("Loading inner sectors\n");
+ LoadInnerSectors();
+ printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
+ fSectors = fOuterSec;
+ fN=fkNOS;
+
+
+ //find track seeds
+ MakeSeedsAll();
+ printf("Time for seeding: \t"); timer.Print();timer.Start();
+ Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
+ Int_t nrows=nlow+nup;
+
+ Int_t gap=Int_t(0.3*nrows);
+ Int_t i;
+ Int_t nseed=fSeeds->GetEntriesFast();
+ // outer sectors parallel tracking
+ ParallelTracking(fSectors->GetNRows()-gap-1,0);
+ printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
+
+ RemoveOverlap(fSeeds, 0.6,3);
+ printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
+ //parallel tracking
+ fSectors = fInnerSec;
+ fN=fkNIS;
+ ParallelTracking(fSectors->GetNRows()-1,0);
+ printf("Number of tracks after inner tracking %d\n",fNtracks);
+ printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
+ //
+ RemoveOverlap(fSeeds,0.6,5,kTRUE); // remove overlap - shared points signed
+ printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
+ //
+ //
+
+ ioseed = (AliTPCseed*)(fSeeds->UncheckedAt(0));
+ AliTPCseed * vseed = new AliTPCseed;
+ vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
+ vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
+ vseed->fPoints->ExpandCreateFast(2);
+
+ TBranch * seedbranch = seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
+ //delete vseed;
+ nseed=fSeeds->GetEntriesFast();
+
+ Int_t found = 0;
+ for (i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
+ if (!pt) continue;
+ Int_t nc=t.GetNumberOfClusters();
+ t.CookdEdx(0.02,0.6);
+ CookLabel(pt,0.1); //For comparison only
+ // if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
+ if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (nc>Int_t(0.3*nrows)))){
+ iotrack=pt;
+ tracktree.Fill();
+ cerr<<found++<<'\r';
+ }
+ else
+ if ( (pt->IsActive())) fNtracks--;
+ pt->RebuildSeed();
+ seedbranch->SetAddress(&pt);
+
+ seedtree.Fill();
+ for (Int_t j=0;j<160;j++){
+ delete pt->fPoints->RemoveAt(j);
+ }
+ delete pt->fPoints;
+ pt->fPoints =0;
+ delete fSeeds->RemoveAt(i);
+ }
+ printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
+
+ UnloadClusters();
+ printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
+
+ tracktree.Write();
+ seedtree.Write();
+ cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
+
+ savedir->cd();
+
+ return 0;
+}
+
+
+void AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
+{
+ //
+ // try to track in parralel
+
+ Int_t nseed=fSeeds->GetEntriesFast();
+ //prepare seeds for tracking
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
+ if (!pt) continue;
+ if (!t.IsActive()) continue;
+ // follow prolongation to the first layer
+ FollowProlongation(t, rfirst+1);
+ }
+
+
+ //
+ for (Int_t nr=rfirst; nr>=rlast; nr--){
+ // make indexes with the cluster tracks for given
+ // for (Int_t i = 0;i<fN;i++)
+ // fSectors[i][nr].MakeClusterTracks();
+
+ // find nearest cluster
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
+ if (!pt) continue;
+ if (!pt->IsActive()) continue;
+ UpdateClusters(t,i,nr);
+ }
+ // prolonagate to the nearest cluster - if founded
+ for (Int_t i=0; i<nseed; i++) {
+ AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
+ if (!pt) continue;
+ if (!pt->IsActive()) continue;
+ FollowToNextCluster(i,nr);
+ }
+ // for (Int_t i= 0;i<fN;i++)
+ // fSectors[i][nr].ClearClusterTracks();
+ }
+
+}
+
+Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
+{
+ //
+ //
+ Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
+ Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
+ Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
+ Float_t angular = seed->GetSnp();
+ angular = angular*angular/(1-angular*angular);
+ // angular*=angular;
+ //angular = TMath::Sqrt(angular/(1-angular));
+ Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
+ return res;
+}
+Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
+{
+ //
+ //
+ Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
+ Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
+ Float_t sres = fParam->GetZSigma();
+ Float_t angular = seed->GetTgl();
+ Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
+ return res;
+}
+
+
+
+//_________________________________________________________________________
+AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
+ //--------------------------------------------------------------------
+ // Return pointer to a given cluster
+ //--------------------------------------------------------------------
+ Int_t sec=(index&0xff000000)>>24;
+ Int_t row=(index&0x00ff0000)>>16;
+ Int_t ncl=(index&0x0000ffff)>>00;
+
+ AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
+ return (AliTPCclusterMI*)(*clrow)[ncl];
+}
+
+//__________________________________________________________________________
+void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
+ //--------------------------------------------------------------------
+ //This function "cooks" a track label. If label<0, this track is fake.
+ //--------------------------------------------------------------------
+ Int_t noc=t->GetNumberOfClusters();
+ Int_t *lb=new Int_t[noc];
+ Int_t *mx=new Int_t[noc];
+ AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
+
+ Int_t i;
+ for (i=0; i<noc; i++) {
+ lb[i]=mx[i]=0;
+ Int_t index=t->GetClusterIndex(i);
+ clusters[i]=GetClusterMI(index);
+ }
+
+ Int_t lab=123456789;
+ for (i=0; i<noc; i++) {
+ AliTPCclusterMI *c=clusters[i];
+ lab=TMath::Abs(c->GetLabel(0));
+ Int_t j;
+ for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
+ lb[j]=lab;
+ (mx[j])++;
+ }
+
+ Int_t max=0;
+ for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
+
+ for (i=0; i<noc; i++) {
+ AliTPCclusterMI *c=clusters[i];
+ if (TMath::Abs(c->GetLabel(1)) == lab ||
+ TMath::Abs(c->GetLabel(2)) == lab ) max++;
+ }
+
+ if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+
+ else {
+ Int_t tail=Int_t(0.10*noc);
+ max=0;
+ for (i=1; i<=tail; i++) {
+ AliTPCclusterMI *c=clusters[noc-i];
+ if (lab == TMath::Abs(c->GetLabel(0)) ||
+ lab == TMath::Abs(c->GetLabel(1)) ||
+ lab == TMath::Abs(c->GetLabel(2))) max++;
+ }
+ if (max < Int_t(0.5*tail)) lab=-lab;
+ }
+
+ t->SetLabel(lab);
+
+ delete[] lb;
+ delete[] mx;
+ delete[] clusters;
+}
+
+//_________________________________________________________________________
+void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
+ //-----------------------------------------------------------------------
+ // Setup inner sector
+ //-----------------------------------------------------------------------
+ if (f==0) {
+ fAlpha=par->GetInnerAngle();
+ fAlphaShift=par->GetInnerAngleShift();
+ fPadPitchWidth=par->GetInnerPadPitchWidth();
+ fPadPitchLength=par->GetInnerPadPitchLength();
+ fN=par->GetNRowLow();
+ fRow=new AliTPCRow[fN];
+ for (Int_t i=0; i<fN; i++) {
+ fRow[i].SetX(par->GetPadRowRadiiLow(i));
+ fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
+ }
+ } else {
+ fAlpha=par->GetOuterAngle();
+ fAlphaShift=par->GetOuterAngleShift();
+ fPadPitchWidth = par->GetOuterPadPitchWidth();
+ fPadPitchLength = par->GetOuter1PadPitchLength();
+ f1PadPitchLength = par->GetOuter1PadPitchLength();
+ f2PadPitchLength = par->GetOuter2PadPitchLength();
+
+ fN=par->GetNRowUp();
+ fRow=new AliTPCRow[fN];
+ for (Int_t i=0; i<fN; i++) {
+ fRow[i].SetX(par->GetPadRowRadiiUp(i));
+ fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
+ }
+ }
+}
+
+
+AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
+ //
+ if (fClusterTracks) delete [] fClusterTracks;
+ fClusterTracks = 0;
+}
+
+void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
+ //create cluster tracks
+ if (fN>0)
+ fClusterTracks = new AliTPCclusterTracks[fN];
+}
+
+void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
+ if (fClusterTracks) delete[] fClusterTracks;
+ fClusterTracks =0;
+}
+
+
+
+void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
+ //
+ //
+ // update information of the cluster tracks - if track is nearer then other tracks to the
+ // given track
+ const AliTPCclusterMI * cl = (*this)[clindex];
+ AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
+ // find the distance of the cluster to the track
+ Float_t dy2 = (cl->GetY()- seed->GetY());
+ dy2*=dy2;
+ Float_t dz2 = (cl->GetZ()- seed->GetZ());
+ dz2*=dz2;
+ //
+ Float_t distance = TMath::Sqrt(dy2+dz2);
+ if (distance> 3)
+ return; // MI - to be changed - AliTPCtrackerParam
+
+ if ( distance < cltracks->fDistance[0]){
+ cltracks->fDistance[2] =cltracks->fDistance[1];
+ cltracks->fDistance[1] =cltracks->fDistance[0];
+ cltracks->fDistance[0] =distance;
+ cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
+ cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
+ cltracks->fTrackIndex[0] =trindex;
+ }
+ else
+ if ( distance < cltracks->fDistance[1]){
+ cltracks->fDistance[2] =cltracks->fDistance[1];
+ cltracks->fDistance[1] =distance;
+ cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
+ cltracks->fTrackIndex[1] =trindex;
+ } else
+ if (distance < cltracks->fDistance[2]){
+ cltracks->fDistance[2] =distance;
+ cltracks->fTrackIndex[2] =trindex;
+ }
+}
+
+
+//_________________________________________________________________________
+void
+AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
+ //-----------------------------------------------------------------------
+ // Insert a cluster into this pad row in accordence with its y-coordinate
+ //-----------------------------------------------------------------------
+ if (fN==kMaxClusterPerRow) {
+ cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
+ }
+ if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
+ Int_t i=Find(c->GetZ());
+ memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
+ memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
+ fIndex[i]=index; fClusters[i]=c; fN++;
+}
+
+//___________________________________________________________________
+Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
+ //-----------------------------------------------------------------------
+ // Return the index of the nearest cluster
+ //-----------------------------------------------------------------------
+ if (fN==0) return 0;
+ if (z <= fClusters[0]->GetZ()) return 0;
+ if (z > fClusters[fN-1]->GetZ()) return fN;
+ Int_t b=0, e=fN-1, m=(b+e)/2;
+ for (; b<e; m=(b+e)/2) {
+ if (z > fClusters[m]->GetZ()) b=m+1;
+ else e=m;
+ }
+ return m;
+}
+
+
+
+AliTPCseed::AliTPCseed():AliTPCtrack(){
+ //
+ fRow=0;
+ fRemoval =0;
+ memset(fClusterIndex,0,sizeof(Int_t)*200);
+ fPoints = 0;
+ fEPoints = 0;
+ fNFoundable =0;
+ fNShared =0;
+ fTrackPoints =0;
+ fRemoval = 0;
+}
+
+AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
+ fPoints = 0;
+ fEPoints = 0;
+ fNShared =0;
+ fTrackPoints =0;
+ fRemoval =0;
+}
+
+AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
+ fRow=0;
+ memset(fClusterIndex,0,sizeof(Int_t)*200);
+ fPoints = 0;
+ fEPoints = 0;
+ fNFoundable =0;
+ fNShared =0;
+ fTrackPoints =0;
+ fRemoval =0;
+}
+
+AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
+ Double_t xr, Double_t alpha):
+ AliTPCtrack(index, xx, cc, xr, alpha) {
+ //
+ //
+ fRow =0;
+ memset(fClusterIndex,0,sizeof(Int_t)*200);
+ fPoints = 0;
+ fEPoints = 0;
+ fNFoundable =0;
+ fNShared = 0;
+ fTrackPoints =0;
+ fRemoval =0;
+}
+
+AliTPCseed::~AliTPCseed(){
+ if (fPoints) delete fPoints;
+ fPoints =0;
+ fEPoints = 0;
+ if (fTrackPoints){
+ for (Int_t i=0;i<8;i++){
+ delete [] fTrackPoints[i];
+ }
+ delete fTrackPoints;
+ fTrackPoints =0;
+ }
+
+}
+
+AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
+{
+ //
+ //
+ if (!fTrackPoints) {
+ fTrackPoints = new AliTPCTrackPoint*[8];
+ for ( Int_t i=0;i<8;i++)
+ fTrackPoints[i]=0;
+ }
+ Int_t index1 = i/20;
+ if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
+ return &(fTrackPoints[index1][i%20]);
+}
+
+void AliTPCseed::RebuildSeed()
+{
+ //
+ // rebuild seed to be ready for storing
+ fPoints = new TClonesArray("AliTPCTrackPoint",160);
+ fPoints->ExpandCreateFast(160);
+ fEPoints = new TClonesArray("AliTPCExactPoint",1);
+ for (Int_t i=0;i<160;i++){
+ AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
+ *trpoint = *(GetTrackPoint(i));
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
+ //-----------------------------------------------------------------
+ // This funtion calculates dE/dX within the "low" and "up" cuts.
+ //-----------------------------------------------------------------
+
+ Float_t amp[200];
+ Float_t angular[200];
+ Float_t weight[200];
+ Int_t index[200];
+ //Int_t nc = 0;
+ // TClonesArray & arr = *fPoints;
+ Float_t meanlog = 100.;
+
+ Float_t mean[4] = {0,0,0,0};
+ Float_t sigma[4] = {1000,1000,1000,1000};
+ Int_t nc[4] = {0,0,0,0};
+ Float_t norm[4] = {1000,1000,1000,1000};
+ //
+ //
+ fNShared =0;
+
+ for (Int_t of =0; of<4; of++){
+ for (Int_t i=of;i<160;i+=4)
+ {
+ //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
+ AliTPCTrackPoint * point = GetTrackPoint(i);
+ if (point==0) continue;
+ if (point->fIsShared){
+ fNShared++;
+ continue;
+ }
+ if (point->GetCPoint().GetMax()<5) continue;
+ Float_t angley = point->GetTPoint().GetAngleY();
+ Float_t anglez = point->GetTPoint().GetAngleZ();
+ Int_t type = point->GetCPoint().GetType();
+ Float_t rsigmay = point->GetCPoint().GetSigmaY();
+ Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
+ Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
+
+ Float_t ampc = 0; // normalization to the number of electrons
+ if (i>64){
+ ampc = 1.*point->GetCPoint().GetMax();
+ //ampc = 1.*point->GetCPoint().GetQ();
+ // AliTPCClusterPoint & p = point->GetCPoint();
+ // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
+ // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+ //Float_t dz =
+ // TMath::Abs( Int_t(iz) - iz + 0.5);
+ //ampc *= 1.15*(1-0.3*dy);
+ //ampc *= 1.15*(1-0.3*dz);
+
+ }
+ else{
+ ampc = 1.0*point->GetCPoint().GetMax();
+ //ampc = 1.0*point->GetCPoint().GetQ();
+ //AliTPCClusterPoint & p = point->GetCPoint();
+ // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
+ //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+ //Float_t dz =
+ // TMath::Abs( Int_t(iz) - iz + 0.5);
+
+ //ampc *= 1.15*(1-0.3*dy);
+ //ampc *= 1.15*(1-0.3*dz);
+
+ }
+ ampc *= 2.0; // put mean value to channel 50
+ //ampc *= 0.58; // put mean value to channel 50
+ Float_t w = 1.;
+ // if (type>0) w = 1./(type/2.-0.5);
+ if (i<64) {
+ ampc /= 0.6;
+ }
+ if (i>128){
+ ampc /=1.5;
+ }
+ if (type<0) { //amp at the border - lower weight
+ // w*= 2.;
+ continue;
+ }
+ if (rsigma>1.5) ampc/=1.3; // if big backround
+ amp[nc[of]] = ampc;
+ angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
+ weight[nc[of]] = w;
+ nc[of]++;
+ }
+
+ TMath::Sort(nc[of],amp,index,kFALSE);
+ Float_t sumamp=0;
+ Float_t sumamp2=0;
+ Float_t sumw=0;
+ //meanlog = amp[index[Int_t(nc[of]*0.33)]];
+ meanlog = 200;
+ for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
+ Float_t ampl = amp[index[i]]/angular[index[i]];
+ ampl = meanlog*TMath::Log(1.+ampl/meanlog);
+ //
+ sumw += weight[index[i]];
+ sumamp += weight[index[i]]*ampl;
+ sumamp2 += weight[index[i]]*ampl*ampl;
+ norm[of] += angular[index[i]]*weight[index[i]];
+ }
+ if (sumw<1){
+ SetdEdx(0);
+ }
+ else {
+ norm[of] /= sumw;
+ mean[of] = sumamp/sumw;
+ sigma[of] = sumamp2/sumw-mean[of]*mean[of];
+ if (sigma[of]>0.1)
+ sigma[of] = TMath::Sqrt(sigma[of]);
+ else
+ sigma[of] = 1000;
+
+ mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
+ //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
+ //mean *=(1-0.1*(norm-1.));
+ }
+ }
+
+ Float_t dedx =0;
+ fSdEdx =0;
+ fMAngular =0;
+ // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
+ // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
+
+
+ // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
+ // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
+
+ Int_t norm2 = 0;
+ Int_t norm3 = 0;
+ for (Int_t i =0;i<4;i++){
+ if (nc[i]>2&&nc[i]<1000){
+ dedx += mean[i] *nc[i];
+ fSdEdx += sigma[i]*(nc[i]-2);
+ fMAngular += norm[i] *nc[i];
+ norm2 += nc[i];
+ norm3 += nc[i]-2;
+ }
+ fDEDX[i] = mean[i];
+ fSDEDX[i] = sigma[i];
+ fNCDEDX[i]= nc[i];
+ }
+
+ if (norm3>0){
+ dedx /=norm2;
+ fSdEdx /=norm3;
+ fMAngular/=norm2;
+ }
+ else{
+ SetdEdx(0);
+ return;
+ }
+ // Float_t dedx1 =dedx;
+
+ dedx =0;
+ for (Int_t i =0;i<4;i++){
+ if (nc[i]>2&&nc[i]<1000){
+ mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
+ dedx += mean[i] *nc[i];
+ }
+ fDEDX[i] = mean[i];
+ }
+ dedx /= norm2;
+
+
+
+ SetdEdx(dedx);
+
+ //mi deDX
+
+
+
+ //Very rough PID
+ Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
+
+ if (p<0.6) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
+ SetMass(0.93827); return;
+ }
+
+ if (p<1.2) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ SetMass(0.93827); return;
+ }
+
+ SetMass(0.13957); return;
+
+}
+
+
--- /dev/null
+#ifndef ALITPCTRACKERMI_H
+#define ALITPCTRACKERMI_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id$ */
+
+//-------------------------------------------------------
+// TPC trackerMI
+//
+// Origin:
+//-------------------------------------------------------
+#include "AliTracker.h"
+#include "AliTPCtrack.h"
+#include "AliTPCClustersArray.h"
+
+#include "AliTPCreco.h"
+#include "Rtypes.h"
+
+class TFile;
+class AliTPCParam;
+class AliTPCseed;
+class AliTPCclusterMI;
+class AliTPCTrackPoint;
+
+
+
+class AliTPCclusterTracks {
+ public:
+ AliTPCclusterTracks();
+ Float_t fDistance[3]; // distance to the 3 nerest track if there overlap with cluster
+ Short_t fTrackIndex[3]; // indexes of the tracks overlapped with clusters
+};
+
+
+class AliTPCseed : public AliTPCtrack {
+ public:
+ AliTPCseed();
+ virtual ~AliTPCseed();
+ AliTPCseed(const AliTPCtrack &t);
+ AliTPCseed(const AliKalmanTrack &t, Double_t a);
+ Int_t Compare(const TObject *o) const;
+
+ Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const;
+ virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const;
+ virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i);
+ AliTPCTrackPoint * GetTrackPoint(Int_t i);
+ void RebuildSeed(); // rebuild seed to be ready for storing
+ AliTPCseed(UInt_t index, const Double_t xx[5],
+ const Double_t cc[15], Double_t xr, Double_t alpha);
+ void SetClusterIndex(Int_t index){
+ fClusterIndex[fRow] = index;
+ }
+ void SetErrorY2(Float_t sy2){fErrorY2=sy2;}
+ void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;}
+ void CookdEdx(Double_t low=0.05, Double_t up=0.70);
+ Bool_t IsActive(){ return !(fRemoval);}
+ void Desactivate(Int_t reason){ fRemoval = reason;}
+ Int_t fRelativeSector; // ! index of current relative sector
+ Int_t fClusterIndex[200]; //array of cluster indexes
+
+ Int_t fRemoval; //reason - why was track removed - 0 - means still active
+ TClonesArray * fPoints; // array with points along the track
+ TClonesArray * fEPoints; // array with exact points - calculated in special macro not used in tracking
+ Int_t fRow; //! current row number
+ Int_t fSector; //!current sector number
+ Float_t fCurrentSigmaY; //!expected current cluster sigma Y
+ Float_t fCurrentSigmaZ; //!expected current cluster sigma Z
+ AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
+ Int_t fCurrentClusterIndex1; //! index of the current cluster
+ Int_t fCurrentClusterIndex2; //! index of the current cluster
+
+ Float_t fErrorY2; //!sigma of current cluster
+ Float_t fErrorZ2; //!sigma of current cluster
+ Int_t fNFoundable; //number of foundable clusters - dead zone taken to the account
+ Bool_t fInDead; // indicate if the track is in dead zone
+ Int_t fFirstPoint; // first cluster position
+ Int_t fLastPoint; // last cluster position
+ Int_t fNShared; // number of shared points
+ Bool_t fIsSeeding; //indicates if it is proces of seeading
+ private:
+ Float_t fSdEdx; // sigma of dedx
+ Float_t fMAngular; // mean angular factor
+ AliTPCTrackPoint ** fTrackPoints; //!track points - array track points
+ Float_t fDEDX[4]; // dedx according padrows
+ Float_t fSDEDX[4]; // sdedx according padrows
+ Int_t fNCDEDX[4]; // number of clusters for dedx measurment
+
+ ClassDef(AliTPCseed,1)
+};
+
+
+
+
+class AliTPCtrackerMI : public AliTracker {
+public:
+ AliTPCtrackerMI():AliTracker(),fkNIS(0),fkNOS(0) {
+ fInnerSec=fOuterSec=0; fSeeds=0;
+ }
+ AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn=0);
+ ~AliTPCtrackerMI();
+
+ Int_t ReadSeeds(const TFile *in);
+ void LoadClusters();
+ void UnloadClusters();
+
+ void LoadInnerSectors();
+ void LoadOuterSectors();
+ AliCluster * GetCluster (int) const {return 0;}
+ AliTPCclusterMI *GetClusterMI(Int_t index) const;
+ Int_t Clusters2Tracks(const TFile *in, TFile *out);
+ // Int_t PropagateBack(const TFile *in, TFile *out);
+
+ virtual void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
+
+ virtual Double_t ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl = 0);
+ virtual Double_t ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl = 0);
+
+ Double_t f1(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t x3,Double_t y3);
+ Double_t f2(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t x3,Double_t y3);
+ Double_t f3(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t z1,Double_t z2);
+public:
+//**************** Internal tracker class **********************
+ class AliTPCRow {
+ public:
+ AliTPCRow() {fN=0; fClusterTracks=0;}
+ ~AliTPCRow();
+ void InsertCluster(const AliTPCclusterMI *c, UInt_t index);
+ operator int() const {return fN;}
+ const AliTPCclusterMI* operator[](Int_t i) const {return fClusters[i];}
+ UInt_t GetIndex(Int_t i) const {return fIndex[i];}
+ Int_t Find(Double_t z) const;
+ void SetX(Double_t x) {fX=x;}
+ Double_t GetX() const {return fX;}
+ AliTPCclusterTracks * GetClusterTracks(Int_t index){ return ( (index<fN) && fClusterTracks!=0)? &(fClusterTracks[index]):0;}
+ void UpdateClusterTrack(Int_t clindex, Int_t trindex,AliTPCseed * seed);
+ void MakeClusterTracks();
+ void ClearClusterTracks();
+ Float_t fDeadZone; // the width of the dead zone
+ private:
+ Int_t fN; //number of clusters
+ const AliTPCclusterMI *fClusters[kMaxClusterPerRow]; //pointers to clusters
+ UInt_t fIndex[kMaxClusterPerRow]; //indeces of clusters
+ Double_t fX; //X-coordinate of this row
+ AliTPCclusterTracks * fClusterTracks; // array of cluster tracks - for overlap calculation
+ private:
+ AliTPCRow(const AliTPCRow& r); //dummy copy constructor
+ AliTPCRow &operator=(const AliTPCRow& r); //dummy assignment operator
+ };
+
+//**************** Internal tracker class **********************
+ class AliTPCSector {
+ public:
+ AliTPCSector() { fN=0; fRow = 0; }
+ ~AliTPCSector() { delete[] fRow; }
+ AliTPCRow& operator[](Int_t i) const { return *(fRow+i); }
+ Int_t GetNRows() const { return fN; }
+ void Setup(const AliTPCParam *par, Int_t flag);
+ Double_t GetX(Int_t l) const {return fRow[l].GetX();}
+ Double_t GetMaxY(Int_t l) const {
+ return GetX(l)*TMath::Tan(0.5*GetAlpha());
+ }
+ Double_t GetAlpha() const {return fAlpha;}
+ Double_t GetAlphaShift() const {return fAlphaShift;}
+ Int_t GetRowNumber(Double_t x) const {
+ //return pad row number for this x
+ Double_t r;
+ if (fN < 64){
+ r=fRow[fN-1].GetX();
+ if (x > r) return fN;
+ r=fRow[0].GetX();
+ if (x < r) return -1;
+ return Int_t((x-r)/fPadPitchLength + 0.5);}
+ else{
+ r=fRow[fN-1].GetX();
+ if (x > r) return fN;
+ r=fRow[0].GetX();
+ if (x < r) return -1;
+ Double_t r1=fRow[64].GetX();
+ if(x<r1){
+ return Int_t((x-r)/f1PadPitchLength + 0.5);}
+ else{
+ return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
+ }
+ }
+ Double_t GetPadPitchWidth() const {return fPadPitchWidth;}
+ Double_t GetPadPitchLength() const {return fPadPitchLength;}
+ Double_t GetPadPitchLength(Float_t x) const {return (x<200) ? fPadPitchLength:f2PadPitchLength ;}
+
+ private:
+ Int_t fN; //number of pad rows
+ AliTPCRow *fRow; //array of pad rows
+ Double_t fAlpha; //opening angle
+ Double_t fAlphaShift; //shift angle;
+ Double_t fPadPitchWidth; //pad pitch width
+ Double_t fPadPitchLength; //pad pitch length
+ Double_t f1PadPitchLength; //pad pitch length
+ Double_t f2PadPitchLength; //pad pitch length
+
+ private:
+ AliTPCSector(const AliTPCSector &s); //dummy copy contructor
+ AliTPCSector& operator=(const AliTPCSector &s);//dummy assignment operator
+ };
+
+ Float_t OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t &sum2);
+ void SignShared(AliTPCseed * s1, AliTPCseed * s2);
+ void RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex, Bool_t shared=kFALSE);
+private:
+ Float_t GetSigmaY(AliTPCseed * seed);
+ Float_t GetSigmaZ(AliTPCseed * seed);
+
+ void MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2);
+ TObjArray * MakeSeedsSectors(Int_t sec1, Int_t sec2); // make seeds from all sectors
+ void MakeSeedsAll();
+ Int_t FollowProlongation(AliTPCseed& t, Int_t rf=0);
+ //Int_t FollowProlongationFast(AliTPCseed& t, Int_t step);
+ Int_t FollowBackProlongation(AliTPCseed& t, Int_t rf=0);
+
+ Int_t FollowToNext(AliTPCseed& t, Int_t nr);
+ Int_t UpdateClusters(AliTPCseed& t, Int_t trindex, Int_t nr);
+ Int_t FollowToNextCluster( Int_t trindex, Int_t nr);
+
+ virtual Int_t PropagateBack (const TFile *, TFile *){return 0;}
+ void ParallelTracking(Int_t rfirst, Int_t rlast);
+ void SetSampledEdx(AliTPCseed *t, Float_t q, Int_t i) {;}
+ Int_t UpdateTrack(AliTPCseed *t,AliTPCclusterMI* c, Double_t chi2, UInt_t i); //update trackinfo
+
+ // Int_t FollowBackProlongation(AliTPCseed &s, const AliTPCtrack &t);
+
+ AliTPCtrackerMI(const AliTPCtrackerMI& r); //dummy copy constructor
+ AliTPCtrackerMI &operator=(const AliTPCtrackerMI& r);//dummy assignment operator
+
+ const Int_t fkNIS; //number of inner sectors
+ AliTPCSector *fInnerSec; //array of inner sectors;
+ const Int_t fkNOS; //number of outer sectors
+ AliTPCSector *fOuterSec; //array of outer sectors;
+
+ Int_t fN; //number of loaded sectors
+ AliTPCSector *fSectors; //pointer to loaded sectors;
+
+ Int_t fEventN; //event number
+ AliTPCClustersArray fClustersArray; //array of TPC clusters
+ Int_t fNtracks; //current number of tracks
+ TObjArray *fSeeds; //array of track seeds
+ // TObjArray * fTrackPointPool; // ! pool with track points
+ const AliTPCParam *fParam; //pointer to the parameters
+};
+
+#endif
+
+
--- /dev/null
+#include "TH1.h"
+#include "TH2.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include "TRandom.h"
+#include "TPad.h"
+#include "TCanvas.h"
+
+
+class TLandauMean: public TObject {
+public:
+ void Init(Int_t n, Float_t mean, Float_t sigma); // initial parameters
+ void Gener(); // gener sample
+ // void Anal();
+ //
+ Int_t fNSample; // number of samples
+ Float_t fLMean; // landau mean
+ Float_t fLSigma; // landau sigma
+ //
+ Float_t fTM_0_6[3]; // truncated method - first 3 momenta
+ Float_t fTM_0_7[3]; // truncated method - first 3 momenta
+ Float_t fTM_0_8[3]; // truncated method - first 3 momenta
+ Float_t fTM_0_10[3]; // truncated method - first 3 momenta
+ //
+ Float_t fLM_0_6[3]; // truncated log. method - first 3 momenta
+ Float_t fLM_0_7[3]; // truncated log. method - first 3 momenta
+ Float_t fLM_0_8[3]; // truncated log. method - first 3 momenta
+ Float_t fLM_0_10[3]; // truncated log. method - first 3 momenta
+
+ Float_t fMedian3; // median 3 value
+private:
+ Float_t Moment3(Float_t sum1, Float_t sum2, Float_t sum3, Int_t n, Float_t m[3]);
+ ClassDef(TLandauMean,1)
+};
+
+ClassImp(TLandauMean)
+
+void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
+{
+ //
+ //init parameters
+ fNSample = n;
+ fLMean = mean;
+ fLSigma = sigma;
+}
+
+Float_t TLandauMean::Moment3(Float_t sumi1, Float_t sumi2, Float_t sumi3, Int_t sum, Float_t m[3])
+{
+ Float_t m3=0;
+
+ // m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi-pos*pos*pos*sum)/sum;
+ Float_t pos = sumi1/sum;
+ m[0] = pos;
+ m[1] = sumi2/sum-pos*pos;
+ if (m[1]<=0){
+ printf("pici pici\n");
+ }
+ else
+ m[1] = TMath::Sqrt(m[1]);
+ m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi1-pos*pos*pos*sum)/sum;
+ Float_t sign = m3/TMath::Abs(m3);
+ m3 = TMath::Power(sign*m3,1/3.);
+ m3*=sign;
+
+ m[2] = m3;
+ return m3;
+}
+
+void TLandauMean::Gener()
+{
+ //
+ // generate sample
+ Float_t * buffer = new Float_t[fNSample];
+
+ for (Int_t i=0;i<fNSample;i++) {
+ buffer[i] = gRandom->Landau(fLMean,fLSigma);
+ if (buffer[i]>1000) buffer[i]=1000;
+ }
+
+ Int_t *index = new Int_t[fNSample];
+ TMath::Sort(fNSample,buffer,index,kFALSE);
+
+ //
+ Float_t median = buffer[index[fNSample/3]];
+ //
+ Float_t sum06[4] = {0.,0.,0.,0.};
+ Float_t sum07[4] = {0.,0.,0.,0.};
+ Float_t sum08[4] = {0.,0.,0.,0.};
+ Float_t sum010[4] = {0.,0.,0.,0.};
+ //
+ Float_t suml06[4] = {0.,0.,0.,0.};
+ Float_t suml07[4] = {0.,0.,0.,0.};
+ Float_t suml08[4] = {0.,0.,0.,0.};
+ Float_t suml010[4] = {0.,0.,0.,0.};
+ //
+
+ for (Int_t i =0; i<fNSample; i++){
+ Float_t amp = buffer[index[i]];
+ Float_t lamp = median*TMath::Log(1.+amp/median);
+ if (i<0.6*fNSample){
+ sum06[0]+= amp;
+ sum06[1]+= amp*amp;
+ sum06[2]+= amp*amp*amp;
+ sum06[3]++;
+ suml06[0]+= lamp;
+ suml06[1]+= lamp*lamp;
+ suml06[2]+= lamp*lamp*lamp;
+ suml06[3]++;
+ }
+
+ if (i<0.7*fNSample){
+ sum07[0]+= amp;
+ sum07[1]+= amp*amp;
+ sum07[2]+= amp*amp*amp;
+ sum07[3]++;
+ suml07[0]+= lamp;
+ suml07[1]+= lamp*lamp;
+ suml07[2]+= lamp*lamp*lamp;
+ suml07[3]++;
+ }
+ if (i<0.8*fNSample){
+ sum08[0]+= amp;
+ sum08[1]+= amp*amp;
+ sum08[2]+= amp*amp*amp;
+ sum08[3]++;
+ suml08[0]+= lamp;
+ suml08[1]+= lamp*lamp;
+ suml08[2]+= lamp*lamp*lamp;
+ suml08[3]++;
+ }
+ if (i<1*fNSample){
+ sum010[0]+= amp;
+ sum010[1]+= amp*amp;
+ sum010[2]+= amp*amp*amp;
+ sum010[3]++;
+ suml010[0]+= lamp;
+ suml010[1]+= lamp*lamp;
+ suml010[2]+= lamp*lamp*lamp;
+ suml010[3]++;
+
+ }
+ }
+ //
+ fMedian3 = median;
+ //
+ Moment3(sum06[0],sum06[1],sum06[2],sum06[3],fTM_0_6);
+ Moment3(sum07[0],sum07[1],sum07[2],sum07[3],fTM_0_7);
+ Moment3(sum08[0],sum08[1],sum08[2],sum08[3],fTM_0_8);
+ Moment3(sum010[0],sum010[1],sum010[2],sum010[3],fTM_0_10);
+ //
+
+ Moment3(suml06[0],suml06[1],suml06[2],suml06[3],fLM_0_6);
+ Moment3(suml07[0],suml07[1],suml07[2],suml07[3],fLM_0_7);
+ Moment3(suml08[0],suml08[1],suml08[2],suml08[3],fLM_0_8);
+ Moment3(suml010[0],suml010[1],suml010[2],suml010[3],fLM_0_10);
+ //
+ fLM_0_6[0] = (TMath::Exp(fLM_0_6[0]/median)-1.)*median;
+ fLM_0_7[0] = (TMath::Exp(fLM_0_7[0]/median)-1.)*median;
+ fLM_0_8[0] = (TMath::Exp(fLM_0_8[0]/median)-1.)*median;
+ fLM_0_10[0] = (TMath::Exp(fLM_0_10[0]/median)-1.)*median;
+ //
+ delete [] buffer;
+}
+
+
+void GenerLandau(Int_t nsamples)
+{
+ TLandauMean * landau = new TLandauMean;
+ TFile f("Landau.root","recreate");
+ TTree * tree = new TTree("Landau","Landau");
+ tree->Branch("Landau","TLandauMean",&landau);
+
+ for (Int_t i=0;i<nsamples;i++){
+ Int_t n = 20 + Int_t(gRandom->Rndm()*150);
+ Float_t mean = 40. +gRandom->Rndm()*50.;
+ Float_t sigma = 5. +gRandom->Rndm()*15.;
+ landau->Init(n, mean, sigma);
+ landau->Gener();
+ tree->Fill();
+ }
+ tree->Write();
+ f.Close();
+
+}
+
+
+
+
+
+TH1F * LandauTest(Float_t meano, Float_t sigma, Float_t meanlog0, Int_t n,Float_t ratio)
+{
+ //
+ // test for different approach of de dx resolution
+ // meano, sigma - mean value of Landau distribution and sigma
+ // meanlog0 - scaling factor for logarithmic mean value
+ // n - number of used layers
+ // ratio - ratio of used amplitudes for truncated mean
+ //
+
+
+ TCanvas * pad = new TCanvas("Landau test");
+ pad->Divide(2,2);
+ TH1F * h1 = new TH1F("h1","Logarithmic mean",300,0,4*meano);
+ TH1F * h2 = new TH1F("h2","Logarithmic amplitudes",300,0,8*meano);
+ TH1F * h3 = new TH1F("h3","Mean",300,0,4*meano);
+ TH1F * h4 = new TH1F("h4","Amplitudes",300,0,8*meano);
+
+ for(Int_t j=0;j<10000;j++){
+ //generate sample and sort it
+ Float_t * buffer = new Float_t[n];
+ Float_t * buffer2= new Float_t[n];
+
+ for (Int_t i=0;i<n;i++) {
+ buffer[i] = gRandom->Landau(meano,sigma);
+ buffer2[i] = buffer[i];
+ }
+ //add crosstalk
+ for (Int_t i=1;i<n-1;i++) {
+ buffer[i] = buffer2[i]*1.0+ buffer2[i-1]*0.0+ buffer2[i+1]*0.0;
+ buffer[i] = TMath::Min(buffer[i],1000.);
+ }
+ Int_t *index = new Int_t[n];
+ TMath::Sort(n,buffer,index,kFALSE);
+
+ //calculate mean
+ Float_t sum;
+ sum=0;
+ Float_t mean;
+ Float_t used = 0;
+ for (Int_t i=0;i<n*ratio;i++) {
+ if (buffer[index[i]]<1000.){
+ Float_t amp = meanlog0*TMath::Log(1+buffer[index[i]]/meanlog0);
+ sum += amp;
+ used++;
+ }
+ }
+ mean = sum/used;
+ //
+ sum=0;
+ used=0;
+ Float_t sum2=0;
+ Float_t meanlog =meanlog0;
+ for (Int_t i=0;i<n*ratio;i++) {
+ if (buffer[index[i]]<1000.){
+ Float_t amp = meanlog*TMath::Log(1.+buffer[index[i]]/(meanlog));
+ sum +=amp;
+ sum2+=buffer[index[i]];
+ used++;
+ h2->Fill(amp);
+ h4->Fill(buffer[index[i]]);
+ }
+ }
+ mean = sum/used;
+ mean = (TMath::Exp(mean/meanlog)-1)*meanlog;
+ Float_t mean2 = sum2/used;
+ //mean2 = (mean+mean2)/2.;
+ h1->Fill(mean);
+ h3->Fill(mean2);
+ }
+
+ pad->cd(1);
+ h1->Draw();
+ pad->cd(2);
+ h2->Draw();
+ pad->cd(3);
+ h3->Draw();
+ pad->cd(4);
+ h4->Draw();
+
+
+ return h1;
+
+}
+
#pragma link C++ class AliTPCtrack+;
#pragma link C++ class AliTPCtracker+;
+
#pragma link C++ class AliTPCParam+;
#pragma link C++ class AliTPCParamSR-;
#pragma link C++ class AliTPCParamCR-;
#pragma link C++ class AliComplexCluster+;
#pragma link C++ class AliDigitCluster+;
-#pragma link C++ class AliDifCluster+;
#pragma link C++ class AliClusters+;
#pragma link C++ class AliClustersArray+;
#pragma link C++ class AliTPCDigitsArray+;
-#pragma link C++ class AliTPCClusterFinder+;
#pragma link C++ class AliH2F+;
#pragma link C++ class AliTPCTrackHits+;
#pragma link C++ class AliTPCtrackerParam;
#pragma link C++ class AliTPCkineGrid;
+// points used in new cluster finder
+#pragma link C++ class AliTPCExactPoint+;
+#pragma link C++ class AliTPCTrackerPoint+;
+#pragma link C++ class AliTPCClusterPoint+;
+#pragma link C++ class AliTPCTrackPoint+;
+#pragma link C++ class AliTPCTrackPointRef+;
+
+#pragma link C++ class AliTPCclustererMI+;
+#pragma link C++ class AliTPCclusterMI+;
+#pragma link C++ class AliTPCclusterLMI+;
+#pragma link C++ class AliTPCtrackerMI+;
+#pragma link C++ class AliTPCseed+;
+
+
+
#endif
#include "AliClusters.h"
#include "AliTPCClustersRow.h"
#include "AliTPCClustersArray.h"
-#include "AliTPCClusterFinder.h"
#include "AliTPCcluster.h"
-
+#include "AliH2F.h"
#include "TMinuit.h"
#include "AliTPC.h"
#include "AliTPCv1.h"
-SRCS:= AliTPCClusterFinder.cxx AliTPC.cxx AliTPCv0.cxx AliTPCv1.cxx AliTPCv2.cxx \
+SRCS:= AliTPC.cxx AliTPCv0.cxx AliTPCv1.cxx AliTPCv2.cxx \
AliTPCv3.cxx AliDetectorParam.cxx AliTPCParam.cxx \
AliTPCParamSR.cxx AliTPCParamCR.cxx AliTPCPRF2D.cxx \
AliTPCRF1D.cxx \
AliTPCtrack.cxx AliTPCtracker.cxx \
AliTPCTrackHits.cxx\
AliTPCDigitizer.cxx\
- AliTPCTrackHitsV2.cxx AliTPCtrackerParam.cxx AliTPCkineGrid.cxx
-
+ AliTPCTrackHitsV2.cxx AliTPCtrackerParam.cxx AliTPCkineGrid.cxx \
+ AliTPCclustererMI.cxx AliTPCclusterMI.cxx AliTPCtrackerMI.cxx
HDRS:= $(SRCS:.cxx=.h)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff