[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDtrackingResolution.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercialf purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

/* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Reconstruction QA                                                     //
//                                                                        //
//  Authors:                                                              //
//    Markus Fasel <M.Fasel@gsi.de>                                       //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <cstring>


#include <TObjArray.h>
#include <TList.h>
#include <TH2.h>
#include <TH1.h>
#include <TF1.h>
#include <TProfile.h>
#include <TGraphErrors.h>
#include <TMath.h>
#include "TTreeStream.h"
#include "TGeoManager.h"

#include "AliAnalysisManager.h"
#include "AliTrackReference.h"
#include "AliTrackPointArray.h"
#include "AliCDBManager.h"

#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDtrackerV1.h"
#include "AliTRDReconstructor.h"
#include "AliTRDrecoParam.h"

#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
#include "AliTRDtrackingResolution.h"

ClassImp(AliTRDtrackingResolution)

//________________________________________________________
AliTRDtrackingResolution::AliTRDtrackingResolution(const char * name):
  AliAnalysisTask(name, "")
  ,fTracks(0x0)
  ,fHistos(0x0)
  ,fReconstructor(0x0)
  ,fHasMCdata(kTRUE)
  ,fDebugLevel(0)
  ,fDebugStream(0x0)
{
  TGeoManager::Import("geometry.root");
  fReconstructor = new AliTRDReconstructor();
  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());

  DefineInput(0, TObjArray::Class());
  DefineOutput(0, TList::Class());
}

//________________________________________________________
AliTRDtrackingResolution::~AliTRDtrackingResolution()
{
  fHistos->Delete(); delete fHistos;
  delete fReconstructor;
}

//________________________________________________________
void AliTRDtrackingResolution::ConnectInputData(Option_t *){
  fTracks = dynamic_cast<TObjArray *>(GetInputData(0));
}

//________________________________________________________
void AliTRDtrackingResolution::CreateOutputObjects()
{
  // spatial resolution
  OpenFile(0, "RECREATE");

  fHistos = new TList();

  // Resolution histos
  if(HasMCdata()){
    // tracklet resolution [0]
    fHistos->AddAt(new TH2I("fY", "", 21, -21., 21., 100, -.5, .5), 0);
    // tracklet angular resolution [1]
    fHistos->AddAt(new TH2I("fPhi", "", 21, -21., 21., 100, -10., 10.), 1);
  }
  // Residual histos
  // cluster to tracklet residuals [2]
  Int_t position = HasMCdata() ? 2 : 0;
  fHistos->AddAt(new TH2I("fYClRes", "", 21, -21., 21., 100, -.5, .5), position);
}

//________________________________________________________
void AliTRDtrackingResolution::Exec(Option_t *)
{
  // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
  // angular Resolution: res = Tracklet angle - TrackRef Angle

  Int_t nTrackInfos = fTracks->GetEntriesFast();
  if(fDebugLevel>=2) printf("Number of Histograms: %d\n", fHistos->GetEntries());

  Double_t dy, dz;
  AliTrackPoint tr[kNLayers], tk[kNLayers];
  AliTRDtrackV1 *fTrack = 0x0;
  AliTRDtrackInfo *fInfo = 0x0;
  if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
  for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
    // check if ESD and MC-Information are available
    if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
    if(!(fTrack = fInfo->GetTRDtrack())) continue;
    if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());


    Int_t npts = 0;
    AliTRDseedV1 *fTracklet = 0x0;
    for(Int_t iplane = 0; iplane < kNLayers; iplane++){
      if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
      if(!fTracklet->IsOK()) continue;

      // Book point arrays
      tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
      tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
      npts++;

      if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());

      Float_t phi = TMath::ATan(fTracklet->GetYref(1));

      // RESOLUTION (compare to MC)
      if(HasMCdata()){
        if(fInfo->GetNTrackRefs() >= 2){ 
        Double_t phiMC;
        if(Resolution(fTracklet, fInfo, phiMC)) phi = phiMC;
        }
      }

      // Do clusters residuals
      if(!fTracklet->Fit(kFALSE)) continue;
      Int_t histpos = HasMCdata() ? 2 : 0;
      AliTRDcluster *c = 0x0;
      for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
        if(!(c = fTracklet->GetClusters(ic))) continue;
        
        dy = fTracklet->GetYat(c->GetX()) - c->GetY();
        ((TH2I*)fHistos->At(histpos))->Fill(phi*TMath::RadToDeg(), dy);
        if(fDebugLevel>=1){
          Float_t q = c->GetQ();
          (*fDebugStream) << "ClsTrkltResidual"
            << "plane="	 	<< iplane
            << "phi="     << phi
            << "q="       << q
            << "dy="      << dy
            << "\n";
        }
      }
    }


    // this protection we might drop TODO
    if(fTrack->GetNumberOfTracklets() < 6) continue;

    AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
    for(Int_t ip=0; ip<npts; ip++){
      dy = tk[ip].GetY() - tr[ip].GetY();
      dz = tk[ip].GetZ() - tr[ip].GetZ();
      if(fDebugLevel>=1){
        (*fDebugStream) << "ResidualsRT"
          << "dy="		  << dy
          << "dz="	 	  << dz
/*          << "phi="			<< phi
          << "theta="		<< theta
          << "dphi="		<< dphi*/
          << "\n";
      }
    }

//     AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
//     for(Int_t ip=0; ip<nc; ip++){
//       dy = cl[ip].GetY() - tr[ip].GetY();
//       dz = cl[ip].GetZ() - tr[ip].GetZ();
//       if(fDebugLevel>=1){
//         (*fDebugStream) << "ResidualsKF"
//           << "dy="		  << dy
//           << "dz="	 	  << dz
// /*          << "phi="			<< phi
//           << "theta="		<< theta
//           << "dphi="		<< dphi*/
//           << "\n";
//       }
//     }    


  }
  PostData(0, fHistos);
}


//________________________________________________________
Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &phi)
{

  AliTrackReference *fTrackRefs[2] = {0x0, 0x0},   *tempTrackRef = 0x0;

  // check for 2 track ref where the radial position has a distance less than 3.7mm
  Int_t nFound = 0;
  for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
    if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
    if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
    if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
    fTrackRefs[nFound++] = tempTrackRef;
    if(nFound == 2) break;
  }
  if(nFound < 2){ 
    if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
    return kFALSE;
  }
  // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation


  // RESOLUTION
  Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
  if(dx <= 0.){
    if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
    return kFALSE;
  }
  Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
  Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
  Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
  Double_t ymc =  fTrackRefs[1]->LocalY() - dydx*dx0;
  Double_t zmc =  fTrackRefs[1]->Z() - dzdx*dx0;
  
  // recalculate tracklet based on the MC info
  tracklet->SetZref(0, zmc);
  tracklet->SetZref(1, dzdx);
  tracklet->Fit();
  Double_t dy = tracklet->GetYfit(0) - ymc;
  Double_t dz = tracklet->GetZfit(0) - zmc;
      
  //res_y *= 100; // in mm
  Double_t momentum = fTrackRefs[0]->P();

  phi   = TMath::ATan(dydx);
  Double_t theta = TMath::ATan(dzdx);
  Double_t dphi   = TMath::ATan(tracklet->GetYfit(1)) - phi;
  if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
  
  // Fill Histograms
  if(TMath::Abs(dx-3.7)<1.E-3){
    ((TH2I*)fHistos->At(0))->Fill(phi*TMath::RadToDeg(), dy);
    ((TH2I*)fHistos->At(1))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
  }        
  // Fill Debug Tree
  if(fDebugLevel>=1){
    Int_t iplane = tracklet->GetPlane();
    (*fDebugStream) << "TrkltResolution"
      << "plane="	 	<< iplane
      << "p="       << momentum
      << "dx="      << dx
      << "dy="		  << dy
      << "dz="	 	  << dz
      << "phi="			<< phi
      << "theta="		<< theta
      << "dphi="		<< dphi
      << "ymc="     << ymc
      << "zmc="     << zmc
      << "tracklet.="<<tracklet 
      << "\n";
  }

  return kTRUE;
}


//________________________________________________________
void AliTRDtrackingResolution::Terminate(Option_t *)
{
  if(fDebugStream) delete fDebugStream;

  TH2I *h2 = 0x0;
  TH1D *h = 0x0;
  TF1 f("f1", "gaus", -.5, .5);  
  if(HasMCdata()){
    //process distributions
    fHistos = dynamic_cast<TList*>(GetOutputData(0));
    if (!fHistos) {
      Printf("ERROR: list not available");
      return;
    }

  // y resolution
    h2 = (TH2I*)fHistos->At(0);
    TGraphErrors *gm = new TGraphErrors(h2->GetNbinsX());
    gm->SetNameTitle("meany", "Mean dy");
    TGraphErrors *gs = new TGraphErrors(h2->GetNbinsX());
    gs->SetNameTitle("sigmy", "Sigma y");
    for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
      Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
      f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
      h = h2->ProjectionY("py", iphi, iphi);
      h->Fit(&f, "QN", "", -.5, .5);
      Int_t jphi = iphi -1;
      gm->SetPoint(jphi, phi, f.GetParameter(1));
      gm->SetPointError(jphi, 0., f.GetParError(1));
      gs->SetPoint(jphi, phi, f.GetParameter(2));
      gs->SetPointError(jphi, 0., f.GetParError(2));
    }
    fHistos->Add(gm);
    fHistos->Add(gs);
  
    // phi resolution
    h2 = (TH2I*)fHistos->At(1);
    gm = new TGraphErrors(h2->GetNbinsX());
    gm->SetNameTitle("meanphi", "Mean Phi");
    gs = new TGraphErrors(h2->GetNbinsX());
    gs->SetNameTitle("sigmphi", "Sigma Phi");
    for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
      Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
      f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
      h = h2->ProjectionY("py", iphi, iphi);
      h->Fit(&f, "QN", "", -.5, .5);
      Int_t jphi = iphi -1;
      gm->SetPoint(jphi, phi, f.GetParameter(1));
      gm->SetPointError(jphi, 0., f.GetParError(1));
      gs->SetPoint(jphi, phi, f.GetParameter(2));
      gs->SetPointError(jphi, 0., f.GetParError(2));
    }
    fHistos->Add(gm);
    fHistos->Add(gs);
  }

  // Fit clusters residuals
  Int_t position_residuals = fHasMCdata ? 2 : 0;
  h2 = (TH2I *)(fHistos->At(position_residuals));
  TGraphErrors *residuals_mean = new TGraphErrors(h2->GetNbinsX());
  residuals_mean->SetLineColor(kGreen);
  residuals_mean->SetMarkerStyle(22);
  residuals_mean->SetMarkerColor(kGreen);
  TGraphErrors *residuals_sigma = new TGraphErrors(h2->GetNbinsX());
  residuals_mean->SetNameTitle("residuals_mean", "Residuals Mean Phi");
  residuals_sigma->SetNameTitle("residuals_sigma", "Residuals Sigma Phi");
  residuals_sigma->SetLineColor(kRed);
  residuals_sigma->SetMarkerStyle(23);
  residuals_sigma->SetMarkerColor(kRed);
  for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
    Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
    Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
    h = h2->ProjectionY("py", ibin, ibin);
    h->Fit(&f, "QN", "", -0.5, 0.5);
    residuals_mean->SetPoint(ibin - 1, phi, f.GetParameter(1));
    residuals_mean->SetPointError(ibin - 1, dphi, f.GetParError(1));
    residuals_sigma->SetPoint(ibin - 1, phi, f.GetParameter(2));
    residuals_sigma->SetPointError(ibin - 1, dphi, f.GetParError(2));
  }
  fHistos->Add(residuals_mean);
  fHistos->Add(residuals_sigma);
}

//________________________________________________________
void AliTRDtrackingResolution::SetDebugLevel(Int_t level){
  fDebugLevel = level;
  if(!fDebugLevel) return;
  if(fDebugStream) return;
  fDebugStream = new TTreeSRedirector("TRD.Resolution.root");
}

//________________________________________________________
void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
{
  fReconstructor->SetRecoParam(r);
}
Commit	Line	Data
77203477	1	/**************************************************************************
d2381af5	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercialf purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
77203477	15
	16	/* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Reconstruction QA //
	21	// //
	22	// Authors: //
	23	// Markus Fasel <M.Fasel@gsi.de> //
	24	// //
	25	////////////////////////////////////////////////////////////////////////////
	26
aaf47b30	27	#include <cstring>
	28
	29
77203477	30	#include <TObjArray.h>
77203477	31	#include <TList.h>
7102d1b1	32	#include <TH2.h>
	33	#include <TH1.h>
	34	#include <TF1.h>
77203477	35	#include <TProfile.h>
7102d1b1	36	#include <TGraphErrors.h>
77203477	37	#include <TMath.h>
aaf47b30	38	#include "TTreeStream.h"
aaf47b30	39	#include "TGeoManager.h"
77203477	40
77203477	41	#include "AliAnalysisManager.h"
77203477	42	#include "AliTrackReference.h"
aaf47b30	43	#include "AliTrackPointArray.h"
	44	#include "AliCDBManager.h"
	45
	46	#include "AliTRDcluster.h"
	47	#include "AliTRDseedV1.h"
	48	#include "AliTRDtrackV1.h"
	49	#include "AliTRDtrackerV1.h"
	50	#include "AliTRDReconstructor.h"
	51	#include "AliTRDrecoParam.h"
77203477	52
	53	#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
	54	#include "AliTRDtrackingResolution.h"
	55
77203477	56	ClassImp(AliTRDtrackingResolution)
	57
	58	//________________________________________________________
	59	AliTRDtrackingResolution::AliTRDtrackingResolution(const char * name):
aaf47b30	60	AliAnalysisTask(name, "")
39779ce6	61	,fTracks(0x0)
39779ce6	62	,fHistos(0x0)
aaf47b30	63	,fReconstructor(0x0)
4b8f8a35	64	,fHasMCdata(kTRUE)
aaf47b30	65	,fDebugLevel(0)
aaf47b30	66	,fDebugStream(0x0)
77203477	67	{
aaf47b30	68	TGeoManager::Import("geometry.root");
	69	fReconstructor = new AliTRDReconstructor();
	70	fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
aaf47b30	71
77203477	72	DefineInput(0, TObjArray::Class());
	73	DefineOutput(0, TList::Class());
	74	}
	75
ed383798	76	//________________________________________________________
	77	AliTRDtrackingResolution::~AliTRDtrackingResolution()
	78	{
	79	fHistos->Delete(); delete fHistos;
	80	delete fReconstructor;
	81	}
	82
77203477	83	//________________________________________________________
77203477	84	void AliTRDtrackingResolution::ConnectInputData(Option_t *){
39779ce6	85	fTracks = dynamic_cast<TObjArray *>(GetInputData(0));
77203477	86	}
	87
	88	//________________________________________________________
	89	void AliTRDtrackingResolution::CreateOutputObjects()
	90	{
	91	// spatial resolution
77203477	92	OpenFile(0, "RECREATE");
39779ce6	93
	94	fHistos = new TList();
	95
	96	// Resolution histos
4b8f8a35	97	if(HasMCdata()){
d2381af5	98	// tracklet resolution [0]
	99	fHistos->AddAt(new TH2I("fY", "", 21, -21., 21., 100, -.5, .5), 0);
	100	// tracklet angular resolution [1]
	101	fHistos->AddAt(new TH2I("fPhi", "", 21, -21., 21., 100, -10., 10.), 1);
	102	}
39779ce6	103	// Residual histos
39779ce6	104	// cluster to tracklet residuals [2]
4b8f8a35	105	Int_t position = HasMCdata() ? 2 : 0;
4b8f8a35	106	fHistos->AddAt(new TH2I("fYClRes", "", 21, -21., 21., 100, -.5, .5), position);
77203477	107	}
	108
	109	//________________________________________________________
7102d1b1	110	void AliTRDtrackingResolution::Exec(Option_t *)
7102d1b1	111	{
77203477	112	// spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
77203477	113	// angular Resolution: res = Tracklet angle - TrackRef Angle
7102d1b1	114
39779ce6	115	Int_t nTrackInfos = fTracks->GetEntriesFast();
39779ce6	116	if(fDebugLevel>=2) printf("Number of Histograms: %d\n", fHistos->GetEntries());
aaf47b30	117
aaf47b30	118	Double_t dy, dz;
39779ce6	119	AliTrackPoint tr[kNLayers], tk[kNLayers];
aaf47b30	120	AliTRDtrackV1 *fTrack = 0x0;
77203477	121	AliTRDtrackInfo *fInfo = 0x0;
	122	if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
	123	for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
77203477	124	// check if ESD and MC-Information are available
39779ce6	125	if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
aaf47b30	126	if(!(fTrack = fInfo->GetTRDtrack())) continue;
7102d1b1	127	if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
7102d1b1	128
39779ce6	129
39779ce6	130	Int_t npts = 0;
aaf47b30	131	AliTRDseedV1 *fTracklet = 0x0;
77203477	132	for(Int_t iplane = 0; iplane < kNLayers; iplane++){
aaf47b30	133	if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
	134	if(!fTracklet->IsOK()) continue;
	135
39779ce6	136	// Book point arrays
	137	tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
	138	tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
	139	npts++;
	140
7102d1b1	141	if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
7102d1b1	142
39779ce6	143	Float_t phi = TMath::ATan(fTracklet->GetYref(1));
7102d1b1	144
39779ce6	145	// RESOLUTION (compare to MC)
4b8f8a35	146	if(HasMCdata()){
d2381af5	147	if(fInfo->GetNTrackRefs() >= 2){
84b3ecda	148	Double_t phiMC;
d2381af5	149	if(Resolution(fTracklet, fInfo, phiMC)) phi = phiMC;
	150	}
	151	}
aaf47b30	152
39779ce6	153	// Do clusters residuals
d2381af5	154	if(!fTracklet->Fit(kFALSE)) continue;
4b8f8a35	155	Int_t histpos = HasMCdata() ? 2 : 0;
aaf47b30	156	AliTRDcluster *c = 0x0;
	157	for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
	158	if(!(c = fTracklet->GetClusters(ic))) continue;
39779ce6	159
39779ce6	160	dy = fTracklet->GetYat(c->GetX()) - c->GetY();
4b8f8a35	161	((TH2I)fHistos->At(histpos))->Fill(phiTMath::RadToDeg(), dy);
39779ce6	162	if(fDebugLevel>=1){
	163	Float_t q = c->GetQ();
	164	(*fDebugStream) << "ClsTrkltResidual"
	165	<< "plane=" << iplane
	166	<< "phi=" << phi
	167	<< "q=" << q
	168	<< "dy=" << dy
	169	<< "\n";
	170	}
aaf47b30	171	}
aaf47b30	172	}
39779ce6	173
d2381af5	174
39779ce6	175	// this protection we might drop TODO
	176	if(fTrack->GetNumberOfTracklets() < 6) continue;
	177
	178	AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
	179	for(Int_t ip=0; ip<npts; ip++){
	180	dy = tk[ip].GetY() - tr[ip].GetY();
	181	dz = tk[ip].GetZ() - tr[ip].GetZ();
aaf47b30	182	if(fDebugLevel>=1){
	183	(*fDebugStream) << "ResidualsRT"
	184	<< "dy=" << dy
	185	<< "dz=" << dz
	186	/* << "phi=" << phi
	187	<< "theta=" << theta
	188	<< "dphi=" << dphi*/
	189	<< "\n";
	190	}
	191	}
	192
	193	// AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
	194	// for(Int_t ip=0; ip<nc; ip++){
	195	// dy = cl[ip].GetY() - tr[ip].GetY();
	196	// dz = cl[ip].GetZ() - tr[ip].GetZ();
	197	// if(fDebugLevel>=1){
	198	// (*fDebugStream) << "ResidualsKF"
	199	// << "dy=" << dy
	200	// << "dz=" << dz
	201	// /* << "phi=" << phi
	202	// << "theta=" << theta
	203	// << "dphi=" << dphi*/
	204	// << "\n";
	205	// }
	206	// }
39779ce6	207
39779ce6	208
77203477	209	}
39779ce6	210	PostData(0, fHistos);
77203477	211	}
77203477	212
39779ce6	213
39779ce6	214	//________________________________________________________
84b3ecda	215	Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 tracklet, AliTRDtrackInfo fInfo, Double_t &phi)
39779ce6	216	{
	217
	218	AliTrackReference fTrackRefs[2] = {0x0, 0x0}, tempTrackRef = 0x0;
	219
	220	// check for 2 track ref where the radial position has a distance less than 3.7mm
	221	Int_t nFound = 0;
	222	for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
	223	if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
	224	if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
	225	if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
	226	fTrackRefs[nFound++] = tempTrackRef;
	227	if(nFound == 2) break;
	228	}
	229	if(nFound < 2){
	230	if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
	231	return kFALSE;
	232	}
	233	// We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
	234
	235
	236	// RESOLUTION
	237	Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
	238	if(dx <= 0.){
	239	if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
	240	return kFALSE;
	241	}
	242	Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
	243	Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
	244	Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
	245	Double_t ymc = fTrackRefs[1]->LocalY() - dydx*dx0;
	246	Double_t zmc = fTrackRefs[1]->Z() - dzdx*dx0;
	247
	248	// recalculate tracklet based on the MC info
84b3ecda	249	tracklet->SetZref(0, zmc);
84b3ecda	250	tracklet->SetZref(1, dzdx);
39779ce6	251	tracklet->Fit();
	252	Double_t dy = tracklet->GetYfit(0) - ymc;
	253	Double_t dz = tracklet->GetZfit(0) - zmc;
	254
	255	//res_y *= 100; // in mm
	256	Double_t momentum = fTrackRefs[0]->P();
	257
	258	phi = TMath::ATan(dydx);
	259	Double_t theta = TMath::ATan(dzdx);
	260	Double_t dphi = TMath::ATan(tracklet->GetYfit(1)) - phi;
	261	if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
	262
	263	// Fill Histograms
	264	if(TMath::Abs(dx-3.7)<1.E-3){
	265	((TH2I)fHistos->At(0))->Fill(phiTMath::RadToDeg(), dy);
	266	((TH2I)fHistos->At(1))->Fill(phiTMath::RadToDeg(), dphi*TMath::RadToDeg());
	267	}
	268	// Fill Debug Tree
	269	if(fDebugLevel>=1){
	270	Int_t iplane = tracklet->GetPlane();
	271	(*fDebugStream) << "TrkltResolution"
	272	<< "plane=" << iplane
	273	<< "p=" << momentum
	274	<< "dx=" << dx
	275	<< "dy=" << dy
	276	<< "dz=" << dz
	277	<< "phi=" << phi
	278	<< "theta=" << theta
	279	<< "dphi=" << dphi
84b3ecda	280	<< "ymc=" << ymc
	281	<< "zmc=" << zmc
	282	<< "tracklet.="<<tracklet
39779ce6	283	<< "\n";
	284	}
	285
	286	return kTRUE;
	287	}
	288
	289
77203477	290	//________________________________________________________
7102d1b1	291	void AliTRDtrackingResolution::Terminate(Option_t *)
7102d1b1	292	{
77203477	293	if(fDebugStream) delete fDebugStream;
7102d1b1	294
d2381af5	295	TH2I *h2 = 0x0;
	296	TH1D *h = 0x0;
	297	TF1 f("f1", "gaus", -.5, .5);
	298	if(HasMCdata()){
	299	//process distributions
	300	fHistos = dynamic_cast<TList*>(GetOutputData(0));
	301	if (!fHistos) {
	302	Printf("ERROR: list not available");
	303	return;
	304	}
	305
	306	// y resolution
	307	h2 = (TH2I*)fHistos->At(0);
	308	TGraphErrors *gm = new TGraphErrors(h2->GetNbinsX());
	309	gm->SetNameTitle("meany", "Mean dy");
	310	TGraphErrors *gs = new TGraphErrors(h2->GetNbinsX());
	311	gs->SetNameTitle("sigmy", "Sigma y");
	312	for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
	313	Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
	314	f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
	315	h = h2->ProjectionY("py", iphi, iphi);
	316	h->Fit(&f, "QN", "", -.5, .5);
	317	Int_t jphi = iphi -1;
	318	gm->SetPoint(jphi, phi, f.GetParameter(1));
	319	gm->SetPointError(jphi, 0., f.GetParError(1));
	320	gs->SetPoint(jphi, phi, f.GetParameter(2));
	321	gs->SetPointError(jphi, 0., f.GetParError(2));
	322	}
	323	fHistos->Add(gm);
	324	fHistos->Add(gs);
	325
	326	// phi resolution
	327	h2 = (TH2I*)fHistos->At(1);
	328	gm = new TGraphErrors(h2->GetNbinsX());
	329	gm->SetNameTitle("meanphi", "Mean Phi");
	330	gs = new TGraphErrors(h2->GetNbinsX());
	331	gs->SetNameTitle("sigmphi", "Sigma Phi");
	332	for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
	333	Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
	334	f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
	335	h = h2->ProjectionY("py", iphi, iphi);
	336	h->Fit(&f, "QN", "", -.5, .5);
	337	Int_t jphi = iphi -1;
	338	gm->SetPoint(jphi, phi, f.GetParameter(1));
	339	gm->SetPointError(jphi, 0., f.GetParError(1));
	340	gs->SetPoint(jphi, phi, f.GetParameter(2));
	341	gs->SetPointError(jphi, 0., f.GetParError(2));
	342	}
	343	fHistos->Add(gm);
	344	fHistos->Add(gs);
	345	}
	346
	347	// Fit clusters residuals
	348	Int_t position_residuals = fHasMCdata ? 2 : 0;
	349	h2 = (TH2I *)(fHistos->At(position_residuals));
	350	TGraphErrors *residuals_mean = new TGraphErrors(h2->GetNbinsX());
	351	residuals_mean->SetLineColor(kGreen);
	352	residuals_mean->SetMarkerStyle(22);
	353	residuals_mean->SetMarkerColor(kGreen);
	354	TGraphErrors *residuals_sigma = new TGraphErrors(h2->GetNbinsX());
	355	residuals_mean->SetNameTitle("residuals_mean", "Residuals Mean Phi");
	356	residuals_sigma->SetNameTitle("residuals_sigma", "Residuals Sigma Phi");
	357	residuals_sigma->SetLineColor(kRed);
	358	residuals_sigma->SetMarkerStyle(23);
359	residuals_sigma->SetMarkerColor(kRed);
360	for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
361	Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
362	Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
363	h = h2->ProjectionY("py", ibin, ibin);
364	h->Fit(&f, "QN", "", -0.5, 0.5);
365	residuals_mean->SetPoint(ibin - 1, phi, f.GetParameter(1));
366	residuals_mean->SetPointError(ibin - 1, dphi, f.GetParError(1));
367	residuals_sigma->SetPoint(ibin - 1, phi, f.GetParameter(2));
368	residuals_sigma->SetPointError(ibin - 1, dphi, f.GetParError(2));
7102d1b1	369	}
d2381af5	370	fHistos->Add(residuals_mean);
d2381af5	371	fHistos->Add(residuals_sigma);
77203477	372	}
	373
	374	//________________________________________________________
	375	void AliTRDtrackingResolution::SetDebugLevel(Int_t level){
	376	fDebugLevel = level;
	377	if(!fDebugLevel) return;
	378	if(fDebugStream) return;
	379	fDebugStream = new TTreeSRedirector("TRD.Resolution.root");
	380	}
aaf47b30	381
	382	//________________________________________________________
	383	void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
	384	{
	385	fReconstructor->SetRecoParam(r);
	386	}