[u/mrichter/AliRoot.git] / TPC / AliL1Delay.C

/// \file AliL1Delay.C
/// \brief Macro that measures the time delay of L1 trigger in TPC
///
/// It reads the ESD tracks and fills histograms with the distance
/// between the primary vertex and the Z position
/// of the point of closest approach between the track and the beam axis.
/// The processing of forward and backward tracks is done separately.
/// The corresponding histograms are fitted with gaussian and the
/// shifts of the Z positions are extracted.
/// The L1 time delay is taken as an average of the forward and backward
/// (with opposite sign) shifts divided by the drift velocity.
/// We assume that the ESD tracks are obtained by the TPC tracking stand
/// alone.
/// The propagation to the vertex can be done either with or without
/// TGeoManager.
///
/// \author Cvetan Cheshkov <Cvetan.Cheshkov@cern.ch>

#if !defined(__CINT__) || defined(__MAKECINT__)
  #include <TMath.h>
  #include <TError.h>
  #include <Riostream.h>
  #include <TH1F.h>
  #include <TF1.h>
  #include <TTree.h>
  #include <TCanvas.h>
  #include <TFile.h>
  #include <TStyle.h>
  #include <TGeoManager.h>
  #include "AliRunLoader.h"
  #include "AliRun.h"
  #include "AliESD.h"
  #include "AliTracker.h"
  #include "AliTPCParam.h"
  #include "AliMagF.h"
  #include "AliITStrackV2.h"
#endif

Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t);
Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ);
Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ);

void AliL1Delay(const Char_t *esdfilename = "./AliESDs.root", const Char_t *galicefilename = "./galice.root", Bool_t withTGeo = kFALSE, const Char_t *geomfilename = "./geometry.root")
{
  const Int_t kMinNClusters = 0;
  const Double_t kPtMin = 0.5;
  const Double_t kMaxImpact = 0.3;
  const Double_t kDeltaZRes = 0.3;
  const Double_t kZRange = 20.;

  // Book histograms
  Int_t NBins = (Int_t)(2.*kZRange/(0.5*kDeltaZRes));
  TH1F *hDeltaZForward = new TH1F("hDeltaZForward","Longitudinal distance to the vertex (forward tracks)",NBins,-kZRange,kZRange); 
  TH1F *hDeltaZBackward = new TH1F("hDeltaZBackward","Longitudinal distance to the vertex (backward tracks)",NBins,-kZRange,kZRange);
  TH1F *hTrImpact = new TH1F("hTrImpact","Transverse impact parameter (all tracks)",100,0,30.*kMaxImpact); 

  // Open RunLoader
  AliRunLoader *rl = AliRunLoader::Open(galicefilename);
  if (!rl) {
    ::Error("AliL1Delay.C","Can't start session !");
    return;
  }
  rl->LoadgAlice();
  gAlice = rl->GetAliRun();

  // Set magnetic field
  AliMagF* field = TGeoGlobalMagField::Instance()->GetField();
  AliExternalTrackParam::SetNonuniformFieldTracking();
  const Float_t sfield = field->SolenoidField();

  // GetGeo manager
  if(withTGeo)
    {
      TFile *geoFile = TFile::Open(geomfilename);
      if (!geoFile || !geoFile->IsOpen()) {
	::Error("AliL1Delay.C","Can't open %s !",geomfilename);
	return;
      }
      gGeoManager = (TGeoManager *)geoFile->Get("Geo");
      if(!gGeoManager) {
	::Error("AliL1Delay.C","Can't find TGeoManager !");
	return;
      }
    }

  // Get X positions of the pad-rows
  TDirectory* saveDir = gDirectory;
  rl->CdGAFile();
  AliTPCParam* param = (AliTPCParam*) gDirectory->Get("75x40_100x60_150x60");
  if (!param) {
    ::Error("AliL1Delay.C","No TPC parameters found !");
    return;
  }
  Int_t nRowLow = param->GetNRowLow();
  Int_t nRowUp = param->GetNRowUp();
  Float_t xRow[159];
  for (Int_t i=0;i<nRowLow;i++){
    xRow[i] = param->GetPadRowRadiiLow(i);
  }
  for (Int_t i=0;i<nRowUp;i++){
    xRow[i+nRowLow] = param->GetPadRowRadiiUp(i);
  }
  saveDir->cd();

  // Open file with ESDs
  TFile *esdFile = TFile::Open(esdfilename);
  if (!esdFile || !esdFile->IsOpen()) {
    ::Error("AliL1Delay.C","Can't open %s !",esdfilename);
    return;
  }
  AliESD* event = new AliESD;
  TTree* esdTree = (TTree*) esdFile->Get("esdTree");
  if (!esdTree) {
    ::Error("AliL1Delay.C","No ESD tree found !");
    return;
  }
  esdTree->SetBranchAddress("ESD", &event);

  // Init PID
  AliPID pid;

  Int_t nEvent = esdTree->GetEntries();

  // Loop over events
  for(Int_t iEvent = 0; iEvent<nEvent; iEvent++)
    {
      esdTree->GetEvent(iEvent);
      const AliESDVertex *esdVtx = event->GetVertex();
      Double_t zVtx = esdVtx->GetZ();
      if(zVtx == 0.) continue; // Vertex is not found. Skip the event.
      Int_t nTrk=event->GetNumberOfTracks();

      ::Info("AliL1Delay.C","Event %d : Z vertex = %f  |  Number of tracks = %d",iEvent,zVtx,nTrk);

      for(Int_t iTrk=0; iTrk<nTrk; iTrk++)
	{
	  AliESDtrack *track = event->GetTrack(iTrk);
	  if(!track) continue;

	  // Track selection 
	  if((track->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
	  if(track->GetTPCclusters((Int_t)0x0)<(kMinNClusters*3)) continue;
	  if(track->GetP()<kPtMin) continue;

	  Int_t firstRow = (Int_t)(track->GetTPCPoints(0)+0.5);
	  Int_t lastRow = (Int_t)(track->GetTPCPoints(2)+0.5);
	  Double_t firstXYZ[3];
	  Double_t lastXYZ[3];
	  track->GetXYZAt(xRow[firstRow], sfield, firstXYZ);
	  track->GetXYZAt(xRow[lastRow], sfield, lastXYZ);

	  // Select if the track is forward or backward one
	  Bool_t IsForward;
	  if(firstXYZ[2] > 0 && lastXYZ[2] > 0)
	    IsForward = kTRUE;
	  else if(firstXYZ[2] < 0 && lastXYZ[2] < 0)
	    IsForward = kFALSE;
	  else
	    continue;

	  // Propagate the track to the vertex
	  Float_t vtxXYZ[3];
	  if(withTGeo) {
	    if(PropagateToVertexG(track,vtxXYZ) == 0) continue;
	  }
	  else {
	    if(PropagateToVertex(track,vtxXYZ) == 0) continue;
	  }

	  // Get the position at point of closest approach to the vertex
	  Float_t impact = TMath::Sqrt(vtxXYZ[0]*vtxXYZ[0]+vtxXYZ[1]*vtxXYZ[1]);
	  hTrImpact->Fill(impact);
	  // Select primary tracks
	  if(impact > kMaxImpact) continue;

	  if (IsForward)
	    hDeltaZForward->Fill(vtxXYZ[2]-zVtx);
	  else
	    hDeltaZBackward->Fill(vtxXYZ[2]-zVtx);
	}
    }

  //  delete event;
  //  esdFile->Close();

  // Check track propagation
  TF1 *fexpo = new TF1("fexpo","expo");
  hTrImpact->Fit("fexpo","E");
  Double_t slope = fexpo->GetParameter(1);
  if(slope > 3.*kMaxImpact)
    ::Warning("AliL1Delay.C","The transverse impact parameter distribution is too broad: %f +- %f",slope,fexpo->GetParError(1));
  delete fexpo;

  // Fit histograms and extract L1 time delay
  Double_t shifts[2],shifterrs2[2];
  {
    Double_t params[3];
    params[0] = hDeltaZForward->GetMaximumBin();
    params[1] = hDeltaZForward->GetBinCenter(hDeltaZForward->GetMaximumBin());
    params[2] = kDeltaZRes;
    TF1 *fgaus = new TF1("fgaus","gaus");
    fgaus->SetParameters(params);
    hDeltaZForward->Fit("fgaus","E");
    shifts[0] = fgaus->GetParameter(1);
    shifterrs2[0] = fgaus->GetParError(1)*fgaus->GetParError(1);
    delete fgaus;
  }
  {
    Double_t params[3];
    params[0] = hDeltaZBackward->GetMaximumBin();
    params[1] = hDeltaZBackward->GetBinCenter(hDeltaZBackward->GetMaximumBin());
    params[2] = kDeltaZRes;
    TF1 *fgaus = new TF1("fgaus","gaus");
    fgaus->SetParameters(params);
    hDeltaZBackward->Fit("fgaus","E");
    shifts[1] = -fgaus->GetParameter(1);
    shifterrs2[1] = fgaus->GetParError(1)*fgaus->GetParError(1);
    delete fgaus;
  }

  // Check the consistency of the two time delays
  if(TMath::Abs(shifts[0]-shifts[1])>3.*TMath::Sqrt(shifterrs2[0]+shifterrs2[1]))
    ::Warning("AliL1Delay.C","The extracted delays are more than 3 sigma away from each other: %f %f",shifts[0],shifts[1]);

  // Calculated the overall time delay
  Double_t shifterr = 1./TMath::Sqrt(1./shifterrs2[0]+1./shifterrs2[1]);
  Double_t shift = (shifts[0]/shifterrs2[0]+shifts[1]/shifterrs2[1])*shifterr*shifterr;

  ::Info("AliL1Delay.C","");
  ::Info("AliL1Delay.C","=====================================================");
  ::Info("AliL1Delay.C","The measured L1 delay is (%f +- %f) cm",shift,shifterr);
  ::Info("AliL1Delay.C","The measured L1 delay is (%f +- %f) ns",shift*1.e9/param->GetDriftV(),shifterr*1.e9/param->GetDriftV());
  ::Info("AliL1Delay.C","=====================================================");

  delete event;
  esdFile->Close();

  gStyle->SetOptFit(1);

  TCanvas *c1 = new TCanvas("c1","",0,0,700,850);

  TPad *p1 = new TPad("p1","",0,0,1,0.5); p1->Draw();
  p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
  p1->SetLogy(1);
  hTrImpact->SetFillColor(4);  hTrImpact->SetXTitle("(cm)");
  hTrImpact->SetStats(kTRUE); hTrImpact->Draw(); c1->cd();

  TPad *p2 = new TPad("p2","",0,0.5,0.5,1); p2->Draw();
  p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10); 
  hDeltaZForward->SetFillColor(4);  hDeltaZForward->SetXTitle("(cm)");
  hDeltaZForward->SetStats(kTRUE); hDeltaZForward->Draw(); c1->cd();

  TPad *p3 = new TPad("p3","",0.5,0.5,1,1); p3->Draw();
  p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10); 
  hDeltaZBackward->SetFillColor(4);  hDeltaZBackward->SetXTitle("(cm)");
  hDeltaZBackward->SetStats(kTRUE); hDeltaZBackward->Draw(); c1->cd();

  TFile f("AliL1Delay.root","RECREATE");
  c1->Write();
  f.Close();

  delete rl;
}

Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) {
  /// Correction for the material between the TPC and the ITS
  /// (should it belong to the TPC code ?)

  Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ? 
  Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
  Double_t yr=12.8, dr=0.03; // rods ?
  Double_t zm=0.2, dm=0.40;  // membrane
  //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
  Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)

  if (t->GetX() > riw) {
     if (!t->PropagateTo(riw,diw,x0iw)) return 1;
     if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
     if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
     if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
     //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
     //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r); 
     if (!t->PropagateTo(rs,ds)) return 1;
  } else if (t->GetX() < rs) {
     if (!t->PropagateTo(rs,-ds)) return 1;
     //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
     //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r); 
     if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
     if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
  } else {
  ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
    return 1;
  }
  
  return 0;
}

Bool_t PropagateToVertex(AliESDtrack *track, Float_t *vtxXYZ)
{
  /// Make an ITS track and propagate it to the vertex

  AliITStrackV2 itstrack(*track);
  if (CorrectForDeadZoneMaterial(&itstrack) != 0) return 0;
  if (!itstrack.PropagateTo(3.,0.0028,65.19)) return 0;
  if (!itstrack.PropagateToVertex()) return 0;

  itstrack.GetXYZ(vtxXYZ);

  return 1;
}

Bool_t PropagateToVertexG(AliESDtrack *track, Float_t *vtxXYZ)
{
  /// Make an ITS track and propagate it to the vertex using TGeoManager

  AliITStrackV2 itstrack(*track);
  AliExternalTrackParam etrack(itstrack);
  Double_t r = 3.;
  Double_t xyz0[3],xyz1[3],y,z;
  Double_t param[7];
  Double_t step = 2.;
  for (Double_t x = etrack.X()-step; x > r; x -= step) {
    etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
    Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
    etrack.RotateTo(alpha);
    etrack.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);      
    if (!etrack.GetProlongationAt(x,y,z)) return 0;
    xyz1[0] = x*TMath::Cos(alpha)-y*TMath::Sin(alpha);
    xyz1[1] = x*TMath::Sin(alpha)+y*TMath::Cos(alpha);
    xyz1[2] = z;
    AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
    if (!etrack.PropagateTo(x,param[1],param[0])) return 0;
  }
  etrack.PropagateTo(r,1,0);
  if (!etrack.PropagateToDCA(0,0,param[1],0)) return 0;

  etrack.GetXYZ(vtxXYZ);

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