from Alex Kalweit:
[u/mrichter/AliRoot.git] / TPC / AliL1Delay.C
1
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 "AliMagF.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 = TGeoGlobalMagField::Instance()->GetField();
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;
334     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
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 }