--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+ #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; i<n; i++) {
+ seedTree->GetEvent(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 "<<nRecTracks<<" entries in the track tree"<<endl;
+
+ for (Int_t i=0; i<nRecTracks; i++) {
+ AliTRDtrack *iotrack=new AliTRDtrack();
+ tbranch->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 "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
+ for (Int_t i=0; i<nMCtracks; i++) {
+ AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack;
+ mctbranch->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; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+ }
+ else {
+ min_tb = Tracker->GetGlobalTimeBin(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; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
+ max_tb = 0;
+ mcX = Tracker->GetX(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();
+ */
+
+}
+
+
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+
+ #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<No_of_tracks_to_analyze; ii++) {
+ printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze);
+
+ for(Int_t plane = 0; plane < nPlanes; plane++) {
+ for(Int_t ltb = 14; ltb > -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();
+
+}
+
+
+
+
+
+
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+ #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 <prtcl>
+
+ Int_t b=0, e=mctarray->GetEntriesFast(), m=(b+e)/2;
+ AliTRDmcTrack *mct = 0;
+ while ((b<e)&&(e!=0)) {
+ m=(b+e)/2;
+ mct = (AliTRDmcTrack*) mctarray->UncheckedAt(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; ii<nparticles; ii++) {
+
+ printf("\r particle %d out of %d",ii,nparticles);
+
+ TParticle *p = gAlice->Particle(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; i<ntracks; i++) {
+ AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(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;
+}
+
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+
+ #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 "<<nRecTracks<<" entries in the track tree "<<tname<<endl;
+
+ for (Int_t i=0; i<nRecTracks; i++) {
+ AliTPCtrack *iotrack=new AliTPCtrack();
+ tbranch->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();
+}
--- /dev/null
+#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'<<name<<"...\n";
+ gBenchmark->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;i<n;i++){
+ printf("Processing event %d\n",i);
+ tracker->SetEventNumber(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'<<name<<"...\n";
+ gBenchmark->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;i<n;i++){
+ printf("Processing event %d\n",i);
+ tracker->SetEventNumber(i);
+ rc=tracker->Clusters2Tracks(in,out);
+ }
+
+ delete tracker;
+ in->Close();
+ out->Close();
+ geofile->Close();
+ clrsfile->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return rc;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff