2 /////////////////////////////////////////////////////////////////////////////
4 // This is a macro which measures the time delay of L1 trigger in TPC //
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 //
16 // The propagation to the vertex can be done either with or without //
19 // Cvetan.Cheshkov@cern.ch //
21 /////////////////////////////////////////////////////////////////////////////
23 #if !defined(__CINT__) || defined(__MAKECINT__)
26 #include <Riostream.h>
33 #include <TGeoManager.h>
34 #include "AliRunLoader.h"
37 #include "AliTracker.h"
38 #include "AliTPCParam.h"
39 #include "AliMagFMaps.h"
40 #include "AliITStrackV2.h"
43 Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t);
44 Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ);
45 Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ);
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")
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.;
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);
62 AliRunLoader *rl = AliRunLoader::Open(galicefilename);
64 ::Error("AliL1Delay.C","Can't start session !");
68 gAlice = rl->GetAliRun();
71 AliMagF* field = gAlice->Field();
72 AliTracker::SetFieldMap(field,kTRUE);
73 AliExternalTrackParam::SetFieldMap(field);
74 AliExternalTrackParam::SetNonuniformFieldTracking();
75 const Float_t sfield = field->SolenoidField();
80 TFile *geoFile = TFile::Open(geomfilename);
81 if (!geoFile || !geoFile->IsOpen()) {
82 ::Error("AliL1Delay.C","Can't open %s !",geomfilename);
85 gGeoManager = (TGeoManager *)geoFile->Get("Geo");
87 ::Error("AliL1Delay.C","Can't find TGeoManager !");
92 // Get X positions of the pad-rows
93 TDirectory* saveDir = gDirectory;
95 AliTPCParam* param = (AliTPCParam*) gDirectory->Get("75x40_100x60_150x60");
97 ::Error("AliL1Delay.C","No TPC parameters found !");
100 Int_t nRowLow = param->GetNRowLow();
101 Int_t nRowUp = param->GetNRowUp();
103 for (Int_t i=0;i<nRowLow;i++){
104 xRow[i] = param->GetPadRowRadiiLow(i);
106 for (Int_t i=0;i<nRowUp;i++){
107 xRow[i+nRowLow] = param->GetPadRowRadiiUp(i);
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);
117 AliESD* event = new AliESD;
118 TTree* esdTree = (TTree*) esdFile->Get("esdTree");
120 ::Error("AliL1Delay.C","No ESD tree found !");
123 esdTree->SetBranchAddress("ESD", &event);
128 Int_t nEvent = esdTree->GetEntries();
131 for(Int_t iEvent = 0; iEvent<nEvent; iEvent++)
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();
139 ::Info("AliL1Delay.C","Event %d : Z vertex = %f | Number of tracks = %d",iEvent,zVtx,nTrk);
141 for(Int_t iTrk=0; iTrk<nTrk; iTrk++)
143 AliESDtrack *track = event->GetTrack(iTrk);
147 if((track->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
148 if(track->GetTPCclusters((Int_t)0x0)<(kMinNClusters*3)) continue;
149 if(track->GetP()<kPtMin) continue;
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];
155 track->GetXYZAt(xRow[firstRow], sfield, firstXYZ);
156 track->GetXYZAt(xRow[lastRow], sfield, lastXYZ);
158 // Select if the track is forward or backward one
160 if(firstXYZ[2] > 0 && lastXYZ[2] > 0)
162 else if(firstXYZ[2] < 0 && lastXYZ[2] < 0)
167 // Propagate the track to the vertex
170 if(PropagateToVertexG(track,vtxXYZ) == 0) continue;
173 if(PropagateToVertex(track,vtxXYZ) == 0) continue;
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;
183 hDeltaZForward->Fill(vtxXYZ[2]-zVtx);
185 hDeltaZBackward->Fill(vtxXYZ[2]-zVtx);
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));
200 // Fit histograms and extract L1 time delay
201 Double_t shifts[2],shifterrs2[2];
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);
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);
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]);
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;
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","=====================================================");
244 gStyle->SetOptFit(1);
246 TCanvas *c1 = new TCanvas("c1","",0,0,700,850);
248 TPad *p1 = new TPad("p1","",0,0,1,0.5); p1->Draw();
249 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
251 hTrImpact->SetFillColor(4); hTrImpact->SetXTitle("(cm)");
252 hTrImpact->SetStats(kTRUE); hTrImpact->Draw(); c1->cd();
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();
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();
264 TFile f("AliL1Delay.root","RECREATE");
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 ?)
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;
298 ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
305 Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ)
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;
313 itstrack.GetXYZ(vtxXYZ);
318 Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ)
320 // Make an ITS track and propagate it to the vertex using TGeoManager
321 AliITStrackV2 itstrack(*track);
322 AliExternalTrackParam etrack(itstrack);
324 Double_t xyz0[3],xyz1[3],y,z;
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);
336 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
337 if (!etrack.PropagateTo(x,param[1],param[0])) return 0;
339 etrack.PropagateTo(r,1,0);
340 if (!etrack.PropagateToDCA(0,0,param[1],0)) return 0;
342 etrack.GetXYZ(vtxXYZ);