--- /dev/null
+
+#include "AliTRDtrackingAnalysis.h"
+
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TObjArray.h"
+#include "TCanvas.h"
+#include "TGeoMatrix.h"
+#include "TStyle.h"
+#include "TGraphErrors.h"
+#include "TF1.h"
+
+#include "AliRunLoader.h"
+#include "AliTRDgeometry.h"
+#include "AliRun.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliTrackReference.h"
+
+#include "AliTRDcluster.h"
+#include "AliTRDCommonParam.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDcalibDB.h"
+#include "AliTracker.h"
+#include "AliTRDtracker.h"
+//#include "AliTRDtracklet.h"
+
+#include "TGeoManager.h"
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+AliTRDtrackingAnalysis::AliTRDtrackingAnalysis():
+ TObject(),
+ fLoader(0),
+ fEsdTree(0),
+ fESD(0),
+ fRefTPC(0),
+ fRefTRD(0)
+{
+
+ fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
+ fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
+
+ fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
+ fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
+
+ fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
+ fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
+
+ fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
+ fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
+
+ fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
+ fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
+
+
+ // cluster studies
+ fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
+ fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
+
+ for(int i=0; i<12; i++) // bewere hidden constants in the code
+ fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
+
+ fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
+ fGrResTgPhi = new TGraphErrors();
+ fGrMeanTgPhi = new TGraphErrors();
+
+ //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
+ //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
+
+
+ fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
+ fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
+
+ fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
+ fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
+
+ fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
+ fBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
+
+ fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
+ fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
+
+ fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
+ fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
+
+ //fGeo = new AliTRDgeometry();
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void AliTRDtrackingAnalysis::DrawResolutionPt(int startEvent, int stopEvent) {
+
+ CheckFiles();
+
+ // loop over ESD events
+ int nevents = fEsdTree->GetEntries();
+
+ for(int iEvent=startEvent; iEvent<nevents && iEvent < stopEvent; iEvent++) {
+
+ Info("Draw", "Event = %d", iEvent);
+
+ fEsdTree->GetEvent(iEvent);
+ fLoader->GetEvent(iEvent);
+ LoadRefs();
+
+ int nTracks = fESD->GetNumberOfTracks();
+ for(int iTrack=0; iTrack<nTracks; iTrack++) {
+
+ //Info("Track", "Track = %d", iTrack);
+ AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
+ if (!esdTrack->GetInnerParam()) continue;
+ const AliExternalTrackParam *param = esdTrack->GetOuterParam();
+ int status = esdTrack->GetStatus();
+
+ if (!(status & AliESDtrack::kTRDout)) continue;
+ if (!(status & AliESDtrack::kTRDrefit)) continue;
+ if (fabs(esdTrack->GetOuterParam()->GetPt()) < 1.0) continue;
+
+ int ch=0;
+ while(param->GetX() > fGeo->GetTime0(ch)+2) ch++;
+ fRefSpace->Fill(2.*ch+0.5, param->GetX() - fGeo->GetTime0(ch));
+ //if (ch < 5) continue;
+
+ double lastX = 0;
+ int label = abs(esdTrack->GetTRDLabel());
+ int ntr = 0;
+ int ngood = 0;
+
+ for(int iPoint=GetReference(label); iPoint<fRefTRD->GetEntries(); iPoint++) {
+
+ AliTrackReference *aRef = (AliTrackReference*)(*fRefTRD)[iPoint];
+ if (aRef->GetTrack() != label) break;
+ ntr++;
+
+ lastX = (aRef->LocalX() < lastX)? lastX : aRef->LocalX();
+ double dx = aRef->LocalX() - param->GetX();
+ if (fabs(dx) > 1.) continue;
+ ngood++;
+
+ double bz=fESD->GetMagneticField();
+ AliExternalTrackParam out(*param);
+ out.PropagateTo(aRef->LocalX(),bz);
+
+ double dp = aRef->Pt() + out.GetPt();
+ double dy = 10. * (aRef->LocalY() - out.GetY()); // in mm
+
+ fDeltaPt->Fill(100. * dp / aRef->Pt());
+ //fDeltaPt->Fill(out.GetPt() / aRef->Pt());
+ fDeltaX->Fill(dx);
+
+ if (esdTrack->GetSign() > 0) fDeltaYPos->Fill(dy);
+ else fDeltaYNeg->Fill(dy);
+
+ fDeltaZ->Fill(aRef->Z() - out.GetZ());
+ }
+
+ //if (ngood == 0) Info("X", "N = %d, X = %f, DX = %f", ntr, param->GetX(), param->GetX()-lastX);
+
+ fNPoints->Fill(ntr);
+ fNGood->Fill(ngood);
+ }
+ }
+
+ new TCanvas();
+ fDeltaPt->Draw();
+
+ new TCanvas();
+ fDeltaX->Draw();
+
+ /*
+ new TCanvas();
+ fNPoints->Draw();
+
+ new TCanvas();
+ fNGood->Draw();
+ */
+
+ new TCanvas();
+ fDeltaYPos->Draw();
+
+ new TCanvas();
+ fDeltaYNeg->Draw();
+
+ new TCanvas();
+ fDeltaZ->Draw();
+
+ new TCanvas();
+ fRefSpace->Draw();
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+// void AliTRDtrackingAnalysis::DrawTrackletResolution(int startEvent, int stopEvent) {
+
+// LoadRecPointsFile();
+
+// TFile *file = new TFile(Form("%s/TRD.Tracklets.root", fPath), "read");
+// TTree *tree = (TTree*)file->Get("TRDtracklets");
+
+// AliTRDtracklet *tracklet = new AliTRDtracklet();
+// tree->SetBranchAddress("tracklets", &tracklet);
+
+// for(int ev=startEvent; ev<stopEvent; ev++) {
+
+// gAlice->GetEvent(ev);
+// LoadRefs();
+
+// int N = tree->GetEntries();
+// for(int i=0; i<N; i++) {
+
+// tree->GetEntry(i);
+
+// Double_t yref, zref, tgphi;
+// int stat = GetMCPosition(tracklet->GetLabel(), tracklet->GetX(), yref, zref, tgphi);
+// if (stat < 0) continue;
+
+// int plane = tracklet->GetPlane();
+// Double_t h01 = tracklet->GetTilt();
+
+// //printf("Tile = %f\tcorrection = %f um \n", h01, 1e4 * h01 * (tracklet->GetZ()-zref));
+// //double dz = zref - tracklet->GetZ() - cls->GetZ();
+
+// fTrklY->Fill(10 * (tracklet->GetY() - yref));
+// fTrklZ->Fill(tracklet->GetZ() - zref);
+
+// int ch=0;
+// while(tracklet->GetX() > fGeo->GetTime0(ch)+2) ch++;
+// fRefSpace->Fill(tracklet->GetY(), tracklet->GetX() - fGeo->GetTime0(ch));
+// }
+// }
+
+// new TCanvas();
+// fTrklZ->Draw();
+
+// new TCanvas();
+// fRefSpace->Draw();
+
+// gStyle->SetOptFit(1);
+// new TCanvas();
+// fTrklY->Draw();
+// fTrklY->Fit("gaus");
+
+// //new TCanvas();
+// //fClZ->Draw();
+// }
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void AliTRDtrackingAnalysis::DrawRecPointResolution(int startEvent, int stopEvent) {
+
+ LoadRecPointsFile();
+ TObjArray *module = new TObjArray();
+
+ int nEvents = gAlice->GetEventsPerRun();
+
+ for(int ev=startEvent; ev<nEvents && ev < stopEvent; ev++) {
+
+ gAlice->GetEvent(ev);
+ LoadRefs();
+
+ TTree *tree = fLoader->GetTreeR("TRD", 0);
+ tree->SetBranchAddress("TRDcluster", &module);
+
+ Info("Res", "Refs Loaded");
+
+ int N = tree->GetEntries();
+ for(int i=0; i<N; i++) {
+
+ tree->GetEntry(i);
+ int m = module->GetEntries();
+
+ for(int j=0; j<m; j++) {
+
+ AliTRDcluster *cls = (AliTRDcluster*)module->At(j);
+ if (cls->GetQ() < 10) continue;
+ fTracker->Transform(cls);
+ fClPos->Fill(cls->GetZ(), cls->GetY());
+
+ int plane = fGeo->GetPlane(cls->GetDetector());
+
+ int nl = 0;
+ for(int k=0; k<3; k++) if (cls->GetLabel(k) > -1) nl++;
+ fNLabels->Fill(nl);
+
+ Double_t yref, zref, tgphi;
+ int stat = GetMCPosition(cls->GetLabel(0), cls->GetX(), yref, zref, tgphi);
+ if (stat < 0) continue;
+
+ fClZXcl->Fill(cls->GetZ(), cls->GetX());
+ fClZXref->Fill(zref, cls->GetX());
+
+ AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
+ Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
+
+ //double dz = zref - padPlane->GetRow0();
+ double dz = zref - cls->GetZ();
+ double dy = dz * h01;
+ double yy = cls->GetY() - dy;
+
+ if (cls->From2pad()) fClY2->Fill(10 * (yy - yref));
+ if (cls->From3pad()) fClY3->Fill(10 * (yy - yref));
+
+ int idx = GetPhiBin(tgphi);
+ if (idx >= 0 && idx < 12) fClYTgPhi[idx]->Fill(10 * (yy - yref));
+
+ fClZZ->Fill(zref, cls->GetZ());
+ fClZ->Fill(dz);
+ fTgPhi->Fill(tgphi);
+ fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(plane));
+ }
+ }
+ }
+
+ new TCanvas();
+ fClYX->Draw();
+
+ //new TCanvas();
+ //fNLabels->Draw();
+ new TCanvas();
+ fClZ->Draw();
+
+ new TCanvas();
+ fTgPhi->Draw();
+
+ gStyle->SetOptFit(1);
+
+ new TCanvas();
+ gPad->SetLogy();
+ fClY2->Draw();
+ fClY2->Fit("gaus", "", "", -2*fClY2->GetRMS(), 2*fClY2->GetRMS());
+
+ new TCanvas();
+ gPad->SetLogy();
+ fClY3->Draw();
+ fClY3->Fit("gaus", "", "", -2*fClY3->GetRMS(), 2*fClY3->GetRMS());
+
+ //new TCanvas();
+ //fRefSpace->Draw();
+
+ //new TCanvas();
+ //fClZXcl->Draw();
+
+ //new TCanvas();
+ //fClZXref->Draw();
+
+ /**/
+ TCanvas *c = new TCanvas();
+ c->Divide(4,3);
+
+ for(int i=0; i<12; i++) {
+
+ c->cd(i+1);
+ fClYTgPhi[i]->Draw();
+ if (fClYTgPhi[i]->GetSum() < 100) continue;
+
+ double mean = fClYTgPhi[i]->GetMean();
+ double rms = fClYTgPhi[i]->GetRMS();
+
+ fClYTgPhi[i]->Fit("gaus", "", "", mean-2*rms, mean+2*rms);
+ TF1 *f = fClYTgPhi[i]->GetFunction("gaus");
+
+ int n = fGrResTgPhi->GetN();
+ fGrResTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(2));
+ fGrResTgPhi->SetPointError(n, 0, f->GetParError(2));
+
+ fGrMeanTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(1));
+ fGrMeanTgPhi->SetPointError(n, 0, f->GetParError(1));
+ }
+
+ //gSystem->Sleep(1000);
+ gStyle->SetOptStat(0);
+
+ c = new TCanvas();
+ TH1D *dummy = new TH1D("dummy", "", 100, -0.3, 0.3);
+
+ dummy->SetTitle(";tg(#phi);resolution (mm)");
+ dummy->SetMinimum(0);
+ dummy->SetMaximum(1);
+
+ //c->cd();
+ (dummy->Clone("dummy1"))->Draw();
+
+ fGrResTgPhi->Draw("PL");
+ fGrResTgPhi->GetHistogram()->SetTitle(";tg(#phi);resolution (mm)");
+ fGrResTgPhi->SetMarkerStyle(20);
+
+ c = new TCanvas();
+ dummy->SetTitle(";tg(#phi);mean value (mm)");
+ dummy->SetMinimum(-0.3);
+ dummy->SetMaximum(0.3);
+
+ //c->cd();
+ (dummy->Clone("dummy2"))->Draw();
+
+ fGrMeanTgPhi->Draw("PL");
+ fGrMeanTgPhi->GetHistogram()->SetTitle(";tg(#phi);mean value (mm)");
+ fGrMeanTgPhi->SetMarkerStyle(20);
+ /**/
+
+
+ //new TCanvas();
+ //fClZZ->Draw("colz");
+
+ //new TCanvas();
+ //fClPos->Draw("colz");
+
+ //new TCanvas();
+ //fBits->Draw();
+
+ //new TCanvas();
+ //fRefDx->Draw();
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void AliTRDtrackingAnalysis::LoadRecPointsFile() {
+
+ char filename[256];
+ sprintf(filename, "%s/galice.root", fPath);
+
+ fLoader = AliRunLoader::Open(filename);
+ if (!fLoader) {
+ Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
+ return;
+ }
+
+ fLoader->LoadgAlice();
+ gAlice = fLoader->GetAliRun();
+
+ if (!gAlice) {
+ Error("CheckFiles", "no galice object found");
+ return;
+ }
+
+ fLoader->LoadKinematics();
+ fLoader->LoadHeader();
+ fLoader->LoadTrackRefs();
+
+ //TGeoManager::Import("/data/alice_u/radomski/condor/run_0/geometry.root");
+ TGeoManager::Import(Form("%s/geometry.root", fPath));
+
+
+ fLoader->CdGAFile();
+ fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
+ fTracker = new AliTRDtracker(gFile);
+
+ fLoader->LoadRecPoints("TRD");
+
+ AliTracker::SetFieldMap(gAlice->Field(), 1);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void AliTRDtrackingAnalysis::CheckFiles() {
+
+ // MC info
+
+ char filename[256];
+ sprintf(filename, "%s/galice.root", fPath);
+
+ fLoader = AliRunLoader::Open(filename);
+ if (!fLoader) {
+ Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
+ return;
+ }
+
+ fLoader->LoadgAlice();
+ gAlice = fLoader->GetAliRun();
+
+ if (!gAlice) {
+ Error("CheckFiles", "no galice object found");
+ return;
+ }
+
+ fLoader->LoadKinematics();
+ fLoader->LoadHeader();
+ fLoader->LoadTrackRefs();
+
+ fLoader->CdGAFile();
+ fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
+ //fGeo->ReadGeoMatrices();
+
+ // ESD
+
+ sprintf(filename,"%s/AliESDs.root", fPath);
+ TFile *esdFile = new TFile(filename, "READ");
+
+ if (esdFile->IsZombie()) {
+ Error("CheckFiles", "file not present: AliESDs.root");
+ return;
+ }
+
+ fEsdTree = (TTree*)esdFile->Get("esdTree");
+ fESD = new AliESD();
+
+ fEsdTree->SetBranchAddress("ESD", &fESD);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+void AliTRDtrackingAnalysis::LoadRefs() {
+
+ if (fRefTPC) delete fRefTPC;
+ if (fRefTRD) delete fRefTRD;
+
+ fRefTPC = new TObjArray();
+ fRefTRD = new TObjArray();
+
+ //fLoader->GetEvent(event);
+ //AliStack* stack = gAlice->Stack();
+ TTree *refTree = fLoader->TreeTR();
+
+ const int nBranch = 2;
+ const char *brName[] = {"TPC", "TRD"};
+ TClonesArray *clRefs = new TClonesArray("AliTrackReference");
+
+ for(int b=0; b<nBranch; b++) {
+
+ TBranch *branch = refTree->GetBranch(brName[b]);
+ refTree->SetBranchAddress(brName[b],&clRefs);
+
+ int nEntries = branch->GetEntries();
+ for(int iTrack = 0; iTrack < nEntries; iTrack++) {
+
+ refTree->GetEvent(iTrack);
+ int nPoints = clRefs->GetEntries();
+ for(int iPoint=0; iPoint<nPoints; iPoint++) {
+ AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
+ if (b == 0) fRefTPC->Add(new AliTrackReference(*ref));
+ if (b == 1) fRefTRD->Add(new AliTrackReference(*ref));
+ }
+ }
+ }
+
+ fRefTPC->Sort();
+ fRefTRD->Sort();
+
+ for(int i=0; i<fRefTRD->GetEntries(); i++) {
+ AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
+ fLabels[i] = ref->GetTrack();
+
+ int p=0;
+ while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
+ fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
+
+ //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fBits->Fill(bit);
+ }
+
+ delete clRefs;
+ Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+Int_t AliTRDtrackingAnalysis::GetReference(Int_t label) {
+
+ int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
+
+ while (start >= 0) {
+ AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
+ if (ref->GetTrack() != label) return start+1;
+ start--;
+ }
+
+ return 0;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi) {
+
+ double lowX = 100.;
+ double highX = 100.;
+ int idLow = -1;
+ int idHigh = -1;
+
+ int nref= 0;
+ int idx = GetReference(label);
+ for(int i=idx; i<fRefTRD->GetEntries(); i++) {
+
+ AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
+ if (ref->GetTrack() != label) break;
+ nref++;
+
+ //int p=0;
+ //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
+ //if (p != layer) continue;
+
+ double dX = ref->LocalX()-x;
+ if ( dX > 0 ) {
+ if (dX < highX) {
+ idHigh = i;
+ highX = dX;
+ }
+ } else {
+ dX = fabs(dX);
+ if (dX < lowX) {
+ idLow = i;
+ lowX = dX;
+ }
+ }
+ }
+
+ if (idLow == -1 || idHigh == -1) return -1;
+
+ AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
+ AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
+
+
+ double dx = refO->LocalX() - refI->LocalX();
+ double dy = refO->LocalY() - refI->LocalY();
+ double dz = refO->Z() - refI->Z();
+ double ddx = (x - refI->LocalX())/dx;
+
+ fRefDx->Fill(dx);
+
+ Y = refI->LocalY() + ddx * dy;
+ Z = refI->Z() + ddx * dz;
+
+ tgphi = dy/dx;
+
+ return 0;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi) {
+ return (int)((phi+0.3)/0.05);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////
+
+Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) {
+ return bin * 0.05 - 0.3 + 0.025;
+}
+//////////////////////////////////////////////////////////////////////////////////////////
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff