TFile *tf=TFile::Open("AliITStracksV2.root");
if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
TObjArray tarray(2000);
- TTree *tracktree=(TTree*)tf->Get("TreeT");
+ TTree *tracktree=(TTree*)tf->Get("ITSf");
+ if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
TBranch *tbranch=tracktree->GetBranch("tracks");
Int_t nentr=(Int_t)tracktree->GetEntries(),i;
for (i=0; i<nentr; i++) {
TH1F *hz=new TH1F("hz","Longitudinal impact parameter",30,-300,300);
//hmpt->SetFillColor(6);
- TH1F *hgood=new TH1F("hgood","Good tracks",30,0.1,6.1);
- TH1F *hfound=new TH1F("hfound","Found tracks",30,0.1,6.1);
- TH1F *hfake=new TH1F("hfake","Fake tracks",30,0.1,6.1);
- TH1F *hg=new TH1F("hg","",30,0.1,6.1); //efficiency for good tracks
+ AliITStrackV2 *trk=(AliITStrackV2*)tarray.UncheckedAt(0);
+ Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
+ Double_t pmax=6.0+pmin;
+
+ TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
+ TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
+ TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
+ TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
hg->SetLineColor(4); hg->SetLineWidth(2);
- TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.1,6.1);
+ TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax);
hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
- TH1F *hptw=new TH1F("hptw","Weghted pt",30,0.1,6.1);
+ TH1F *hptw=new TH1F("hptw","Weghted pt",30,pmax,pmin);
while (ngood--) {
Int_t lab=gt[ngood].lab, tlab=-1;
Double_t pxg=gt[ngood].px, pyg=gt[ngood].py, pzg=gt[ngood].pz;
Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ if (ptg<pmin) continue;
+
hgood->Fill(ptg);
AliITStrackV2 *track=0;
hg->SetXTitle("Pt (GeV/c)");
hg->Draw();
- TLine *line1 = new TLine(0,1.0,7,1.0); line1->SetLineStyle(4);
+ TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
line1->Draw("same");
- TLine *line2 = new TLine(0,0.9,7,0.9); line2->SetLineStyle(4);
+ TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
line2->Draw("same");
hf->SetFillColor(1);
Int_t entr=(Int_t)cTree->GetEntries();
for (k=0; k<entr; k++) {
if (!cTree->GetEvent(k)) continue;
- Int_t lay,lad,det; geom->GetModuleId(k-1,lay,lad,det);
+ Int_t lay,lad,det; geom->GetModuleId(k,lay,lad,det);
if (lay<1 || lay>6) {
cerr<<"wrong layer !\n"; exit(10);
}
SetLabel(t.GetLabel());
SetChi2(0.);
SetNumberOfClusters(0);
+ //SetConvConst(t.GetConvConst());
+
fdEdx = 0.;
fAlpha = t.GetAlpha();
if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
if (fC22<=0) {cout<<"fC22="<<fC22<<endl; return 0;}
if (fC33<=0) {cout<<"fC33="<<fC33<<endl; return 0;}
if (fC44<=0) {cout<<"fC44="<<fC44<<endl; return 0;}
-
- TMatrixD m(5,5);
- m(0,0)=fC00;
- m(1,0)=fC10; m(1,1)=fC11;
- m(2,0)=fC20; m(2,1)=fC21; m(2,2)=fC22;
- m(3,0)=fC30; m(3,1)=fC31; m(3,2)=fC32; m(3,3)=fC33;
- m(4,0)=fC40; m(4,1)=fC41; m(4,2)=fC42; m(4,3)=fC43; m(4,4)=fC44;
-
- m(0,1)=m(1,0);
- m(0,2)=m(2,0); m(1,2)=m(2,1);
- m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
- m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
- /*
- Double_t det=m.Determinant();
-
- if (det <= 0) {
- cout<<" bad determinant "<<det<<endl;
- m.Print();
- return 0;
- }
- */
- return 1;
+ /*
+ TMatrixD m(5,5);
+ m(0,0)=fC00;
+ m(1,0)=fC10; m(1,1)=fC11;
+ m(2,0)=fC20; m(2,1)=fC21; m(2,2)=fC22;
+ m(3,0)=fC30; m(3,1)=fC31; m(3,2)=fC32; m(3,3)=fC33;
+ m(4,0)=fC40; m(4,1)=fC41; m(4,2)=fC42; m(4,3)=fC43; m(4,4)=fC44;
+
+ m(0,1)=m(1,0);
+ m(0,2)=m(2,0); m(1,2)=m(2,1);
+ m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
+ m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
+
+ Double_t det=m.Determinant();
+
+ if (det <= 0) {
+ cout<<" bad determinant "<<det<<endl;
+ m.Print();
+ return 0;
+ }
+ */
+ return 1;
}
//____________________________________________________________________________
}
*/
+void AliITStrackV2::ResetCovariance() {
+ //------------------------------------------------------------------
+ //This function makes a track forget its history :)
+ //------------------------------------------------------------------
+
+ fC00*=10.;
+ fC10=0.; fC11*=10.;
+ fC20=0.; fC21=0.; fC22*=10.;
+ fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
+ fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
+
+}
Int_t nentr=(Int_t)cTree->GetEntries();
for (i=0; i<nentr; i++) {
if (!cTree->GetEvent(i)) continue;
- Int_t lay,lad,det; g->GetModuleId(i-1,lay,lad,det);
+ //Int_t lay,lad,det; g->GetModuleId(i-1,lay,lad,det);
+ Int_t lay,lad,det; g->GetModuleId(i,lay,lad,det);
Int_t ncl=clusters->GetEntriesFast();
while (ncl--) {
AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
#endif
- fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
+ AliITSdetector &det=fLayers[lay-1].GetDetector(c->GetDetectorIndex());
+ Double_t r=det.GetR(); r=TMath::Sqrt(r*r + c->GetY()*c->GetY());
+ if (r > TMath::Abs(c->GetZ())-0.5)
+ fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
}
clusters->Delete();
}
fI=kMaxLayer;
}
+#ifdef DEBUG
static Int_t lbl;
+#endif
Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
//--------------------------------------------------------------------
}
in->cd();
- TTree *tpcTree=(TTree*)in->Get("TreeT");
+ TTree *tpcTree=(TTree*)in->Get("TPCf");
if (!tpcTree) {
cerr<<"AliITStrackerV2::Clusters2Tracks() ";
cerr<<"can't get a tree with TPC tracks !\n";
tpcTree->SetBranchAddress("tracks",&itrack);
out->cd();
- TTree itsTree("TreeT","Tree with ITS tracks");
+ TTree itsTree("ITSf","Tree with ITS tracks");
AliITStrackV2 *otrack=&fBestTrack;
itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
if (!tpcTree->GetEvent(i)) continue;
Int_t tpcLabel=itrack->GetLabel(); //save the TPC track label
-lbl=tpcLabel;
#ifdef DEBUG
+lbl=tpcLabel;
if (TMath::Abs(tpcLabel)!=LAB) continue;
cout<<tpcLabel<<" *****************\n";
#endif
Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
fTrackToFollow.Improve(x0,fYV,fZV);
- //Double_t xk=77.2;
- Double_t xk=88.;
+ Double_t xk=80.;
fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
xk-=0.005;
xk-=0.005;
fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk-=14.5;
+ xk=61.;
fTrackToFollow.PropagateTo(xk,0.,0.); //C02
xk -=0.005;
return 0;
}
+Int_t AliITStrackerV2::PropagateBack(const TFile *inp, TFile *out) {
+ //--------------------------------------------------------------------
+ //This functions propagates reconstructed ITS tracks back
+ //--------------------------------------------------------------------
+ TFile *in=(TFile*)inp;
+ TDirectory *savedir=gDirectory;
+
+ if (!in->IsOpen()) {
+ cerr<<"AliITStrackerV2::PropagateBack(): ";
+ cerr<<"file with ITS tracks is not open !\n";
+ return 1;
+ }
+
+ if (!out->IsOpen()) {
+ cerr<<"AliITStrackerV2::PropagateBack(): ";
+ cerr<<"file for back propagated ITS tracks is not open !\n";
+ return 2;
+ }
+
+ in->cd();
+ TTree *itsTree=(TTree*)in->Get("ITSf");
+ if (!itsTree) {
+ cerr<<"AliITStrackerV2::PropagateBack() ";
+ cerr<<"can't get a tree with ITS tracks !\n";
+ return 3;
+ }
+ AliITStrackV2 *itrack=new AliITStrackV2;
+ itsTree->SetBranchAddress("tracks",&itrack);
+
+ out->cd();
+ TTree backTree("ITSb","Tree with back propagated ITS tracks");
+ AliTPCtrack *otrack=0;
+ backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
+
+ Int_t ntrk=0;
+
+ Int_t nentr=(Int_t)itsTree->GetEntries();
+ for (Int_t i=0; i<nentr; i++) {
+ itsTree->GetEvent(i);
+ ResetTrackToFollow(*itrack);
+ fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
+ Int_t itsLabel=fTrackToFollow.GetLabel(); //save the ITS track label
+
+#ifdef DEBUG
+if (TMath::Abs(itsLabel)!=LAB) continue;
+cout<<itsLabel<<" *****************\n";
+#endif
+
+ try {
+ Int_t nc=itrack->GetNumberOfClusters();
+#ifdef DEBUG
+for (Int_t k=0; k<nc; k++) {
+ Int_t index=itrack->GetClusterIndex(k);
+ AliITSclusterV2 *c=(AliITSclusterV2*)GetCluster(index);
+ cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
+}
+#endif
+ Int_t idx=0, l=0;
+ const AliITSclusterV2 *c=0;
+ if (nc--) {
+ idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
+ c=(AliITSclusterV2*)GetCluster(idx);
+ }
+ for (fI=0; fI<kMaxLayer; fI++) {
+ AliITSlayer &layer=fLayers[fI];
+ Double_t r=layer.GetR();
+ if (fI==2 || fI==4) {
+ Double_t rs=0.5*(fLayers[fI-1].GetR() + r);
+ Double_t ds=0.034; if (fI==4) ds=0.039;
+ Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
+ if (!fTrackToFollow.PropagateTo(rs,1*dx0r,dr)) throw "";
+ }
+
+ Double_t x,y,z;
+ if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z))
+ throw "AliITStrackerV2::PropagateBack: failed to estimate track !";
+ Double_t phi=TMath::ATan2(y,x);
+ Double_t d=layer.GetThickness(phi,z);
+ Int_t idet=layer.FindDetectorIndex(phi,z);
+ if (idet<0)
+ throw "AliITStrackerV2::PropagateBack: failed to find a detector !\n";
+ const AliITSdetector &det=layer.GetDetector(idet);
+ r=det.GetR(); phi=det.GetPhi();
+ if (!fTrackToFollow.Propagate(phi,r,d/21.82*2.33,d*2.33)) throw "";
+
+ const AliITSclusterV2 *cl=0;
+ Int_t index=0;
+ Double_t maxchi2=kMaxChi2;
+
+ if (l==fI) {
+ if (idet != c->GetDetectorIndex()) {
+ idet=c->GetDetectorIndex();
+ const AliITSdetector &det=layer.GetDetector(idet);
+ r=det.GetR(); phi=det.GetPhi();
+ if (!fTrackToFollow.Propagate(phi,r,0.,0.)) throw "";
+ }
+ Double_t chi2=fTrackToFollow.GetPredictedChi2(c);
+ if (chi2<kMaxChi2) {
+ cl=c; maxchi2=chi2; index=idx;
+ }
+ if (nc--) {
+ idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
+ c=(AliITSclusterV2*)GetCluster(idx);
+ }
+ }
+
+ if (fTrackToFollow.GetNumberOfClusters()>2) {
+ Double_t dz=3*TMath::Sqrt(fTrackToFollow.GetSigmaZ2()+kSigmaZ2[fI]);
+ Double_t dy=3*TMath::Sqrt(fTrackToFollow.GetSigmaY2()+kSigmaY2[fI]);
+ Double_t zmin=fTrackToFollow.GetZ() - dz;
+ Double_t zmax=fTrackToFollow.GetZ() + dz;
+ Double_t ymin=fTrackToFollow.GetY() + phi*r - dy;
+ Double_t ymax=fTrackToFollow.GetY() + phi*r + dy;
+ layer.SelectClusters(zmin,zmax,ymin,ymax);
+
+ const AliITSclusterV2 *cc=0; Int_t ci;
+ while ((cc=layer.GetNextCluster(ci))!=0) {
+ if (idet != cc->GetDetectorIndex()) continue;
+ Double_t chi2=fTrackToFollow.GetPredictedChi2(cc);
+ if (chi2<maxchi2) {
+ cl=cc; index=(fI<<28)+ci; maxchi2=chi2;
+ }
+ }
+ }
+
+ if (cl) {
+ if (!fTrackToFollow.Update(cl,maxchi2,index))
+ cerr<<"AliITStrackerV2::PropagateBack: filtering failed !\n";
+ continue;
+ }
+ }
+
+ Double_t xk=61.;
+ fTrackToFollow.PropagateTo(xk,0.,0.); //Air
+
+ xk +=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
+ xk +=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk +=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk +=0.5;
+ fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex
+ xk +=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk +=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk +=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
+
+ xk=80.;
+ fTrackToFollow.PropagateTo(xk,0.,0.); //CO2
+
+ xk+=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ xk+=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk+=2.0;
+ fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
+ xk+=0.02;
+ fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
+ xk+=0.005;
+ fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+
+ fTrackToFollow.SetLabel(itsLabel);
+ otrack=new AliTPCtrack(fTrackToFollow,fTrackToFollow.GetAlpha());
+ backTree.Fill(); delete otrack;
+ UseClusters(&fTrackToFollow);
+ cerr << ++ntrk << " \r";
+ }
+ catch (const Char_t *msg) {
+ cerr<<msg<<endl;
+ }
+ }
+
+ backTree.Write();
+ savedir->cd();
+ cerr<<"Number of ITS tracks: "<<nentr<<endl;
+ cerr<<"Number of back propagated ITS tracks: "<<ntrk<<endl;
+
+ delete itrack;
+
+ return 0;
+}
AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
//--------------------------------------------------------------------
while (fI) {
fI--;
-
#ifdef DEBUG
cout<<fI<<' ';
#endif
-
AliITSlayer &layer=fLayers[fI];
AliITStrackV2 &track=fTracks[fI];
+ Double_t r=layer.GetR();
if (fI==3 || fI==1) {
- Double_t rs=0.5*(fLayers[fI+1].GetR() + layer.GetR());
+ Double_t rs=0.5*(fLayers[fI+1].GetR() + r);
Double_t ds=0.034; if (fI==3) ds=0.039;
Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,1*dx0r,dr);
//find intersection
Double_t x,y,z;
- if (!fTrackToFollow.GetGlobalXYZat(layer.GetR(),x,y,z)) {
+ if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
cerr<<"AliITStrackerV2::FollowProlongation: failed to estimate track !\n";
break;
}
//propagate to the intersection
const AliITSdetector &det=layer.GetDetector(idet);
- if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR(),1*d/21.82*2.33,d*2.33))
- {
+ phi=det.GetPhi();
+ if (!fTrackToFollow.Propagate(phi,det.GetR(),1*d/21.82*2.33,d*2.33)) {
cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
break;
}
//Select possible prolongations and store the current track estimation
track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[fI]);
- if (dz<0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
+ if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
if (dz > kMaxRoad/4) {
//cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
break;
}
Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]);
- if (dy<0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
+ if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
if (dy > kMaxRoad/4) {
//cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
break;
Double_t zmin=track.GetZ() - dz;
Double_t zmax=track.GetZ() + dz;
- Double_t ymin=track.GetY() + layer.GetR()*det.GetPhi() - dy;
- Double_t ymax=track.GetY() + layer.GetR()*det.GetPhi() + dy;
- if (ymax>layer.GetR()*2*TMath::Pi()) {
- ymax-=layer.GetR()*2*TMath::Pi();
- ymin-=layer.GetR()*2*TMath::Pi();
- }
+ Double_t ymin=track.GetY() + r*phi - dy;
+ Double_t ymax=track.GetY() + r*phi + dy;
layer.SelectClusters(zmin,zmax,ymin,ymax);
- //if (1/TMath::Abs(track.Get1Pt())<0.2)
- //cout<<fI<<' '<<1/TMath::Abs(track.Get1Pt())<<' '
- // <<dy<<' '<<dz<<' '<<layer.InRoad()<<endl;
-
//take another prolongation
if (!TakeNextProlongation()) if (!tryAgain--) break;
tryAgain=kLayersToSkip;
const AliITSclusterV2 *c=0; Int_t ci=-1;
Double_t chi2=12345.;
while ((c=layer.GetNextCluster(ci))!=0) {
- //if (fI==0)
- //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; //88888888888888888888
+
+#ifdef DEBUG
+//if (fI==0)
+//if (c->GetLabel(0)!=TMath::Abs(lbl)) continue;
+#endif
+
Int_t idet=c->GetDetectorIndex();
if (t.GetDetectorIndex()!=idet) {
// This function sets the "window"
//--------------------------------------------------------------------
fI=FindClusterIndex(zmin); fZmax=zmax;
+ Double_t circle=2*TMath::Pi()*fR;
+ if (ymax>circle) { ymax-=circle; ymin-=circle; }
fYmin=ymin; fYmax=ymax;
}