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 "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;
 336     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
 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 }