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