ALIROOT-5769: Better handling of the error (R. Preghenella)
[u/mrichter/AliRoot.git] / TPC / AliL1Delay.C
CommitLineData
a7a1dd76 1/// \file AliL1Delay.C
2/// \brief Macro that measures the time delay of L1 trigger in TPC
3///
4/// It reads the ESD tracks and fills histograms with the distance
5/// between the primary vertex and the Z position
6/// of the point of closest approach between the track and the beam axis.
7/// The processing of forward and backward tracks is done separately.
8/// The corresponding histograms are fitted with gaussian and the
9/// shifts of the Z positions are extracted.
10/// The L1 time delay is taken as an average of the forward and backward
11/// (with opposite sign) shifts divided by the drift velocity.
12/// We assume that the ESD tracks are obtained by the TPC tracking stand
13/// alone.
14/// The propagation to the vertex can be done either with or without
15/// TGeoManager.
16///
17/// \author Cvetan Cheshkov <Cvetan.Cheshkov@cern.ch>
35d4a90a 18
19#if !defined(__CINT__) || defined(__MAKECINT__)
20 #include <TMath.h>
21 #include <TError.h>
22 #include <Riostream.h>
23 #include <TH1F.h>
24 #include <TF1.h>
25 #include <TTree.h>
26 #include <TCanvas.h>
27 #include <TFile.h>
28 #include <TStyle.h>
29 #include <TGeoManager.h>
30 #include "AliRunLoader.h"
31 #include "AliRun.h"
32 #include "AliESD.h"
33 #include "AliTracker.h"
34 #include "AliTPCParam.h"
f7a1cc68 35 #include "AliMagF.h"
35d4a90a 36 #include "AliITStrackV2.h"
37#endif
38
39Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t);
40Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ);
41Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ);
42
43void AliL1Delay(const Char_t *esdfilename = "./AliESDs.root", const Char_t *galicefilename = "./galice.root", Bool_t withTGeo = kFALSE, const Char_t *geomfilename = "./geometry.root")
44{
45 const Int_t kMinNClusters = 0;
46 const Double_t kPtMin = 0.5;
47 const Double_t kMaxImpact = 0.3;
48 const Double_t kDeltaZRes = 0.3;
49 const Double_t kZRange = 20.;
50
51 // Book histograms
52 Int_t NBins = (Int_t)(2.*kZRange/(0.5*kDeltaZRes));
53 TH1F *hDeltaZForward = new TH1F("hDeltaZForward","Longitudinal distance to the vertex (forward tracks)",NBins,-kZRange,kZRange);
54 TH1F *hDeltaZBackward = new TH1F("hDeltaZBackward","Longitudinal distance to the vertex (backward tracks)",NBins,-kZRange,kZRange);
55 TH1F *hTrImpact = new TH1F("hTrImpact","Transverse impact parameter (all tracks)",100,0,30.*kMaxImpact);
56
57 // Open RunLoader
58 AliRunLoader *rl = AliRunLoader::Open(galicefilename);
59 if (!rl) {
60 ::Error("AliL1Delay.C","Can't start session !");
61 return;
62 }
63 rl->LoadgAlice();
64 gAlice = rl->GetAliRun();
65
66 // Set magnetic field
f7a1cc68 67 AliMagF* field = TGeoGlobalMagField::Instance()->GetField();
35d4a90a 68 AliExternalTrackParam::SetNonuniformFieldTracking();
69 const Float_t sfield = field->SolenoidField();
70
71 // GetGeo manager
72 if(withTGeo)
73 {
74 TFile *geoFile = TFile::Open(geomfilename);
75 if (!geoFile || !geoFile->IsOpen()) {
76 ::Error("AliL1Delay.C","Can't open %s !",geomfilename);
77 return;
78 }
79 gGeoManager = (TGeoManager *)geoFile->Get("Geo");
80 if(!gGeoManager) {
81 ::Error("AliL1Delay.C","Can't find TGeoManager !");
82 return;
83 }
84 }
85
86 // Get X positions of the pad-rows
87 TDirectory* saveDir = gDirectory;
88 rl->CdGAFile();
89 AliTPCParam* param = (AliTPCParam*) gDirectory->Get("75x40_100x60_150x60");
90 if (!param) {
91 ::Error("AliL1Delay.C","No TPC parameters found !");
92 return;
93 }
94 Int_t nRowLow = param->GetNRowLow();
95 Int_t nRowUp = param->GetNRowUp();
96 Float_t xRow[159];
97 for (Int_t i=0;i<nRowLow;i++){
98 xRow[i] = param->GetPadRowRadiiLow(i);
99 }
100 for (Int_t i=0;i<nRowUp;i++){
101 xRow[i+nRowLow] = param->GetPadRowRadiiUp(i);
102 }
103 saveDir->cd();
104
105 // Open file with ESDs
106 TFile *esdFile = TFile::Open(esdfilename);
107 if (!esdFile || !esdFile->IsOpen()) {
108 ::Error("AliL1Delay.C","Can't open %s !",esdfilename);
109 return;
110 }
111 AliESD* event = new AliESD;
112 TTree* esdTree = (TTree*) esdFile->Get("esdTree");
113 if (!esdTree) {
114 ::Error("AliL1Delay.C","No ESD tree found !");
115 return;
116 }
117 esdTree->SetBranchAddress("ESD", &event);
118
119 // Init PID
120 AliPID pid;
121
122 Int_t nEvent = esdTree->GetEntries();
123
124 // Loop over events
125 for(Int_t iEvent = 0; iEvent<nEvent; iEvent++)
126 {
127 esdTree->GetEvent(iEvent);
128 const AliESDVertex *esdVtx = event->GetVertex();
e690d4d0 129 Double_t zVtx = esdVtx->GetZ();
35d4a90a 130 if(zVtx == 0.) continue; // Vertex is not found. Skip the event.
131 Int_t nTrk=event->GetNumberOfTracks();
132
133 ::Info("AliL1Delay.C","Event %d : Z vertex = %f | Number of tracks = %d",iEvent,zVtx,nTrk);
134
135 for(Int_t iTrk=0; iTrk<nTrk; iTrk++)
136 {
137 AliESDtrack *track = event->GetTrack(iTrk);
138 if(!track) continue;
139
140 // Track selection
141 if((track->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
142 if(track->GetTPCclusters((Int_t)0x0)<(kMinNClusters*3)) continue;
143 if(track->GetP()<kPtMin) continue;
144
145 Int_t firstRow = (Int_t)(track->GetTPCPoints(0)+0.5);
146 Int_t lastRow = (Int_t)(track->GetTPCPoints(2)+0.5);
147 Double_t firstXYZ[3];
148 Double_t lastXYZ[3];
149 track->GetXYZAt(xRow[firstRow], sfield, firstXYZ);
150 track->GetXYZAt(xRow[lastRow], sfield, lastXYZ);
151
152 // Select if the track is forward or backward one
153 Bool_t IsForward;
154 if(firstXYZ[2] > 0 && lastXYZ[2] > 0)
155 IsForward = kTRUE;
156 else if(firstXYZ[2] < 0 && lastXYZ[2] < 0)
157 IsForward = kFALSE;
158 else
159 continue;
160
161 // Propagate the track to the vertex
162 Float_t vtxXYZ[3];
163 if(withTGeo) {
164 if(PropagateToVertexG(track,vtxXYZ) == 0) continue;
165 }
166 else {
167 if(PropagateToVertex(track,vtxXYZ) == 0) continue;
168 }
169
170 // Get the position at point of closest approach to the vertex
171 Float_t impact = TMath::Sqrt(vtxXYZ[0]*vtxXYZ[0]+vtxXYZ[1]*vtxXYZ[1]);
172 hTrImpact->Fill(impact);
173 // Select primary tracks
174 if(impact > kMaxImpact) continue;
175
176 if (IsForward)
177 hDeltaZForward->Fill(vtxXYZ[2]-zVtx);
178 else
179 hDeltaZBackward->Fill(vtxXYZ[2]-zVtx);
180 }
181 }
182
183 // delete event;
184 // esdFile->Close();
185
186 // Check track propagation
187 TF1 *fexpo = new TF1("fexpo","expo");
188 hTrImpact->Fit("fexpo","E");
189 Double_t slope = fexpo->GetParameter(1);
190 if(slope > 3.*kMaxImpact)
191 ::Warning("AliL1Delay.C","The transverse impact parameter distribution is too broad: %f +- %f",slope,fexpo->GetParError(1));
192 delete fexpo;
193
194 // Fit histograms and extract L1 time delay
195 Double_t shifts[2],shifterrs2[2];
196 {
197 Double_t params[3];
198 params[0] = hDeltaZForward->GetMaximumBin();
199 params[1] = hDeltaZForward->GetBinCenter(hDeltaZForward->GetMaximumBin());
200 params[2] = kDeltaZRes;
201 TF1 *fgaus = new TF1("fgaus","gaus");
202 fgaus->SetParameters(params);
203 hDeltaZForward->Fit("fgaus","E");
204 shifts[0] = fgaus->GetParameter(1);
205 shifterrs2[0] = fgaus->GetParError(1)*fgaus->GetParError(1);
206 delete fgaus;
207 }
208 {
209 Double_t params[3];
210 params[0] = hDeltaZBackward->GetMaximumBin();
211 params[1] = hDeltaZBackward->GetBinCenter(hDeltaZBackward->GetMaximumBin());
212 params[2] = kDeltaZRes;
213 TF1 *fgaus = new TF1("fgaus","gaus");
214 fgaus->SetParameters(params);
215 hDeltaZBackward->Fit("fgaus","E");
216 shifts[1] = -fgaus->GetParameter(1);
217 shifterrs2[1] = fgaus->GetParError(1)*fgaus->GetParError(1);
218 delete fgaus;
219 }
220
221 // Check the consistency of the two time delays
222 if(TMath::Abs(shifts[0]-shifts[1])>3.*TMath::Sqrt(shifterrs2[0]+shifterrs2[1]))
223 ::Warning("AliL1Delay.C","The extracted delays are more than 3 sigma away from each other: %f %f",shifts[0],shifts[1]);
224
225 // Calculated the overall time delay
226 Double_t shifterr = 1./TMath::Sqrt(1./shifterrs2[0]+1./shifterrs2[1]);
227 Double_t shift = (shifts[0]/shifterrs2[0]+shifts[1]/shifterrs2[1])*shifterr*shifterr;
228
229 ::Info("AliL1Delay.C","");
230 ::Info("AliL1Delay.C","=====================================================");
231 ::Info("AliL1Delay.C","The measured L1 delay is (%f +- %f) cm",shift,shifterr);
232 ::Info("AliL1Delay.C","The measured L1 delay is (%f +- %f) ns",shift*1.e9/param->GetDriftV(),shifterr*1.e9/param->GetDriftV());
233 ::Info("AliL1Delay.C","=====================================================");
234
235 delete event;
236 esdFile->Close();
237
238 gStyle->SetOptFit(1);
239
240 TCanvas *c1 = new TCanvas("c1","",0,0,700,850);
241
242 TPad *p1 = new TPad("p1","",0,0,1,0.5); p1->Draw();
243 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
244 p1->SetLogy(1);
245 hTrImpact->SetFillColor(4); hTrImpact->SetXTitle("(cm)");
246 hTrImpact->SetStats(kTRUE); hTrImpact->Draw(); c1->cd();
247
248 TPad *p2 = new TPad("p2","",0,0.5,0.5,1); p2->Draw();
249 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
250 hDeltaZForward->SetFillColor(4); hDeltaZForward->SetXTitle("(cm)");
251 hDeltaZForward->SetStats(kTRUE); hDeltaZForward->Draw(); c1->cd();
252
253 TPad *p3 = new TPad("p3","",0.5,0.5,1,1); p3->Draw();
254 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
255 hDeltaZBackward->SetFillColor(4); hDeltaZBackward->SetXTitle("(cm)");
256 hDeltaZBackward->SetStats(kTRUE); hDeltaZBackward->Draw(); c1->cd();
257
258 TFile f("AliL1Delay.root","RECREATE");
259 c1->Write();
260 f.Close();
261
262 delete rl;
263}
264
265Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) {
266 //--------------------------------------------------------------------
267 // Correction for the material between the TPC and the ITS
268 // (should it belong to the TPC code ?)
269 //--------------------------------------------------------------------
270 Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
271 Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
272 Double_t yr=12.8, dr=0.03; // rods ?
273 Double_t zm=0.2, dm=0.40; // membrane
274 //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
275 Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
276
277 if (t->GetX() > riw) {
278 if (!t->PropagateTo(riw,diw,x0iw)) return 1;
279 if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
280 if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
281 if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
282 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
283 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
284 if (!t->PropagateTo(rs,ds)) return 1;
285 } else if (t->GetX() < rs) {
286 if (!t->PropagateTo(rs,-ds)) return 1;
287 //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
288 //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
289 if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
290 if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
291 } else {
292 ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
293 return 1;
294 }
295
296 return 0;
297}
298
299Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ)
300{
301 // Make an ITS track and propagate it to the vertex
302 AliITStrackV2 itstrack(*track);
303 if (CorrectForDeadZoneMaterial(&itstrack) != 0) return 0;
304 if (!itstrack.PropagateTo(3.,0.0028,65.19)) return 0;
305 if (!itstrack.PropagateToVertex()) return 0;
306
307 itstrack.GetXYZ(vtxXYZ);
308
309 return 1;
310}
311
312Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ)
313{
314 // Make an ITS track and propagate it to the vertex using TGeoManager
315 AliITStrackV2 itstrack(*track);
316 AliExternalTrackParam etrack(itstrack);
317 Double_t r = 3.;
318 Double_t xyz0[3],xyz1[3],y,z;
319 Double_t param[7];
320 Double_t step = 2.;
321 for (Double_t x = etrack.X()-step; x > r; x -= step) {
322 etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
323 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
324 etrack.RotateTo(alpha);
325 etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
326 if (!etrack.GetProlongationAt(x,y,z)) return 0;
327 xyz1[0] = x*TMath::Cos(alpha)-y*TMath::Sin(alpha);
328 xyz1[1] = x*TMath::Sin(alpha)+y*TMath::Cos(alpha);
329 xyz1[2] = z;
70ae41ca 330 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
35d4a90a 331 if (!etrack.PropagateTo(x,param[1],param[0])) return 0;
332 }
333 etrack.PropagateTo(r,1,0);
334 if (!etrack.PropagateToDCA(0,0,param[1],0)) return 0;
335
336 etrack.GetXYZ(vtxXYZ);
337
338 return 1;
339}