]>
Commit | Line | Data |
---|---|---|
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> | |
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" | |
35 | #include "AliMagF.h" | |
36 | #include "AliITStrackV2.h" | |
37 | #endif | |
38 | ||
39 | Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t); | |
40 | Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ); | |
41 | Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ); | |
42 | ||
43 | void 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 | |
67 | AliMagF* field = TGeoGlobalMagField::Instance()->GetField(); | |
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(); | |
129 | Double_t zVtx = esdVtx->GetZ(); | |
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 | ||
265 | Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) { | |
266 | /// Correction for the material between the TPC and the ITS | |
267 | /// (should it belong to the TPC code ?) | |
268 | ||
269 | Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ? | |
270 | Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ? | |
271 | Double_t yr=12.8, dr=0.03; // rods ? | |
272 | Double_t zm=0.2, dm=0.40; // membrane | |
273 | //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails | |
274 | Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?) | |
275 | ||
276 | if (t->GetX() > riw) { | |
277 | if (!t->PropagateTo(riw,diw,x0iw)) return 1; | |
278 | if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr); | |
279 | if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm); | |
280 | if (!t->PropagateTo(rcd,dcd,x0cd)) return 1; | |
281 | //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); | |
282 | //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r); | |
283 | if (!t->PropagateTo(rs,ds)) return 1; | |
284 | } else if (t->GetX() < rs) { | |
285 | if (!t->PropagateTo(rs,-ds)) return 1; | |
286 | //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z); | |
287 | //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r); | |
288 | if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1; | |
289 | if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1; | |
290 | } else { | |
291 | ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !"); | |
292 | return 1; | |
293 | } | |
294 | ||
295 | return 0; | |
296 | } | |
297 | ||
298 | Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ) | |
299 | { | |
300 | /// Make an ITS track and propagate it to the vertex | |
301 | ||
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 | ||
312 | Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ) | |
313 | { | |
314 | /// Make an ITS track and propagate it to the vertex using TGeoManager | |
315 | ||
316 | AliITStrackV2 itstrack(*track); | |
317 | AliExternalTrackParam etrack(itstrack); | |
318 | Double_t r = 3.; | |
319 | Double_t xyz0[3],xyz1[3],y,z; | |
320 | Double_t param[7]; | |
321 | Double_t step = 2.; | |
322 | for (Double_t x = etrack.X()-step; x > r; x -= step) { | |
323 | etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); | |
324 | Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]); | |
325 | etrack.RotateTo(alpha); | |
326 | etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); | |
327 | if (!etrack.GetProlongationAt(x,y,z)) return 0; | |
328 | xyz1[0] = x*TMath::Cos(alpha)-y*TMath::Sin(alpha); | |
329 | xyz1[1] = x*TMath::Sin(alpha)+y*TMath::Cos(alpha); | |
330 | xyz1[2] = z; | |
331 | AliTracker::MeanMaterialBudget(xyz0,xyz1,param); | |
332 | if (!etrack.PropagateTo(x,param[1],param[0])) return 0; | |
333 | } | |
334 | etrack.PropagateTo(r,1,0); | |
335 | if (!etrack.PropagateToDCA(0,0,param[1],0)) return 0; | |
336 | ||
337 | etrack.GetXYZ(vtxXYZ); | |
338 | ||
339 | return 1; | |
340 | } |