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