Track time measurement (S.Radomski)
[u/mrichter/AliRoot.git] / macros / CompareTime.C
1 //
2 // This macro draws time resolution for true hypothesis and fits it with gausian
3 // on a given reference plane in a given momentum range.
4 // The results is then printed to a PostScript file.
5 // The filename is composed out of input parameters.
6 //
7 // Input parameters:
8 // nRefPlane - indicate the reference plane.
9 //   1: in TPC  
10 //   2: out TPC
11 //   3: out TRD
12 // 
13 // minP, maxP - total momentum cuts. momentum from gAlice are used for cuts
14 // 
15 // The macro is using SORTED reference points in a file "trackRef.root".
16 // Track references can be sorted using CopyReferences.C
17 // 
18 // Before using check filenames and names of trees invlolved
19 // 
20 // ITS tracks (nRefPlane == 1) are compared at inner TPC ref Plane
21 // make shure they are properly propagated to this plane
22 //
23 // Sylwester Radomski, e-mail: S.Radomski@gsi.de
24 // Feb 14, 2002 
25 // 
26
27
28 CompareTime (Int_t nRefPlane, Double_t minP, Double_t maxP) 
29 {
30   
31   // Check input data Consistancy  
32   if (nRefPlane < 1 || nRefPlane > 3) {
33     cout << "Wrong Reference Plane Id ( " << nRefPalne << ") [1-3] " << endl;
34     return;
35   }
36
37   if (maxP < minP) {
38     cout << "MaxP lesser than MinP" << endl;
39     return;
40   }
41
42   // Style
43
44   gROOT->SetStyle("Plain");
45   gStyle->SetOptStat(1110);
46   gStyle->SetOptFit(11);
47
48   Double_t cc = 0.299792458;
49   AliKalmanTrack::SetConvConst(100/cc/0.4);
50   
51   const char* gFileName = "galice.root";
52   const char* refFileName = "trackRef.root";
53
54   const char* trackFileName[] = {"AliTPCtracks.root", "AliTPCBackTracks.root", "AliTRDtracks.root"};
55   const char *treeName[] = {"TreeT_ITSb_0", "seedsTPCtoTRD_0", "TRDb_0"};
56   const char *refTreeName[] = {"TPC", "TPC", "TRD"};
57
58   // Histograms
59
60   const Int_t nTypes = 5; 
61   Int_t codes[] = {11, 13, 211, 321, 2212};
62   
63   //Double_t cutoff[] = {20, 20, 20, 100, 100};
64   Double_t cutoff[] = {30, 30, 30, 100, 200};
65
66   const char *names[nTypes] = {"Electrons", "Muons", "Pions", "Kaons", "Protons"};
67   
68   TH1D *hDiffSame[nTypes];;
69
70   for(Int_t i=0; i<nTypes; i++)
71     hDiffSame[i] = new TH1D(names[i], names[i], 50, -cutoff[i], cutoff[i]);
72   
73   TH1D *hLength = new TH1D("length", "Length Resolution", 100, -0.5, 0.5);
74
75   // Particles
76
77   TFile *refFile = new TFile(gFileName, "READ");
78   gAlice = (AliRun*)refFile->Get("gAlice");
79   gAlice->GetEvent(0);
80
81
82   // Reference tracks
83
84   refFile = new TFile(refFileName, "READ");
85   TTree *treeRef = (TTree*)refFile->Get("TreeTR0_Sorted"); 
86   TBranch *trkRef = treeRef->GetBranch( refTreeName[nRefPlane-1] );
87   TClonesArray *trkRefs = new TClonesArray("AliTrackReference", 100);
88   trkRef->SetAddress(&trkRefs);
89   
90   Int_t ntracks = trkRef->GetEntries();
91   cout << ntracks << endl;
92
93   // Tracks
94   
95   TFile *trackFile = new TFile(trackFileName[nRefPlane-1]);
96
97   TTree *tracks = (TTree*)trackFile->Get(treeName[nRefPlane-1]);
98   TBranch *recTracks = tracks->GetBranch("tracks");
99   AliTPCtrack *track = 0;
100   recTracks->SetAddress(&track);
101
102   for(Int_t t=0; t<tracks->GetEntries(); t++) {
103
104     cout << t << "\r";
105     
106     recTracks->GetEvent(t);
107     
108     Int_t lab = track->GetLabel();
109     if (lab < 0 || lab > 100000) continue;
110     
111     Int_t pdg = gAlice->Particle(lab)->GetPdgCode();
112     treeRef->GetEvent(lab);
113     
114     Int_t nEntrRef = trkRefs->GetEntries();
115     
116     if (gAlice->Particle(lab)->P() < minP ) continue; 
117     if (gAlice->Particle(lab)->P() > maxP ) continue;
118     if (gAlice->Particle(lab)->GetFirstMother() != -1) continue;
119     
120     if (!(nEntrRef)) continue;
121       
122     Double_t timeTrue, timeTrack;
123     Double_t trueLength;
124
125     Int_t index =  nEntrRef-1;
126     if (nRefPlane == 1) index = 0;
127
128     timeTrue = ((AliTrackReference*)(*trkRefs)[index])->GetTime() * 1e12;
129     trueLength = ((AliTrackReference*)(*trkRefs)[index])->GetLength();
130     
131     hLength->Fill(trueLength - track->GetIntegratedLength());
132     
133     for(Int_t i=0; i<nTypes; i++) {
134       timeTrack = track->GetIntegratedTime(codes[i]);
135       if (abs(pdg) == codes[i]) hDiffSame[i]->Fill(timeTrue - timeTrack);
136     }
137   }  
138   
139   TCanvas *c = new TCanvas();
140   c->SetCanvasSize(640, 800);  
141   c->SetWindowSize(650, 830);  
142   c->Divide(2,3);
143  
144   c->cd(1);
145   hLength->Draw();
146   hLength->SetXTitle("True Length - Measured Length [cm]");
147   
148   for(Int_t i=0; i<nTypes; i++) {
149     c->cd(2+i);
150     gPad->SetGridx();
151     gPad->SetGridy();
152     hDiffSame[i]->Draw();  
153     hDiffSame[i]->SetXTitle("time true - measured [ps]");    
154     hDiffSame[i]->Fit("gaus");
155   }
156
157   char outFileName[80];
158   const char *refPlaneNames[] = {"inTPC", "outTPC","outTRD"}; 
159   
160   c->Modified();
161   c->ForceUpdate();
162
163   sprintf(outFileName, "time_%s_%3.2f-%3.2f.ps", refPlaneNames[nRefPlane-1], minP, maxP);
164   c->Print(outFileName);
165
166   if (treeRef) delete treeRef;
167   if (tracks) delete tracks;
168   
169 }