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