1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////
20 // Fills a set of QA histograms to check the correctness of //
21 // the TRD reconstruction //
23 ////////////////////////////////////////////////////////////////////
25 #include "AliTRDtrackingAnalysis.h"
31 #include "TObjArray.h"
33 #include "TGeoMatrix.h"
35 #include "TGraphErrors.h"
39 #include "AliRunLoader.h"
40 #include "AliTRDgeometry.h"
42 #include "AliESDEvent.h"
43 #include "AliESDtrack.h"
44 #include "AliTrackReference.h"
45 #include "AliTracker.h"
47 #include "AliTRDcluster.h"
48 #include "AliTRDpadPlane.h"
49 #include "AliTRDcalibDB.h"
50 #include "AliTRDtracker.h"
51 //#include "AliTRDtracklet.h"
53 #include "TGeoManager.h"
55 //////////////////////////////////////////////////////////////////////////////////////////
57 AliTRDtrackingAnalysis::AliTRDtrackingAnalysis():
94 fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
95 fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
97 fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
98 fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
100 fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
101 fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
103 fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
104 fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
106 fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
107 fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
111 fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
112 fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
114 for(int i=0; i<12; i++) // bewere hidden constants in the code
115 fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
117 fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
118 fGrResTgPhi = new TGraphErrors();
119 fGrMeanTgPhi = new TGraphErrors();
121 //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
122 //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
125 fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
126 fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
128 fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
129 fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
131 fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
132 fTestBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
134 fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
135 fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
137 fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
138 fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
140 //fGeo = new AliTRDgeometry();
143 //////////////////////////////////////////////////////////////////////////////////////////
145 AliTRDtrackingAnalysis::AliTRDtrackingAnalysis(const AliTRDtrackingAnalysis &t):
182 fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
183 fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
185 fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
186 fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
188 fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
189 fNGood = new TH1D("nGood", ";np", 40, -0.5, 39.5);
191 fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
192 fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
194 fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
195 fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
199 fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
200 fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
202 for(int i=0; i<12; i++) // bewere hidden constants in the code
203 fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
205 fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
206 fGrResTgPhi = new TGraphErrors();
207 fGrMeanTgPhi = new TGraphErrors();
209 //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
210 //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
213 fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
214 fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
216 fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
217 fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
219 fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
220 fTestBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
222 fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
223 fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
225 fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
226 fClZXcl = new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
228 //fGeo = new AliTRDgeometry();
231 //////////////////////////////////////////////////////////////////////////////////////////
233 void AliTRDtrackingAnalysis::DrawResolutionPt(int startEvent, int stopEvent)
236 // Check the pt resolution
241 // loop over ESD events
242 int nevents = fEsdTree->GetEntries();
244 for(int iEvent=startEvent; iEvent<nevents && iEvent < stopEvent; iEvent++) {
246 Info("Draw", "Event = %d", iEvent);
248 fEsdTree->GetEvent(iEvent);
249 fLoader->GetEvent(iEvent);
252 int nTracks = fESD->GetNumberOfTracks();
253 for(int iTrack=0; iTrack<nTracks; iTrack++) {
255 //Info("Track", "Track = %d", iTrack);
256 AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
257 if (!esdTrack->GetInnerParam()) continue;
258 const AliExternalTrackParam *param = esdTrack->GetOuterParam();
259 int status = esdTrack->GetStatus();
261 if (!(status & AliESDtrack::kTRDout)) continue;
262 if (!(status & AliESDtrack::kTRDrefit)) continue;
263 if (esdTrack->GetOuterParam()->Pt() < 1.0) continue;
266 while(param->GetX() > fGeo->GetTime0(ch)+2) ch++;
267 fRefSpace->Fill(2.*ch+0.5, param->GetX() - fGeo->GetTime0(ch));
268 //if (ch < 5) continue;
271 int label = abs(esdTrack->GetTRDLabel());
275 for(int iPoint=GetReference(label); iPoint<fRefTRD->GetEntries(); iPoint++) {
277 AliTrackReference *aRef = (AliTrackReference*)(*fRefTRD)[iPoint];
278 if (aRef->GetTrack() != label) break;
281 lastX = (aRef->LocalX() < lastX)? lastX : aRef->LocalX();
282 double dx = aRef->LocalX() - param->GetX();
283 if (TMath::Abs(dx) > 1.) continue;
286 double bz=fESD->GetMagneticField();
287 AliExternalTrackParam out(*param);
288 out.PropagateTo(aRef->LocalX(),bz);
290 double dp = aRef->Pt() + out.GetSignedPt();
291 double dy = 10. * (aRef->LocalY() - out.GetY()); // in mm
293 fDeltaPt->Fill(100. * dp / aRef->Pt());
294 //fDeltaPt->Fill(out.GetPt() / aRef->Pt());
297 if (esdTrack->GetSign() > 0) fDeltaYPos->Fill(dy);
298 else fDeltaYNeg->Fill(dy);
300 fDeltaZ->Fill(aRef->Z() - out.GetZ());
303 //if (ngood == 0) Info("X", "N = %d, X = %f, DX = %f", ntr, param->GetX(), param->GetX()-lastX);
338 //////////////////////////////////////////////////////////////////////////////////////////
340 // void AliTRDtrackingAnalysis::DrawTrackletResolution(int startEvent, int stopEvent) {
342 // LoadRecPointsFile();
344 // TFile *file = new TFile(Form("%s/TRD.Tracklets.root", fPath), "read");
345 // TTree *tree = (TTree*)file->Get("TRDtracklets");
347 // AliTRDtracklet *tracklet = new AliTRDtracklet();
348 // tree->SetBranchAddress("tracklets", &tracklet);
350 // for(int ev=startEvent; ev<stopEvent; ev++) {
352 // gAlice->GetEvent(ev);
355 // int N = tree->GetEntries();
356 // for(int i=0; i<N; i++) {
358 // tree->GetEntry(i);
360 // Double_t yref, zref, tgphi;
361 // int stat = GetMCPosition(tracklet->GetLabel(), tracklet->GetX(), yref, zref, tgphi);
362 // if (stat < 0) continue;
364 // int plane = tracklet->GetPlane();
365 // Double_t h01 = tracklet->GetTilt();
367 // //printf("Tile = %f\tcorrection = %f um \n", h01, 1e4 * h01 * (tracklet->GetZ()-zref));
368 // //double dz = zref - tracklet->GetZ() - cls->GetZ();
370 // fTrklY->Fill(10 * (tracklet->GetY() - yref));
371 // fTrklZ->Fill(tracklet->GetZ() - zref);
374 // while(tracklet->GetX() > fGeo->GetTime0(ch)+2) ch++;
375 // fRefSpace->Fill(tracklet->GetY(), tracklet->GetX() - fGeo->GetTime0(ch));
383 // fRefSpace->Draw();
385 // gStyle->SetOptFit(1);
388 // fTrklY->Fit("gaus");
394 //////////////////////////////////////////////////////////////////////////////////////////
396 void AliTRDtrackingAnalysis::DrawRecPointResolution(int startEvent, int stopEvent)
399 // Check the resolution of the reconstructed points
403 TObjArray *module = new TObjArray();
405 int nEvents = gAlice->GetEventsPerRun();
407 for(int ev=startEvent; ev<nEvents && ev < stopEvent; ev++) {
409 gAlice->GetEvent(ev);
412 TTree *tree = fLoader->GetTreeR("TRD", 0);
413 tree->SetBranchAddress("TRDcluster", &module);
415 Info("Res", "Refs Loaded");
417 Int_t nn = tree->GetEntries();
418 for(int i=0; i<nn; i++) {
421 int m = module->GetEntries();
423 for(int j=0; j<m; j++) {
425 AliTRDcluster *cls = (AliTRDcluster*)module->At(j);
426 if (cls->GetQ() < 10) continue;
427 //fTracker->Transform(cls);
428 fClPos->Fill(cls->GetZ(), cls->GetY());
430 int layer = fGeo->GetLayer(cls->GetDetector());
433 for(int k=0; k<3; k++) if (cls->GetLabel(k) > -1) nl++;
436 Double_t yref, zref, tgphi;
437 int stat = GetMCPosition(cls->GetLabel(0), cls->GetX(), yref, zref, tgphi);
438 if (stat < 0) continue;
440 fClZXcl->Fill(cls->GetZ(), cls->GetX());
441 fClZXref->Fill(zref, cls->GetX());
443 AliTRDpadPlane *padPlane = fGeo->GetPadPlane(layer,0);
444 Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
446 //double dz = zref - padPlane->GetRow0();
447 double dz = zref - cls->GetZ();
448 double dy = dz * h01;
449 double yy = cls->GetY() - dy;
451 if (cls->GetNPads() == 2) fClY2->Fill(10 * (yy - yref));
452 if (cls->GetNPads() == 3) fClY3->Fill(10 * (yy - yref));
454 int idx = GetPhiBin(tgphi);
455 if (idx >= 0 && idx < 12) fClYTgPhi[idx]->Fill(10 * (yy - yref));
457 fClZZ->Fill(zref, cls->GetZ());
460 fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(layer));
476 gStyle->SetOptFit(1);
481 fClY2->Fit("gaus", "", "", -2*fClY2->GetRMS(), 2*fClY2->GetRMS());
486 fClY3->Fit("gaus", "", "", -2*fClY3->GetRMS(), 2*fClY3->GetRMS());
498 TCanvas *c = new TCanvas();
501 for(int i=0; i<12; i++) {
504 fClYTgPhi[i]->Draw();
505 if (fClYTgPhi[i]->GetSum() < 100) continue;
507 double mean = fClYTgPhi[i]->GetMean();
508 double rms = fClYTgPhi[i]->GetRMS();
510 fClYTgPhi[i]->Fit("gaus", "", "", mean-2*rms, mean+2*rms);
511 TF1 *f = fClYTgPhi[i]->GetFunction("gaus");
513 int n = fGrResTgPhi->GetN();
514 fGrResTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(2));
515 fGrResTgPhi->SetPointError(n, 0, f->GetParError(2));
517 fGrMeanTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(1));
518 fGrMeanTgPhi->SetPointError(n, 0, f->GetParError(1));
521 //gSystem->Sleep(1000);
522 gStyle->SetOptStat(0);
525 TH1D *dummy = new TH1D("dummy", "", 100, -0.3, 0.3);
527 dummy->SetTitle(";tg(#phi);resolution (mm)");
528 dummy->SetMinimum(0);
529 dummy->SetMaximum(1);
532 (dummy->Clone("dummy1"))->Draw();
534 fGrResTgPhi->Draw("PL");
535 fGrResTgPhi->GetHistogram()->SetTitle(";tg(#phi);resolution (mm)");
536 fGrResTgPhi->SetMarkerStyle(20);
539 dummy->SetTitle(";tg(#phi);mean value (mm)");
540 dummy->SetMinimum(-0.3);
541 dummy->SetMaximum(0.3);
544 (dummy->Clone("dummy2"))->Draw();
546 fGrMeanTgPhi->Draw("PL");
547 fGrMeanTgPhi->GetHistogram()->SetTitle(";tg(#phi);mean value (mm)");
548 fGrMeanTgPhi->SetMarkerStyle(20);
553 //fClZZ->Draw("colz");
556 //fClPos->Draw("colz");
566 //////////////////////////////////////////////////////////////////////////////////////////
568 void AliTRDtrackingAnalysis::LoadRecPointsFile()
571 // Load the clusters from the input file
575 sprintf(filename, "%s/galice.root", fPath);
577 fLoader = AliRunLoader::Open(filename);
579 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
583 fLoader->LoadgAlice();
584 gAlice = fLoader->GetAliRun();
587 Error("CheckFiles", "no galice object found");
591 fLoader->LoadKinematics();
592 fLoader->LoadHeader();
593 fLoader->LoadTrackRefs();
595 //TGeoManager::Import("/data/alice_u/radomski/condor/run_0/geometry.root");
596 TGeoManager::Import(Form("%s/geometry.root", fPath));
600 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
601 fTracker = new AliTRDtracker(gFile);
603 fLoader->LoadRecPoints("TRD");
605 AliTracker::SetFieldMap(gAlice->Field(), 1);
608 //////////////////////////////////////////////////////////////////////////////////////////
610 void AliTRDtrackingAnalysis::CheckFiles()
613 // Check the presence of the input files
619 sprintf(filename, "%s/galice.root", fPath);
621 fLoader = AliRunLoader::Open(filename);
623 Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
627 fLoader->LoadgAlice();
628 gAlice = fLoader->GetAliRun();
631 Error("CheckFiles", "no galice object found");
635 fLoader->LoadKinematics();
636 fLoader->LoadHeader();
637 fLoader->LoadTrackRefs();
640 fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
641 //fGeo->ReadGeoMatrices();
645 sprintf(filename,"%s/AliESDs.root", fPath);
646 TFile *esdFile = new TFile(filename, "READ");
648 if (esdFile->IsZombie()) {
649 Error("CheckFiles", "file not present: AliESDs.root");
653 fEsdTree = (TTree*)esdFile->Get("esdTree");
654 fESD = new AliESDEvent();
655 fESD->ReadFromTree(fEsdTree);
656 //fEsdTree->SetBranchAddress("ESD", &fESD);
659 //////////////////////////////////////////////////////////////////////////////////////////
661 void AliTRDtrackingAnalysis::LoadRefs()
664 // Load the track references
667 if (fRefTPC) delete fRefTPC;
668 if (fRefTRD) delete fRefTRD;
670 fRefTPC = new TObjArray();
671 fRefTRD = new TObjArray();
673 //fLoader->GetEvent(event);
674 //AliStack* stack = gAlice->Stack();
675 TTree *refTree = fLoader->TreeTR();
677 TClonesArray *clRefs = new TClonesArray("AliTrackReference");
679 TBranch *branch = refTree->GetBranch("TrackReferences");
680 refTree->SetBranchAddress("TrackReferences",&clRefs);
682 int nEntries = branch->GetEntries();
683 for(int iTrack = 0; iTrack < nEntries; iTrack++) {
685 refTree->GetEvent(iTrack);
686 int nPoints = clRefs->GetEntries();
687 for(int iPoint=0; iPoint<nPoints; iPoint++) {
688 AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
689 if (ref->DetectorId() == AliTrackReference::kTPC) fRefTPC->Add(new AliTrackReference(*ref));
690 if (ref->DetectorId() == AliTrackReference::kTRD) fRefTRD->Add(new AliTrackReference(*ref));
697 for(int i=0; i<fRefTRD->GetEntries(); i++) {
698 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
699 fLabels[i] = ref->GetTrack();
702 while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
703 fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
705 //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fTestBits->Fill(bit);
709 Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
712 //////////////////////////////////////////////////////////////////////////////////////////
714 Int_t AliTRDtrackingAnalysis::GetReference(Int_t label)
717 // Sort the track references
720 int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
723 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
724 if (ref->GetTrack() != label) return start+1;
731 //////////////////////////////////////////////////////////////////////////////////////////
733 int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi)
736 // Determine the MC positions from the track references
745 int idx = GetReference(label);
746 for(int i=idx; i<fRefTRD->GetEntries(); i++) {
748 AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
749 if (ref->GetTrack() != label) break;
753 //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
754 //if (p != layer) continue;
756 double dX = ref->LocalX()-x;
771 if (idLow == -1 || idHigh == -1) return -1;
773 AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
774 AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
777 double dx = refO->LocalX() - refI->LocalX();
778 double dy = refO->LocalY() - refI->LocalY();
779 double dz = refO->Z() - refI->Z();
780 double ddx = (x - refI->LocalX())/dx;
784 Y = refI->LocalY() + ddx * dy;
785 Z = refI->Z() + ddx * dz;
792 //////////////////////////////////////////////////////////////////////////////////////////
793 Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi) const
796 // Return the phi bin
798 return (int)((phi+0.3)/0.05);
801 //////////////////////////////////////////////////////////////////////////////////////////
802 Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) const
805 // Return phi for a given bin
807 return bin * 0.05 - 0.3 + 0.025;
809 //////////////////////////////////////////////////////////////////////////////////////////