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