From 59dfcaddf72fc4d0990e82f74c88a3afa34f3ebd Mon Sep 17 00:00:00 2001 From: cblume Date: Mon, 10 Jun 2002 14:56:40 +0000 Subject: [PATCH] Add new tracking macros --- TRD/AliTRDbackTrackAnalysis.C | 606 ++++++++++++++++++++++++++++++++ TRD/AliTRDclusterErrors.C | 268 ++++++++++++++ TRD/AliTRDsaveMCtracks.C | 314 +++++++++++++++++ TRD/AliTRDtoTOFanalysis.C | 185 ++++++++++ TRD/AliTRDtrackReconstruction.C | 112 ++++++ 5 files changed, 1485 insertions(+) create mode 100644 TRD/AliTRDbackTrackAnalysis.C create mode 100644 TRD/AliTRDclusterErrors.C create mode 100644 TRD/AliTRDsaveMCtracks.C create mode 100644 TRD/AliTRDtoTOFanalysis.C create mode 100644 TRD/AliTRDtrackReconstruction.C diff --git a/TRD/AliTRDbackTrackAnalysis.C b/TRD/AliTRDbackTrackAnalysis.C new file mode 100644 index 00000000000..c1b82d385a0 --- /dev/null +++ b/TRD/AliTRDbackTrackAnalysis.C @@ -0,0 +1,606 @@ +#ifndef __CINT__ + #include + #include "AliTRDtracker.h" + #include "AliTRDcluster.h" + #include "AliTRDv1.h" + #include "AliTRDgeometry.h" + #include "AliTRDparameter.h" + #include "alles.h" + #include "AliTRDmcTrack.h" + #include "AliTRDtrack.h" + + #include "TFile.h" + #include "TParticle.h" + #include "TStopwatch.h" +#endif + +void AliTRDbackTrackAnalysis() { + + // const Int_t nPrimaries = 84210/400; + const Int_t nPrimaries = 84210/16; + + + Float_t Pt_min = 0.; + Float_t Pt_max = 20.; + + TH1F *hp=new TH1F("hp","PHI resolution",100,-20.,20.); hp->SetFillColor(4); + TH1F *hl=new TH1F("hl","LAMBDA resolution",100,-100,100);hl->SetFillColor(4); + TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.); + + TH1F *hpta=new TH1F("hpta","norm. Pt resolution",50,-5.,5.); + hpta->SetFillColor(17); + hpta->SetXTitle("(%)"); + + TH1F *hla=new TH1F("hla","norm. Lambda resolution",50,-5.,5.); + hla->SetFillColor(17); + TH1F *hya=new TH1F("hya","norm. Y resolution",50,-5.,5.); + hya->SetFillColor(17); + TH1F *hza=new TH1F("hza","norm. Z resolution",50,-5.,5.); + hza->SetFillColor(17); + + hpt->SetFillColor(2); + + TH1F *hy=new TH1F("hy","Y resolution",50,-0.5,0.5);hy->SetLineColor(4); + hy->SetLineWidth(2); + hy->SetXTitle("(cm)"); + + TH1F *hz=new TH1F("hz","Z resolution",80,-4,4);hz->SetLineColor(2); + hz->SetLineWidth(2); + hz->SetXTitle("(cm)"); + + TH1F *hx=new TH1F("hx","X(out)",150,250.,400.); hx->SetFillColor(4); + + const Int_t nPtSteps = 30; + const Float_t maxPt = 3.; + + TH1F *hgood=new TH1F("hgood","long TRD tracks, TPC seeds",nPtSteps,0,maxPt); + hgood->SetYTitle("Counts"); + hgood->SetXTitle("Pt (GeV/c)"); + hgood->SetLineColor(3); + + TH1F *hGoodAndSeed = new TH1F("GoodAndSeed","TPC seed and good",nPtSteps,0,maxPt); + hGoodAndSeed->SetLineColor(2); + + TH2F *hRTvsMC = new TH2F("RTvsMC","RTvsMC",100,-4.5,95.5,100,-4.5,95.5); + hRTvsMC->SetMarkerColor(4); + hRTvsMC->SetMarkerSize(2); + + TH2F *hXvsMCX = new TH2F("XvsMCX","XvsMCX",150,250,400,150,250,400); + hXvsMCX->SetMarkerColor(4); + hXvsMCX->SetMarkerSize(2); + + TH2F *h2seed = new TH2F("2dGood","TPC seeds",60,0.,60,50,0.,3.); + h2seed->SetMarkerColor(4); + h2seed->SetMarkerSize(2); + + TH2F *h2lost = new TH2F("2dSeedAndGood","SeedButNotGood",60,0.,60.,50,0.,3.); + h2lost->SetMarkerColor(2); + h2lost->SetMarkerSize(2); + + + TH1F *hseed=new TH1F("seed","TPC seeds",nPtSteps,0,maxPt); + hseed->SetLineColor(4); + + + TH1F *hfound=new TH1F("hfound","Found tracks",nPtSteps,0,maxPt); + hfound->SetYTitle("Counts"); + hfound->SetXTitle("Pt (GeV/c)"); + + TH1F *heff=new TH1F("heff","Matching Efficiency",nPtSteps,0,maxPt); // efficiency, good tracks + heff->SetLineColor(4); heff->SetLineWidth(2); + heff->SetXTitle("Pt, GeV/c"); + heff->SetYTitle("Efficiency"); + + TH1F *hSeedEff = new TH1F("hSeedEff","TPC Efficiency",nPtSteps,0,maxPt); + hSeedEff->SetLineColor(4); hSeedEff->SetLineWidth(2); + hSeedEff->SetXTitle("Pt, GeV/c"); + hSeedEff->SetYTitle("Efficiency"); + + + TH1F *hol=new TH1F("hol","Overlap fraction",105,-2.5,102.5); + hol->SetLineColor(4); hol->SetLineWidth(2); + hol->SetXTitle("Fraction,(%)"); + hol->SetYTitle("Counts"); + + TH1F *hend=new TH1F("end","missing tail",80,-10.5,69.5); + + TH1F *hFraction=new TH1F("fraction","Fraction of found clusters",110,0,1.1); + TH1F *hCorrect=new TH1F("correct","Fraction of correct clusters",110,0,1.1); + + + Int_t nEvent = 0; + const Int_t maxIndex = nPrimaries; + Bool_t seedLabel[maxIndex]; + Int_t mcIndex[maxIndex]; + Int_t rtIndex[maxIndex]; + + for(Int_t i = 0; i < maxIndex; i++) { + seedLabel[i] = kFALSE; + mcIndex[i] = -1; + rtIndex[i] = -1; + } + + // mark available seeds from TPC + Int_t nSeeds = 0; + Int_t nPrimarySeeds = 0; + + printf("marking found seeds from TPC\n"); + TDirectory *savedir=gDirectory; + + TFile *in=TFile::Open("AliTPCBackTracks.root"); + if (!in->IsOpen()) { + cerr<<"can't open file AliTPCBackTracks.root !\n"; return; + } + + char tname[100]; + sprintf(tname,"seedsTPCtoTRD_%d",nEvent); + TTree *seedTree=(TTree*)in->Get(tname); + if (!seedTree) { + cerr<<"AliTRDtracker::PropagateBack(): "; + cerr<<"can't get a tree with seeds from TPC !\n"; + } + + AliTPCtrack *seed=new AliTPCtrack; + seedTree->SetBranchAddress("tracks",&seed); + + Int_t n=(Int_t)seedTree->GetEntries(); + for (Int_t i=0; iGetEvent(i); + Int_t lbl = seed->GetLabel(); + if(lbl < 0) { + printf("negative seed label %d \n",lbl); + continue; + } + if(lbl >= maxIndex) continue; + seedLabel[lbl] = kTRUE; + if(lbl < nPrimaries) nPrimarySeeds++; + nSeeds++; + } + delete seed; + delete seedTree; + + printf("Found %d seeds from primaries among overall %d seeds \n", + nPrimarySeeds, nSeeds); + + savedir->cd(); + // done with marking TPC seeds + + + TFile *tf=TFile::Open("AliTRDtracks.root"); + + if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;} + TObjArray tarray(2000); + + sprintf(tname,"TRDb_%d",nEvent); + TTree *tracktree=(TTree*)tf->Get(tname); + + TBranch *tbranch=tracktree->GetBranch("tracks"); + + Int_t nRecTracks = (Int_t) tracktree->GetEntries(); + cerr<<"Found "<SetAddress(&iotrack); + tracktree->GetEvent(i); + tarray.AddLast(iotrack); + Int_t trackLabel = iotrack->GetLabel(); + + // printf("rt with %d clusters and label %d \n", + // iotrack->GetNumberOfClusters(), trackLabel); + + if(trackLabel < 0) continue; + if(trackLabel >= maxIndex) continue; + rtIndex[trackLabel] = i; + } + tf->Close(); + + // return; + + // Load MC tracks + TFile *mctf=TFile::Open("AliTRDmcTracks.root"); + if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDmcTracks.root !\n"; return;} + TObjArray mctarray(2000); + TTree *mctracktree=(TTree*)mctf->Get("MCtracks"); + TBranch *mctbranch=mctracktree->GetBranch("MCtracks"); + Int_t nMCtracks = (Int_t) mctracktree->GetEntries(); + cerr<<"Found "<SetAddress(&ioMCtrack); + mctracktree->GetEvent(i); + mctarray.AddLast(ioMCtrack); + Int_t mcLabel = ioMCtrack->GetTrackIndex(); + if(mcLabel < 0) {printf("negative mc label detected!\n"); continue;} + if(mcLabel >= maxIndex) continue; + mcIndex[mcLabel] = i; + } + mctf->Close(); + + + // Load clusters + + TFile *geofile =TFile::Open("AliTRDclusters.root"); + AliTRDtracker *Tracker = new AliTRDtracker(geofile); + Tracker->SetEventNumber(nEvent); + + AliTRDgeometry *fGeom = (AliTRDgeometry*) geofile->Get("TRDgeometry"); + AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter"); + + Char_t *alifile = "AliTRDclusters.root"; + TObjArray carray(2000); + TObjArray *ClustersArray = &carray; + Tracker->ReadClusters(ClustersArray,alifile); + + + // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits + alifile = "galice.root"; + + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + + // Define the objects + AliTRDv1 *TRD = (AliTRDv1*) gAlice->GetDetector("TRD"); + + // Import the Trees for the event nEvent in the file + const Int_t nparticles = gAlice->GetEvent(nEvent); + if (nparticles <= 0) return; + + TParticle *p; + Bool_t electrons[300000] = { kFALSE }; + + Bool_t mark_electrons = kFALSE; + if(mark_electrons) { + printf("mark electrons\n"); + + for(Int_t i = nPrimaries; i < nparticles; i++) { + p = gAlice->Particle(i); + if(p->GetMass() > 0.01) continue; + if(p->GetMass() < 0.00001) continue; + electrons[i] = kTRUE; + } + } + + AliTRDcluster *cl = 0; + Int_t nw, label, index, ti[3], tbwc; + Int_t det, plane, ltb, gtb, gap, max_gap, sector, mc_sector, min_tb, max_tb; + Double_t Pt, Px, Py, Pz; + Double_t mcPt, mcPx, mcPy, mcPz; + Double_t x,y,z, mcX; + Int_t rtClusters, rtCorrect; + + printf("\n"); + + AliTRDmcTrack *mct = 0; + AliTRDtrack *rt = 0; + + Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region + Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region + Double_t dx = (Double_t) fPar->GetTimeBinSize(); + + Int_t tbAmp = fPar->GetTimeBefore(); + Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); + Int_t tbDrift = fPar->GetTimeMax(); + Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); + + tbDrift = TMath::Min(tbDrift,maxDrift); + tbAmp = TMath::Min(tbAmp,maxAmp); + + const Int_t nPlanes = fGeom->Nplan(); + const Int_t tbpp = Tracker->GetTimeBinsPerPlane(); + const Int_t nTB = tbpp * nPlanes; + Int_t mask[nTB]; + + for(Int_t i=0; i < maxIndex; i++) { + + rt = 0; mct = 0; + if(rtIndex[i] >= 0) rt = (AliTRDtrack *) tarray.UncheckedAt(rtIndex[i]); + if(mcIndex[i] >= 0) mct = (AliTRDmcTrack *) mctarray.UncheckedAt(mcIndex[i]); + + if(!mct) continue; + label = mct->GetTrackIndex(); + + // if(TMath::Abs(mct->GetMass()-0.136) > 0.01) continue; + + Int_t ncl = mct->GetNumberOfClusters(); + + // check how many time bins have a cluster + + for(nw = 0; nw < nTB; nw++) mask[nw] = -1; + + for(Int_t j = 0; j < ncl; j++) { + index = mct->GetClusterIndex(j); + cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index); + + for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw); + + if((ti[0] != label) && (ti[1] != label) && (ti[2] != label)) { + printf("wrong track label: %d, %d, %d != %d \n", + ti[0], ti[1], ti[2], label); + } + det=cl->GetDetector(); + plane = fGeom->GetPlane(det); + ltb = cl->GetLocalTimeBin(); + gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); + mask[gtb] = index; + } + + + + + for(plane = 0; plane < nPlanes; plane++) { + for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) { + gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); + if(mask[gtb] > -1) break; + } + if(ltb >= -tbAmp) break; + } + if((plane == nPlanes) && (ltb == -tbAmp-1)) { + printf("warning: for track %d min tb is not found and set to %d!\n", + label, nTB-1); + min_tb = nTB-1; + // for(Int_t tb = 0; tbGetGlobalTimeBin(0,plane,ltb); + } + + for(plane = nPlanes-1 ; plane>=0; plane--) { + for(ltb = -tbAmp; ltb < tbDrift; ltb++) { + gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); + if(mask[gtb] > -1) break; + } + if(ltb < tbDrift) break; + } + if((plane == -1) && (ltb == tbDrift)) { + printf("warning: for track %d max tb is not found and set to 0!\n",label); + // for(Int_t tb = 0; tbGetX(0,0,tbDrift-1); + } + else { + max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb); + mcX = Tracker->GetX(0,plane,ltb); + cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb]); + det=cl->GetDetector(); + sector = fGeom->GetSector(det); + mc_sector = sector; + } + + + tbwc = 0; + max_gap = 0; + gap = -1; + + for(nw = min_tb; nw < max_tb+1; nw++) { + gap++; + if(mask[nw] > -1) { + tbwc++; + if(gap > max_gap) max_gap = gap; + gap = 0; + } + } + + if(tbwc < ((Int_t) (nTB * Tracker->GetMinClustersInTrack()))) continue; + if(gap > Tracker->GetMaxGap()) continue; + p = gAlice->Particle(label); + + printf("good track %d with min_tb %d, max_tb %d, gap %d\n", + label,min_tb,max_tb,gap); + + hgood->Fill(p->Pt()); + + if(!rt) continue; + + if(rt->GetLabel() != label) { + printf("mct vs rt index mismatch: %d != %d \n", + label, rt->GetLabel()); + return; + } + + if(!seedLabel[i]) continue; + + hGoodAndSeed->Fill(p->Pt()); + + hXvsMCX->Fill(mcX,rt->GetX()); + + rtClusters = rt->GetNumberOfClusters(); + + // find number of tb with correct cluster + rtCorrect = 0; + + for(Int_t j = 0; j < rtClusters; j++) { + index = rt->GetClusterIndex(j); + cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index); + + for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw); + + if((ti[0] != label)&&(ti[1] != label)&&(ti[2] != label)) continue; + rtCorrect++; + } + + Float_t foundFraction = ((Float_t) rtCorrect) / (tbwc + 0.00001); + Float_t correctFraction = ((Float_t) rtCorrect) / ((Float_t) rtClusters); + + if(foundFraction > 1) printf("fraction = %f/%f \n", + (Float_t) rtCorrect, (Float_t) tbwc); + + + if(foundFraction <= 1) hFraction->Fill(foundFraction); + hCorrect->Fill(correctFraction); + + hRTvsMC->Fill((Float_t) tbwc, (Float_t) rtClusters); + + if((foundFraction < 0.7) || (correctFraction < 0.7)) { + printf("not found track %d with FrctnFound %f and FrctnCorrect %f\n", + label, foundFraction, correctFraction); + continue; + } + + hfound->Fill(p->Pt()); + + Pt = TMath::Abs(rt->GetPt()); + Double_t cc = TMath::Abs(rt->GetSigmaC2()); + mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1); + rt->GetPxPyPz(Px,Py,Pz); + + printf("\n\ntrack %d \n", label); + printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); + printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); + + mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy); + if(mcPt < 0.0001) mcPt = 0.0001; + + Float_t lamg=TMath::ATan2(mcPz,mcPt); + Float_t lam =TMath::ATan2(Pz,Pt); + if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt); + else hpta->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc)); + + Float_t phig=TMath::ATan2(mcPy,mcPx); + Float_t phi =TMath::ATan2(Py,Px); + + + if(!(rt->PropagateTo(x))) continue; + + + if((mcPt > Pt_min) && (mcPt < Pt_max)) { + hl->Fill((lam - lamg)*1000.); + hla->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2())); + if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt); + else hpt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); + hp->Fill((phi - phig)*1000.); + hy->Fill(rt->GetY() - y); + hz->Fill(rt->GetZ() - z); + hya->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2())); + hza->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2())); + hx->Fill((Float_t) x); + } + } + + + gStyle->SetOptStat(0); + // gStyle->SetOptStat(111110); + gStyle->SetOptFit(1); + + + /* + TCanvas *c1=new TCanvas("c1","",0,0,700,850); + + // gStyle->SetOptStat(0); + // gStyle->SetOptStat(111110); + gStyle->SetOptFit(1); + + TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw(); + p1->cd(); p1->SetFillColor(10); p1->SetFrameFillColor(10); + hgood->Draw(); hGoodAndSeed->Draw("same"); c1->cd(); + + TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw(); + p2->cd(); p2->SetFillColor(10); p2->SetFrameFillColor(10); + hFraction->SetYTitle("Counts"); + hFraction->SetXTitle("Fraction"); + hFraction->Draw(); c1->cd(); + + TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw(); + p3->cd(); p3->SetFillColor(10); p3->SetFrameFillColor(10); + hfound->Draw(); c1->cd(); + + TPad *p4=new TPad("p4","",0.5,0,1.,0.3); p4->Draw(); + p4->cd(); p4->SetFillColor(10); p4->SetFrameFillColor(10); + hCorrect->SetYTitle("Counts"); + hCorrect->SetXTitle("Fraction"); + hCorrect->Draw(); c1->cd(); + + TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd(); + p5->SetFillColor(10); p5->SetFrameFillColor(10); + hgood->Sumw2(); hGoodAndSeed->Sumw2(); // hfake->Sumw2(); + hSeedEff->Divide(hGoodAndSeed,hgood,1,1.,"b"); + // hf->Divide(hfake,hgood,1,1.,"b"); + hSeedEff->SetMaximum(1.4); + hSeedEff->SetYTitle("TPC efficiency"); + hSeedEff->SetXTitle("Pt (GeV/c)"); + hSeedEff->Draw(); + + TLine *line1 = new TLine(0,1.0,maxPt,1.0); line1->SetLineStyle(4); + line1->Draw("same"); + TLine *line2 = new TLine(0,0.9,maxPt,0.9); line2->SetLineStyle(4); + line2->Draw("same"); + */ + + TCanvas *c2=new TCanvas("c2","",0,0,700,850); + + TPad *p12=new TPad("p12","",0,0.3,.5,.6); p12->Draw(); + p12->cd(); p12->SetFillColor(42); p12->SetFrameFillColor(10); + hp->SetFillColor(4); hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c2->cd(); + + TPad *p22=new TPad("p22","",0.5,.3,1,.6); p22->Draw(); + p22->cd(); p22->SetFillColor(42); p22->SetFrameFillColor(10); + hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c2->cd(); + + TPad *p32=new TPad("p32","",0,0,0.5,0.3); p32->Draw(); + p32->cd(); p32->SetFillColor(42); p32->SetFrameFillColor(10); + hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c2->cd(); + + TPad *p42=new TPad("p42","",0.5,0,1.,0.3); p42->Draw(); + p42->cd(); p42->SetFillColor(42); p42->SetFrameFillColor(10); + hgood->Draw(); hGoodAndSeed->Draw("same"); c2->cd(); + + TPad *p52=new TPad("p52","",0,0.6,1,1); p52->Draw(); p52->cd(); + p52->SetFillColor(41); p52->SetFrameFillColor(10); + hfound->Sumw2(); + heff->Divide(hfound,hGoodAndSeed,1,1.,"b"); + // hf->Divide(hfake,hgood,1,1.,"b"); + heff->SetMaximum(1.4); + heff->SetXTitle("Pt (GeV/c)"); + heff->Draw(); + + TLine *line12 = new TLine(0,1.0,maxPt,1.0); line12->SetLineStyle(4); + line12->Draw("same"); + TLine *line22 = new TLine(0,0.9,maxPt,0.9); line22->SetLineStyle(4); + line22->Draw("same"); + + c2->Print("matching.ps","ps"); + + + /* + TCanvas *cxyz = new TCanvas("cxyz","",50,50,750,900); + cxyz->Divide(2,2); + cxyz->cd(1); hx->Draw(); + cxyz->cd(2); hy->Draw(); + cxyz->cd(3); hz->Draw(); + cxyz->cd(4); hXvsMCX->Draw(); + */ + + TCanvas *cs = new TCanvas("cs","",0,0,700,850); + cs->Divide(2,3); + + cs->cd(1); hy->Fit("gaus"); + cs->cd(2); hz->Fit("gaus"); + cs->cd(3); hpta->Fit("gaus"); + cs->cd(4); hla->Fit("gaus"); + cs->cd(5); hya->Fit("gaus"); + cs->cd(6); hza->Fit("gaus"); + + cs->Print("resolution.ps","ps"); + + /* + TCanvas *cvs = new TCanvas("cvs","",0,0,700,850); + cvs->cd(); hRTvsMC->Draw(); + */ + +} + + diff --git a/TRD/AliTRDclusterErrors.C b/TRD/AliTRDclusterErrors.C new file mode 100644 index 00000000000..9a6a5552149 --- /dev/null +++ b/TRD/AliTRDclusterErrors.C @@ -0,0 +1,268 @@ +#ifndef __CINT__ + #include + + #include "AliTRDtracker.h" + #include "AliTRDcluster.h" + #include "AliTRDhit.h" + #include "AliTRDv1.h" + #include "AliTRDgeometry.h" + #include "AliTRDparameter.h" + + #include "alles.h" + #include "TFile.h" + #include "TStopwatch.h" + +#endif + +void AliTRDclusterErrors() { + + Int_t No_of_tracks_to_analyze = 40; + + const Int_t tbpp = 15; + const Int_t nPlanes = 6; + const Int_t ntb = tbpp * nPlanes; + + TH1F *hy = new TH1F("delta R*phi","Cluster displacement in R*phi",200,-1.,1.); + TH1F *hyp = new TH1F("delta R*phi pos","delta R*phi, positive",200,-1.,1.); + TH1F *hym = new TH1F("delta R*phi neg","delta R*phi, negative",200,-1.,1.); + + TH1F *hyn = new TH1F("Norm., d(R*phi)","Norm. cluster displacement in R*phi",400,-8.,8.); + TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.); + TH2F *hy2 = new TH2F("Amp vs delta R*phi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.); + TH2F *herr = new TH2F("sigmaY vs delta R*phi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1); + TH2F *hy3 = new TH2F("Position within pad vs delta R*phi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5); + TH2F *hy4 = new TH2F("local tb vs delta R*phi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5); + + hy->SetXTitle("Displacement, cm"); + hyn->SetXTitle("Displacement, SigmaY"); + hy2->SetXTitle("Displacement, cm"); + hy2->SetYTitle("Amplitude"); + hz->SetXTitle("Displacement, cm"); + hy3->SetXTitle("Displacement, cm"); + hy3->SetYTitle("Position, cm"); + hy4->SetXTitle("Displacement, cm"); + hy4->SetYTitle("local time bin"); + + /* + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + cout << "Loaded shared libraries" << endl; + } + */ + + // Load clusters + Char_t *clusterfile = "AliTRDclusters.root"; + Int_t nEvent = 0; + + TObjArray carray(2000); + TObjArray *ClustersArray = &carray; + TFile *geofile =TFile::Open("AliTRDclusters.root"); + + AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter"); + AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry"); + + AliTRDtracker *Tracker = new AliTRDtracker(geofile); + Tracker->SetEventNumber(nEvent); + Tracker->ReadClusters(ClustersArray,clusterfile); + Int_t nClusters = carray.GetEntriesFast(); + + printf("Total number of clusters %d \n", nClusters); + + Char_t *alifile = "galice.root"; + // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD"); + + // Import the Trees for the event nEvent in the file + Int_t nparticles = gAlice->GetEvent(nEvent); + + TObjArray *particles=gAlice->Particles(); + + TTree *hitTree = gAlice->TreeH(); + Int_t nBytes = 0; + + // Get the number of entries in the hit tree + // (Number of primary particles creating a hit somewhere) + Int_t nTrack = (Int_t) hitTree->GetEntries(); + cout << " Found " << nTrack << " primary particles with hits" << endl; + No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack); + + + // Loop through particles and fill histoes + + Float_t hitY[ntb]; + Float_t hitZ[ntb]; + Float_t hitO[ntb]; + Float_t clusterY[ntb]; + Float_t clusterZ[ntb]; + Float_t clusterQ[ntb]; + Float_t clusterSigmaY[ntb]; + Float_t pos[3]; + Float_t rot[3]; + Float_t global[3]; + Int_t det = 0; + Int_t track_index[3]; + + printf("\n"); + + for (Int_t ii=0; ii -1; ltb--) { + + if(ii >= nTrack) continue; + + TParticle *p = gAlice->Particle(ii); + if (p->GetFirstMother()>=0) continue; + TParticlePDG *pdg = p->GetPDG(); + Float_t charge=pdg->Charge(); + if(TMath::Abs(charge) < 0.5) continue; + + Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); + Double_t x = Tracker->GetX(0,plane,ltb); + + // loop through clusters + + Bool_t cluster_found = kFALSE; + Int_t nw = 0; + + for (Int_t i = 0; i < nClusters; i++) { + + AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i); + + nw = cl->GetDetector(); + + nw = fGeo->GetPlane(nw); + + if(nw != plane) continue; + for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j); + if((track_index[0] != ii) && + (track_index[1] != ii) && + (track_index[2] != ii)) continue; + + nw = cl->GetLocalTimeBin(); + if(nw != ltb) continue; + + clusterY[gtb] = cl->GetY(); + clusterZ[gtb] = cl->GetZ(); + clusterQ[gtb] = cl->GetQ(); + + clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2()); + cluster_found = kTRUE; + break; + } + + if(!cluster_found) continue; + + gAlice->ResetHits(); + + // nBytes += hitTree->GetEvent(nPrimaries - ii - 1); + nBytes += hitTree->GetEvent(nTrack - ii - 1); + + // Loop through the TRD hits + Bool_t found_hit = kFALSE; + + + for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); + hit; + hit = (AliTRDhit *) fTRD->NextHit()) { + nw = hit->Track(); + if(nw != ii) continue; + det = hit->GetDetector(); + nw = fGeo->GetPlane(det); + if(nw != plane) continue; + + + pos[0]=hit->X(); + pos[1]=hit->Y(); + pos[2]=hit->Z(); + fGeo->Rotate(det,pos,rot); + + if(TMath::Abs(rot[0]-x) > 0.01) continue; + hitY[gtb] = rot[1]; + hitZ[gtb] = rot[2]; + + Float_t col0 = par->GetCol0(plane); + Float_t colPadSize = par->GetColPadSize(plane); + Float_t colH = (Int_t ((rot[1] - col0)/colPadSize)) * colPadSize; + hitO[gtb] = (rot[1] - col0) - colH; + found_hit = kTRUE; + break; + } + + if(!found_hit) continue; + + /* + printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n", + gtb, x, rot[0], rot[1], clusterY[gtb]); + printf("\n Zhit - Zcl = %f - %f = %f\n", + rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]); + + + + printf("found hit within dx = %f - %f \n",rot[0],x); + printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]); + printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]); + printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]); + */ + + + hy->Fill(hitY[gtb]-clusterY[gtb]); + if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]); + else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]); + + if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0)) + hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]); + hz->Fill(hitZ[gtb]-clusterZ[gtb]); + hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]); + hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]); + hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp)); + herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]); + } + } + } + + gStyle->SetOptStat(1); + gStyle->SetOptFit(1); + + TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840); + c->SetFillColor(10); + c->Divide(2,2); + c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw(); + c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw(); + c->cd(3); hyn->Draw(); + c->cd(4); hy4->Draw(); + + TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940); + c1->SetFillColor(10); + c1->Divide(2,2); + c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus"); + c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus"); + c1->cd(3); hy3->Draw(); + c1->cd(4); hy2->Draw(); + +} + + + + + + diff --git a/TRD/AliTRDsaveMCtracks.C b/TRD/AliTRDsaveMCtracks.C new file mode 100644 index 00000000000..281405a4e61 --- /dev/null +++ b/TRD/AliTRDsaveMCtracks.C @@ -0,0 +1,314 @@ +#ifndef __CINT__ + #include + #include "AliTRDtracker.h" + + #include "AliTRDcluster.h" + #include "AliTRDhit.h" + #include "AliTRDv1.h" + #include "AliTRDgeometry.h" + #include "AliTRDparameter.h" + #include "alles.h" + #include "AliTRDmcTrack.h" + + #include "AliTRDtrack.h" + + #include "TFile.h" + #include "TParticle.h" + #include "TStopwatch.h" +#endif + +Int_t Find(Int_t prtcl, TObjArray *mctarray) { + +// Returns index of the mct which corresponds to particle + + Int_t b=0, e=mctarray->GetEntriesFast(), m=(b+e)/2; + AliTRDmcTrack *mct = 0; + while ((bUncheckedAt(m); + if (prtcl > mct->GetTrackIndex()) b=m+1; + else e=m; + } + + if(mct->GetTrackIndex() == prtcl) return m; + + if((m+1) < mctarray->GetEntriesFast()) { + mct = (AliTRDmcTrack*) mctarray->UncheckedAt(m+1); + if(mct->GetTrackIndex() == prtcl) return m+1; + } + + return -1; +} + +void AliTRDsaveMCtracks() { + + TObjArray mctracks(2000); + TObjArray *TRDmcTracks = &mctracks; + + Char_t *alifile = "galice.root"; + Int_t nEvent = 0; + Float_t low_pt_cut = 0.05; + Float_t low_p_cut = 0.02; + Int_t min_no_of_clusters = 10; + + // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD"); + + // Import the Trees for the event nEvent in the file + const Int_t nparticles = gAlice->GetEvent(nEvent); + + printf("found %d particles in event %d \n", nparticles, nEvent); + + // Create TRDmcTracks for each particle + Int_t label = -1, charge = 0; + Bool_t primary; + Float_t mass; + Int_t pdg_code; + + printf("\n"); + + for (Int_t ii=0; iiParticle(ii); + if(p->P() < low_p_cut) continue; + + primary = kTRUE; + if (p->GetFirstMother()>=0) primary = kFALSE; + + pdg_code = (Int_t) p->GetPdgCode(); + + if ((pdg_code == 10010020) || + (pdg_code == 10010030) || + (pdg_code == 50000050) || + (pdg_code == 50000051) || + (pdg_code == 10020040)) { + + mass = 0.; + charge = 0; + } + else { + TParticlePDG *pdg = p->GetPDG(); + charge = (Int_t) pdg->Charge(); + mass=pdg->Mass(); + } + if(charge == 0) continue; + + AliTRDmcTrack *t = new AliTRDmcTrack(ii, primary, mass, charge, pdg_code); + TRDmcTracks->AddLast(t); + } + printf("\r\n"); + + // Loop through TRD clusters and assign indexes to MC tracks + + TFile *geofile =TFile::Open("AliTRDclusters.root"); + AliTRDtracker *Tracker = new AliTRDtracker(geofile); + Tracker->SetEventNumber(nEvent); + + AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry"); + AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter"); + + + Char_t *clusterfile = "AliTRDclusters.root"; + TObjArray carray(2000); + TObjArray *ClustersArray = &carray; + Tracker->ReadClusters(ClustersArray,clusterfile); + + Int_t nClusters = carray.GetEntriesFast(); + + Int_t track, index; + AliTRDmcTrack *tpr = NULL; + + for (Int_t i = 0; i < nClusters; i++) { + + printf("\r assigning cluster %d out of %d", i, nClusters); + AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i); + + for(Int_t j=0; j<3; j++) { + track = cl->GetLabel(j); + if(track < 0) continue; + if(track >= nparticles) + printf("Track index %d is larger than total number of particles %d\n" + ,track,nparticles); + else { + index = Find(track, TRDmcTracks); + tpr = 0; + if(index > 0) tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index); + + if(tpr) { + label = tpr->GetTrackIndex(); + if(label != track) printf("Got label %d while expected %d !\n", + label, track); + TParticle *p = gAlice->Particle(track); + + if (p->Pt() > low_pt_cut) tpr->Update(i); + } + } + } + } + + // Loop through the TRD hits and set XYZ and Pin and Pout + + TTree *hitTree = gAlice->TreeH(); + gAlice->ResetHits(); + Int_t preamp, amp, det, plane = 0, nBytes = 0; + Int_t nTrack = (Int_t) hitTree->GetEntries(); + Float_t pos[3], rot[3]; + Float_t prepos[3], prerot[3]; + + Bool_t entrance = kFALSE; + Bool_t exit = kFALSE; + + + Double_t xin[6], xout[6]; + for(Int_t j = 0; j < 6; j++) { + xin[j] = Tracker->GetX(0,j,14); + xout[j] = Tracker->GetX(0,j,0); + } + + AliTRDhit *hit; + + Int_t current_track = -1; + + printf("\n"); + while(nTrack--) { + + gAlice->ResetHits(); + nBytes += hitTree->GetEvent(nTrack); + entrance = kTRUE; + preamp = 1111; + amp = 1111; + for(Int_t j = 0; j < 3; j++) { pos[j] = 1111.; rot[j] = 1111.;} + + for(hit = (AliTRDhit *) fTRD->FirstHit(-1); + hit; + hit = (AliTRDhit *) fTRD->NextHit()) { + + preamp = amp; + for(Int_t j = 0; j < 3; j++) { prepos[j] = pos[j]; prerot[j] = rot[j];} + + amp = hit->GetCharge(); + if(amp != 0) continue; + track = hit->GetTrack(); + + TParticle *p = gAlice->Particle(track); + if (p->Pt() <= low_pt_cut) continue; + + if(track != current_track) { + current_track = track; + // printf("\n\n"); + } + + det = hit->GetDetector(); + plane = fGeo->GetPlane(det); + pos[0] = hit->X(); + pos[1] = hit->Y(); + pos[2] = hit->Z(); + fGeo->Rotate(det,pos,rot); + + if(TMath::Abs(rot[0] - xin[plane]) < 1.2) entrance = kTRUE; + else entrance = kFALSE; + if(TMath::Abs(rot[0] - xout[plane]) < 1.2) exit = kTRUE; + else exit = kFALSE; + + if(entrance && (preamp == 0) && (prerot[0] < 200)) { + index = Find(track, TRDmcTracks); + if(index > 0) { + tpr = 0; + tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index); + if(tpr) { + label = tpr->GetTrackIndex(); + if(label != track) printf("Got label %d while expected %d !\n", + label, track); + else { + /* + printf("\n momentum at plane %d entrance, Pxyz: %f, %f, %f\n", + plane, prepos[0], prepos[1], prepos[2]); + printf(" XYZ at plane %d entrance: %f, %f, %f\n", + plane, rot[0], rot[1], rot[2]); + printf("expected x range: %f .. %f\n",xin[plane]-1.2,xin[plane]+1.2); + */ + tpr->SetPin(plane, prepos[0], prepos[1], prepos[2]); + tpr->SetXYZin(plane, rot[0], rot[1], rot[2]); + } + } + } + } + + + if(exit && (preamp == 0) && (prerot[0] < 200)) { + index = Find(track, TRDmcTracks); + if(index > 0) { + tpr = 0; + tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index); + if(tpr) { + label = tpr->GetTrackIndex(); + if(label != track) printf("Got label %d while expected %d !\n", + label, track); + else { + /* + printf("\n momentum at plane %d exit, Pxyz: %f, %f, %f\n", + plane, prepos[0], prepos[1], prepos[2]); + printf(" XYZ at plane %d exit: %f, %f, %f\n", + plane, rot[0], rot[1], rot[2]); + printf("expected x range: %f .. %f\n",xout[plane]-1.2,xout[plane]+1.2); + */ + tpr->SetPout(plane, prepos[0], prepos[1], prepos[2]); + tpr->SetXYZout(plane, rot[0], rot[1], rot[2]); + } + } + } + } + } + } + + + // Write TRDmcTracks in output file + + TDirectory *savedir=gDirectory; + Char_t *filename = "AliTRDmcTracks.root"; + TFile *out = new TFile(filename,"RECREATE"); + + TTree *tracktree = new TTree("MCtracks","TRD MC tracks"); + + AliTRDmcTrack *iotrack=0; + tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0); + + Int_t ntracks = TRDmcTracks->GetEntriesFast(); + + for (Int_t i=0; iUncheckedAt(i); + + Int_t n = pt->GetNumberOfClusters(); + if(n > min_no_of_clusters) { + iotrack=pt; + tracktree->Fill(); + printf("Put track with label %d and %d clusters in the output tree \n", + pt->GetTrackIndex(),n); + } + } + + tracktree->Write(); + out->Close(); + savedir->cd(); + + return; +} + diff --git a/TRD/AliTRDtoTOFanalysis.C b/TRD/AliTRDtoTOFanalysis.C new file mode 100644 index 00000000000..581181aa373 --- /dev/null +++ b/TRD/AliTRDtoTOFanalysis.C @@ -0,0 +1,185 @@ +#ifndef __CINT__ + #include + + #include "AliTOF.h" + #include "AliTOFhit.h" + #include "AliTPCtrack.h" + #include "AliTRDtrack.h" + #include "AliTRDgeometry.h" + + #include "alles.h" + #include "TFile.h" + #include "TParticle.h" + #include "TStopwatch.h" +#endif + +void AliTRDtoTOFanalysis() +{ + + Int_t nEvent = 0; + const Int_t maxIndex = 80000; // max number of primaries to be analysed + Int_t rtIndex[maxIndex]; + for(Int_t i = 0; i < maxIndex; i++) rtIndex[i] = -1; + + //****** tracking: Load reconstructed tracks + + TFile *tf=TFile::Open("AliTRDtracks.root"); + + if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;} + TObjArray tarray(2000); + + char tname[100]; + sprintf(tname,"seedsTRDtoTOF2_%d",nEvent); + // for other radial positions use different trees: + // sprintf(tname,"seedsTRDtoTOF1_%d",nEvent); + // sprintf(tname,"seedsTRDtoPHOS_%d",nEvent); + // sprintf(tname,"seedsTRDtoRHIC_%d",nEvent); + + TTree *tracktree=(TTree*)tf->Get(tname); + + TBranch *tbranch=tracktree->GetBranch("tracks"); + + Int_t nRecTracks = (Int_t) tracktree->GetEntries(); + cout<<"Found "<SetAddress(&iotrack); + tracktree->GetEvent(i); + tarray.AddLast(iotrack); + Int_t trackLabel = iotrack->GetLabel(); + + // printf("rt with %d clusters and label %d \n", + // iotrack->GetNumberOfClusters(), trackLabel); + + if(trackLabel < 0) continue; + if(trackLabel >= maxIndex) continue; + rtIndex[trackLabel] = i; + } + tf->Close(); + + //******** + + TH1F *hdx = new TH1F("dx","X(hit) - X(track)",40,-10,10); + TH1F *hdy = new TH1F("dy","Y(hit) - Y(track)",100,-20,20); + TH1F *hdz = new TH1F("dz","Z(hit) - Z(track)",100,-20,20); + TH2F *hdt = new TH2F("dt","T(hit) vs S(track)",100,0,500,500,200,700); + + + // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits + Char_t *alifile = "galice.root"; + + TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); + if (!gafl) { + cout << "Open the ALIROOT-file " << alifile << endl; + gafl = new TFile(alifile); + } + else { + cout << alifile << " is already open" << endl; + } + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*) gafl->Get("gAlice"); + if (gAlice) + cout << "AliRun object found on file" << endl; + else + gAlice = new AliRun("gAlice","Alice test program"); + + + // Import the Trees for the event nEvent in the file + const Int_t nparticles = gAlice->GetEvent(nEvent); + if (nparticles <= 0) return; + + Float_t x,y,z,mass,e; + Double_t gX, gY, gZ; // rec. track coordinates in the global system + Double_t Px, Py, Pz; // rec. track momentum components in the global system + Float_t X, Y, Z; // local tracking coordinates for reconstructed track + Int_t nbytes = 0; + Int_t j,hit,ipart; + Int_t nhits; + Float_t tof; + TParticle *particle; + AliTPCtrack *rt; + +// Get pointers to Alice detectors and Hits containers + AliDetector *TOF = gAlice->GetDetector("TOF"); + Int_t ntracks = (Int_t) gAlice->TreeH()->GetEntries(); + +// Start loop on tracks in the hits containers + + for (Int_t track=0; track < ntracks; track++) { + + if(TOF) { + for(AliTOFhit* tofHit = (AliTOFhit*)TOF->FirstHit(track); + tofHit; + tofHit=(AliTOFhit*)TOF->NextHit()) { + + ipart = tofHit->GetTrack(); + if(ipart >= maxIndex) continue; + if(rtIndex[ipart] < 0) continue; + + x = tofHit->X(); + y = tofHit->Y(); + z = tofHit->Z(); + + //******* tracking: extract track coordinates, momentum, etc. + + rt = (AliTPCtrack*) tarray.UncheckedAt(rtIndex[ipart]); + Double_t tr_length = rt->GetLength(); // meaningfull only with ITS + // part fully implemented + + AliTRDtrack *trd_rt = new AliTRDtrack(*rt, rt->GetAlpha()); + trd_rt->GetGlobalXYZ(gX,gY,gZ); // returns running global coordinates + trd_rt->GetPxPyPz(Px,Py,Pz); // returns momentum in global system + + /* + printf("\n hit position - rec. track position:\n"); + printf(" X: %f - %f = %f \n", x, gX, x - gX); + printf(" Y: %f - %f = %f \n", y, gY, y - gY); + printf(" Z: %f - %f = %f \n", z, gZ, z - gZ); + printf("\n"); + */ + + delete trd_rt; + + // conversion from hit coordinates from global coordinate system + // to "local" tracking system + + Float_t phi =TMath::ATan2(y,x); + if (phi > 2.*TMath::Pi()) phi -= 2.*TMath::Pi(); + if (phi < 0. ) phi += 2.*TMath::Pi(); + Int_t sec = Int_t(phi/AliTRDgeometry::GetAlpha()) % 18; + Float_t alpha = (sec + 0.5) * AliTRDgeometry::GetAlpha(); + + Float_t tmp=x*TMath::Cos(alpha) + y*TMath::Sin(alpha); + y= - x*TMath::Sin(alpha) + y*TMath::Cos(alpha); + x=tmp; + + // particle = (TParticle*)gAlice->Particle(ipart); + // if (particle->GetFirstMother() < 0) continue; + + X = (Float_t) rt->GetX(); // radial position in the tracking system + Y = (Float_t) rt->GetY(); // r*phi position within the sector + Z = (Float_t) rt->GetZ(); // Z position + + + if(TMath::Abs(X-x-1) > 2) continue; + + hdx->Fill(X-x); + hdy->Fill(Y-y); + hdz->Fill(Z-z); + + //********** + + } + } + } + + TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940); + c1->SetFillColor(10); + c1->Divide(2,2); + c1->cd(1); hdx->Draw(); + c1->cd(2); hdy->Draw(); + c1->cd(3); hdz->Draw(); + c1->cd(4); hdt->Draw(); +} diff --git a/TRD/AliTRDtrackReconstruction.C b/TRD/AliTRDtrackReconstruction.C new file mode 100644 index 00000000000..fc2d0016746 --- /dev/null +++ b/TRD/AliTRDtrackReconstruction.C @@ -0,0 +1,112 @@ +#ifndef __CINT__ + #include "alles.h" + #include "AliMagF.h" + #include "AliTRDtracker.h" +#endif + +Int_t TRDPropagateBack(const Char_t *geoname, const Char_t *clrsname, + const Char_t *inname, const Char_t *outname, Int_t n); + +Int_t TRDFindTracks(const Char_t *geoname, const Char_t *clrsname, + const Char_t *inname, const Char_t *outname, Int_t n); + +Int_t AliTRDtrackReconstruction(Int_t n=1) { + const Char_t *TRDdigName="galice.root"; + const Char_t *dummyName="dummy.root"; + const Char_t *TRDclsName="AliTRDclusters.root"; + const Char_t *TRDtrkName="AliTRDtracks.root"; + const Char_t *TPCbkTrkName="AliTPCBackTracks.root"; + + AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor()); + + +// ********** Find TRD tracks from TPC back propagated tracks *********** // + + + if (TRDPropagateBack(TRDclsName, TRDclsName, TPCbkTrkName, TRDtrkName, n)) { + cerr<<"Failed to propagate back through TRD !\n"; + return 1; + } + + +// ********** Find TRD tracks and make seeds for TPC *********** // + + /* + if (TRDFindTracks(TRDclsName,TRDclsName, TRDtrkName, TRDtrkName,n)) { + cerr<<"Failed to find TRD tracks !\n"; + return 1; + } + */ + + return 0; +} + + +Int_t TRDPropagateBack(const Char_t *geoname, const Char_t *clrsname, + const Char_t *inname, const Char_t *outname, Int_t n) +{ + Int_t rc=0; + const Char_t *name="TRDPropagateBack"; + cerr<<'\n'<Start(name); + + TFile *geofile =TFile::Open(geoname); + TFile *out=TFile::Open(outname,"update"); + TFile *in =TFile::Open(inname); + TFile *clrsfile =TFile::Open(clrsname); + + AliTRDtracker *tracker=new AliTRDtracker(geofile); + + for (Int_t i=0;iSetEventNumber(i); + rc=tracker->PropagateBack(in,out); + } + + delete tracker; + in->Close(); + out->Close(); + geofile->Close(); + clrsfile->Close(); + + gBenchmark->Stop(name); + gBenchmark->Show(name); + + return rc; +} + + +Int_t TRDFindTracks(const Char_t *geoname, const Char_t *clrsname, + const Char_t *inname, const Char_t *outname, Int_t n) +{ + Int_t rc=0; + const Char_t *name="TRDFindTracks"; + cerr<<'\n'<Start(name); + + TFile *geofile =TFile::Open(geoname); + TFile *out=TFile::Open(outname,"update"); + TFile *in =TFile::Open(inname); + TFile *clrsfile =TFile::Open(clrsname); + + AliTRDtracker *tracker=new AliTRDtracker(geofile); + tracker->SetAddTRDseeds(); + + for (Int_t i=0;iSetEventNumber(i); + rc=tracker->Clusters2Tracks(in,out); + } + + delete tracker; + in->Close(); + out->Close(); + geofile->Close(); + clrsfile->Close(); + + gBenchmark->Stop(name); + gBenchmark->Show(name); + + return rc; +} + -- 2.43.0