[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)
{
  AliCDBManager *cdb = AliCDBManager::Instance();
  cdb->SetDefaultStorage("local://$ALICE_ROOT");
  cdb->SetRun(0);
  TGeoManager::Import("geometry.root");
  fReconstructor = new AliTRDReconstructor();
  fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  AliTRDtrackerV1::SetNTimeBins(24);

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

//________________________________________________________
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){ 
   	     Float_t phiMC;
   	     if(Resolution(fTracklet, fInfo, phiMC)) phi = phiMC;
   	   	}
   	  }

      // Do clusters residuals
      fTracklet->Fit(kFALSE);
      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, Float_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->SetYref(0, zmc);
  tracklet->SetYref(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
      << "\n";
  }

  return kTRUE;
}


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

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

 	 // y resolution
	  TF1 *f = new TF1("f1", "gaus", -.5, .5);  
	  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, "q", "goff", -.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, "q", "goff", -.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);
	
	
	  delete f;
  }
}

//________________________________________________________
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	/**************************************************************************
	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 *
39779ce6	8	* documentation strictly for non-commercialf purposes is hereby granted *
77203477	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	**************************************************************************/
	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	AliCDBManager *cdb = AliCDBManager::Instance();
	69	cdb->SetDefaultStorage("local://$ALICE_ROOT");
	70	cdb->SetRun(0);
	71	TGeoManager::Import("geometry.root");
	72	fReconstructor = new AliTRDReconstructor();
	73	fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
	74	AliTRDtrackerV1::SetNTimeBins(24);
	75
77203477	76	DefineInput(0, TObjArray::Class());
	77	DefineOutput(0, TList::Class());
	78	}
	79
	80	//________________________________________________________
	81	void AliTRDtrackingResolution::ConnectInputData(Option_t *){
39779ce6	82	fTracks = dynamic_cast<TObjArray *>(GetInputData(0));
77203477	83	}
	84
	85	//________________________________________________________
	86	void AliTRDtrackingResolution::CreateOutputObjects()
	87	{
	88	// spatial resolution
77203477	89	OpenFile(0, "RECREATE");
39779ce6	90
	91	fHistos = new TList();
	92
	93	// Resolution histos
4b8f8a35	94	if(HasMCdata()){
	95	// tracklet resolution [0]
	96	fHistos->AddAt(new TH2I("fY", "", 21, -21., 21., 100, -.5, .5), 0);
	97	// tracklet angular resolution [1]
	98	fHistos->AddAt(new TH2I("fPhi", "", 21, -21., 21., 100, -10., 10.), 1);
	99	}
39779ce6	100	// Residual histos
39779ce6	101	// cluster to tracklet residuals [2]
4b8f8a35	102	Int_t position = HasMCdata() ? 2 : 0;
4b8f8a35	103	fHistos->AddAt(new TH2I("fYClRes", "", 21, -21., 21., 100, -.5, .5), position);
77203477	104	}
	105
	106	//________________________________________________________
7102d1b1	107	void AliTRDtrackingResolution::Exec(Option_t *)
7102d1b1	108	{
77203477	109	// spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
77203477	110	// angular Resolution: res = Tracklet angle - TrackRef Angle
7102d1b1	111
39779ce6	112	Int_t nTrackInfos = fTracks->GetEntriesFast();
39779ce6	113	if(fDebugLevel>=2) printf("Number of Histograms: %d\n", fHistos->GetEntries());
aaf47b30	114
aaf47b30	115	Double_t dy, dz;
39779ce6	116	AliTrackPoint tr[kNLayers], tk[kNLayers];
aaf47b30	117	AliTRDtrackV1 *fTrack = 0x0;
77203477	118	AliTRDtrackInfo *fInfo = 0x0;
	119	if(fDebugLevel>=2) printf("Number of TrackInfos: %d\n", nTrackInfos);
	120	for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
77203477	121	// check if ESD and MC-Information are available
39779ce6	122	if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
aaf47b30	123	if(!(fTrack = fInfo->GetTRDtrack())) continue;
7102d1b1	124	if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
7102d1b1	125
39779ce6	126
39779ce6	127	Int_t npts = 0;
aaf47b30	128	AliTRDseedV1 *fTracklet = 0x0;
77203477	129	for(Int_t iplane = 0; iplane < kNLayers; iplane++){
aaf47b30	130	if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
	131	if(!fTracklet->IsOK()) continue;
	132
39779ce6	133	// Book point arrays
	134	tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
	135	tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
	136	npts++;
	137
7102d1b1	138	if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
7102d1b1	139
39779ce6	140	Float_t phi = TMath::ATan(fTracklet->GetYref(1));
7102d1b1	141
39779ce6	142	// RESOLUTION (compare to MC)
4b8f8a35	143	if(HasMCdata()){
	144	if(fInfo->GetNTrackRefs() >= 2){
	145	Float_t phiMC;
	146	if(Resolution(fTracklet, fInfo, phiMC)) phi = phiMC;
	147	}
	148	}
aaf47b30	149
39779ce6	150	// Do clusters residuals
d2b9977a	151	fTracklet->Fit(kFALSE);
4b8f8a35	152	Int_t histpos = HasMCdata() ? 2 : 0;
aaf47b30	153	AliTRDcluster *c = 0x0;
	154	for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
	155	if(!(c = fTracklet->GetClusters(ic))) continue;
39779ce6	156
39779ce6	157	dy = fTracklet->GetYat(c->GetX()) - c->GetY();
4b8f8a35	158	((TH2I)fHistos->At(histpos))->Fill(phiTMath::RadToDeg(), dy);
39779ce6	159	if(fDebugLevel>=1){
	160	Float_t q = c->GetQ();
	161	(*fDebugStream) << "ClsTrkltResidual"
	162	<< "plane=" << iplane
	163	<< "phi=" << phi
	164	<< "q=" << q
	165	<< "dy=" << dy
	166	<< "\n";
	167	}
aaf47b30	168	}
aaf47b30	169	}
39779ce6	170
	171
	172	// this protection we might drop TODO
	173	if(fTrack->GetNumberOfTracklets() < 6) continue;
	174
	175	AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
	176	for(Int_t ip=0; ip<npts; ip++){
	177	dy = tk[ip].GetY() - tr[ip].GetY();
	178	dz = tk[ip].GetZ() - tr[ip].GetZ();
aaf47b30	179	if(fDebugLevel>=1){
	180	(*fDebugStream) << "ResidualsRT"
	181	<< "dy=" << dy
	182	<< "dz=" << dz
	183	/* << "phi=" << phi
	184	<< "theta=" << theta
	185	<< "dphi=" << dphi*/
	186	<< "\n";
	187	}
	188	}
	189
	190	// AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
	191	// for(Int_t ip=0; ip<nc; ip++){
	192	// dy = cl[ip].GetY() - tr[ip].GetY();
	193	// dz = cl[ip].GetZ() - tr[ip].GetZ();
	194	// if(fDebugLevel>=1){
	195	// (*fDebugStream) << "ResidualsKF"
	196	// << "dy=" << dy
	197	// << "dz=" << dz
	198	// /* << "phi=" << phi
	199	// << "theta=" << theta
	200	// << "dphi=" << dphi*/
	201	// << "\n";
	202	// }
	203	// }
39779ce6	204
39779ce6	205
77203477	206	}
39779ce6	207	PostData(0, fHistos);
77203477	208	}
77203477	209
39779ce6	210
	211	//________________________________________________________
	212	Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 tracklet, AliTRDtrackInfo fInfo, Float_t &phi)
	213	{
	214
	215	AliTrackReference fTrackRefs[2] = {0x0, 0x0}, tempTrackRef = 0x0;
	216
	217	// check for 2 track ref where the radial position has a distance less than 3.7mm
	218	Int_t nFound = 0;
	219	for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
	220	if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
	221	if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
	222	if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
	223	fTrackRefs[nFound++] = tempTrackRef;
	224	if(nFound == 2) break;
	225	}
	226	if(nFound < 2){
	227	if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
	228	return kFALSE;
	229	}
	230	// We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
	231
	232
	233	// RESOLUTION
	234	Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
	235	if(dx <= 0.){
	236	if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
	237	return kFALSE;
	238	}
	239	Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
	240	Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
	241	Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
	242	Double_t ymc = fTrackRefs[1]->LocalY() - dydx*dx0;
	243	Double_t zmc = fTrackRefs[1]->Z() - dzdx*dx0;
	244
	245	// recalculate tracklet based on the MC info
	246	tracklet->SetYref(0, zmc);
	247	tracklet->SetYref(1, dzdx);
	248	tracklet->Fit();
	249	Double_t dy = tracklet->GetYfit(0) - ymc;
	250	Double_t dz = tracklet->GetZfit(0) - zmc;
	251
	252	//res_y *= 100; // in mm
	253	Double_t momentum = fTrackRefs[0]->P();
	254
	255	phi = TMath::ATan(dydx);
	256	Double_t theta = TMath::ATan(dzdx);
	257	Double_t dphi = TMath::ATan(tracklet->GetYfit(1)) - phi;
	258	if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
	259
	260	// Fill Histograms
	261	if(TMath::Abs(dx-3.7)<1.E-3){
	262	((TH2I)fHistos->At(0))->Fill(phiTMath::RadToDeg(), dy);
	263	((TH2I)fHistos->At(1))->Fill(phiTMath::RadToDeg(), dphi*TMath::RadToDeg());
	264	}
	265	// Fill Debug Tree
	266	if(fDebugLevel>=1){
	267	Int_t iplane = tracklet->GetPlane();
	268	(*fDebugStream) << "TrkltResolution"
	269	<< "plane=" << iplane
	270	<< "p=" << momentum
	271	<< "dx=" << dx
	272	<< "dy=" << dy
	273	<< "dz=" << dz
274	<< "phi=" << phi
275	<< "theta=" << theta
276	<< "dphi=" << dphi
277	<< "\n";
278	}
279
280	return kTRUE;
281	}
282
283
77203477	284	//________________________________________________________
7102d1b1	285	void AliTRDtrackingResolution::Terminate(Option_t *)
7102d1b1	286	{
77203477	287	if(fDebugStream) delete fDebugStream;
7102d1b1	288
4b8f8a35	289	if(HasMCdata()){
	290	//process distributions
	291	fHistos = dynamic_cast<TList*>(GetOutputData(0));
	292	if (!fHistos) {
	293	Printf("ERROR: list not available");
	294	return;
	295	}
	296	TH2I *h2 = 0x0;
	297	TH1D *h = 0x0;
	298
	299	// y resolution
	300	TF1 *f = new TF1("f1", "gaus", -.5, .5);
	301	h2 = (TH2I*)fHistos->At(0);
	302	TGraphErrors *gm = new TGraphErrors(h2->GetNbinsX());
	303	gm->SetNameTitle("meany", "Mean dy");
	304	TGraphErrors *gs = new TGraphErrors(h2->GetNbinsX());
	305	gs->SetNameTitle("sigmy", "Sigma y");
	306	for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
	307	Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
	308	f->SetParameter(1, 0.);f->SetParameter(2, 2.e-2);
	309	h = h2->ProjectionY("py", iphi, iphi);
	310	h->Fit(f, "q", "goff", -.5, .5);
	311	Int_t jphi = iphi -1;
	312	gm->SetPoint(jphi, phi, f->GetParameter(1));
	313	gm->SetPointError(jphi, 0., f->GetParError(1));
	314	gs->SetPoint(jphi, phi, f->GetParameter(2));
	315	gs->SetPointError(jphi, 0., f->GetParError(2));
	316	}
	317	fHistos->Add(gm);
	318	fHistos->Add(gs);
	319
	320	// phi resolution
	321	h2 = (TH2I*)fHistos->At(1);
	322	gm = new TGraphErrors(h2->GetNbinsX());
	323	gm->SetNameTitle("meanphi", "Mean Phi");
	324	gs = new TGraphErrors(h2->GetNbinsX());
	325	gs->SetNameTitle("sigmphi", "Sigma Phi");
	326	for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
	327	Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
	328	f->SetParameter(1, 0.);f->SetParameter(2, 2.e-2);
	329	h = h2->ProjectionY("py", iphi, iphi);
	330	h->Fit(f, "q", "goff", -.5, .5);
	331	Int_t jphi = iphi -1;
	332	gm->SetPoint(jphi, phi, f->GetParameter(1));
	333	gm->SetPointError(jphi, 0., f->GetParError(1));
	334	gs->SetPoint(jphi, phi, f->GetParameter(2));
	335	gs->SetPointError(jphi, 0., f->GetParError(2));
	336	}
	337	fHistos->Add(gm);
	338	fHistos->Add(gs);
	339
	340
	341	delete f;
7102d1b1	342	}
77203477	343	}
	344
	345	//________________________________________________________
	346	void AliTRDtrackingResolution::SetDebugLevel(Int_t level){
	347	fDebugLevel = level;
	348	if(!fDebugLevel) return;
	349	if(fDebugStream) return;
	350	fDebugStream = new TTreeSRedirector("TRD.Resolution.root");
	351	}
aaf47b30	352
	353	//________________________________________________________
	354	void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
	355	{
	356	fReconstructor->SetRecoParam(r);
	357	}