[u/mrichter/AliRoot.git] / PWGPP / AliTrackComparison.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-commercial 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.                  *
 **************************************************************************/

//-------------------------------------------------------------------
// Base class to test the extrapolation performance from TPC to outer 
// detectors. Several member functions AddTracks() with different
// arguments can be called to execute extrapolation and several 
// THnSparse are filled with residuls and basic track information
//
// Anthor: R.Ma, M.Ivanov 02/04/2011
//--------------------------------------------------------------------
/* $Id:$ */

#include "AliTrackComparison.h"
#include "AliExternalTrackParam.h"
#include "AliTrackReference.h"
#include "AliTracker.h"
#include "AliTrackPointArray.h"
#include "THnSparse.h"
#include "TParticle.h"
#include "TMatrix.h"
#include "TParticlePDG.h"
#include "TParticle.h"
#include "TTreeStream.h"
#include "TAxis.h"
#include "TH1D.h"
#include "TTreeStream.h"

ClassImp(AliTrackComparison)

//________________________________________________________________________
AliTrackComparison::AliTrackComparison() 
  :TNamed()
  , fStep(1)
  , fLowBinDY(-10)
  , fUpBinDY(10)
  , fLowBinDZ(-10)
  , fUpBinDZ(10)
  , fLowBinDSnp(-0.5)
  , fUpBinDSnp(0.5)
  , fLowBinDTheta(-0.5)
  , fUpBinDTheta(0.5)
  , fLowBin1Pt(-3)
  , fUpBin1Pt(3)
  , fLowBin1PtLoss(-2)
  , fUpBin1PtLoss(2)
  , fNBinsDY(50)
  , fNBinsDZ(50)
  , fNBinsDSnp(50)
  , fNBinsDTheta(50)
  , fNBins1Pt(50)
  , fNBins1PtLoss(100)
  , fLayerID(-1)
  , fFillAll(kTRUE)
  , fNCombineBin(1)
{  
  //
  // Default constructor
  //
  for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
}


//________________________________________________________________________
AliTrackComparison::AliTrackComparison(const Text_t *name, const Text_t *title)
  :TNamed(name,title)
  , fStep(1)
  , fLowBinDY(-10)
  , fUpBinDY(10)
  , fLowBinDZ(-10)
  , fUpBinDZ(10)
  , fLowBinDSnp(-0.5)
  , fUpBinDSnp(0.5)
  , fLowBinDTheta(-0.5)
  , fUpBinDTheta(0.5)
  , fLowBin1Pt(-3)
  , fUpBin1Pt(3)
  , fLowBin1PtLoss(-2)
  , fUpBin1PtLoss(2)
  , fNBinsDY(50)
  , fNBinsDZ(50)
  , fNBinsDSnp(50)
  , fNBinsDTheta(50)
  , fNBins1Pt(50)
  , fNBins1PtLoss(100)
  , fLayerID(-1)
  , fFillAll(kTRUE)
  , fNCombineBin(1)
{
  //
  // Non default cosntructor
  //
  for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
}

//________________________________________________________________________
AliTrackComparison::AliTrackComparison(const AliTrackComparison& comp)
  :TNamed(comp)
  , fStep(comp.fStep)
  , fLowBinDY(comp.fLowBinDY)
  , fUpBinDY(comp.fUpBinDY)
  , fLowBinDZ(comp.fLowBinDZ)
  , fUpBinDZ(comp.fUpBinDZ)
  , fLowBinDSnp(comp.fLowBinDSnp)
  , fUpBinDSnp(comp.fUpBinDSnp)
  , fLowBinDTheta(comp.fLowBinDTheta)
  , fUpBinDTheta(comp.fUpBinDTheta)
  , fLowBin1Pt(comp.fLowBin1Pt)
  , fUpBin1Pt(comp.fUpBin1Pt)
  , fLowBin1PtLoss(comp.fLowBin1PtLoss)
  , fUpBin1PtLoss(comp.fUpBin1PtLoss)
  , fNBinsDY(comp.fNBinsDY)
  , fNBinsDZ(comp.fNBinsDZ)
  , fNBinsDSnp(comp.fNBinsDSnp)
  , fNBinsDTheta(comp.fNBinsDTheta)
  , fNBins1Pt(comp.fNBins1Pt)
  , fNBins1PtLoss(comp.fNBins1PtLoss)
  , fLayerID(comp.fLayerID)
  , fFillAll(comp.fFillAll)
  , fNCombineBin(comp.fNCombineBin)
{
  //
  // copy constructor
  //

  for (Int_t i=0;i<6; i++) fResolHisto[i]=comp.fResolHisto[i];
}

//________________________________________________________________________
AliTrackComparison& AliTrackComparison::operator=(const AliTrackComparison& comp)
{
  //
  //
  //
  if(this != &comp) {
    TNamed::operator=(comp);
  }
  return *this;
}

//________________________________________________________________________
void AliTrackComparison::Init(){
  //
  //Initilized the output histograms
  //
  MakeHistos();

}

//________________________________________________________________________
AliTrackComparison::~AliTrackComparison(){
  //
  //
  //
}
 
//________________________________________________________________________
void AliTrackComparison::Analyze() {
  //
  //
  //
}

//________________________________________________________________________
Int_t AliTrackComparison::AddTracks(AliTrackReference *ref0,  AliTrackReference *ref1, Double_t mass, Int_t charge)
{
  // Make track param out of track reference
  // Test propagation from ref0 to ref1
  // Fill the THnSparse with ref0 and ref1

  AliExternalTrackParam *param0 = 0;
  AliExternalTrackParam *param1 = 0;
   
  param0=MakeTrack(ref0,charge);
  param1=MakeTrack(ref1,charge);

  if (!param0 || !param1) return 0;
  
  Double_t tr1Pt = param0->GetSigned1Pt();

  if(!PropagateToPoint(param0,param1,mass)) return 0;

  FillHistos(param0,param1,tr1Pt);
  return 1;
}

//________________________________________________________________________
Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0,  AliExternalTrackParam *param1, Double_t mass)
{
  //Test propagation from param0 to param1
  //Fill the THnSparse with param0 and param1
  //
  Double_t tr1Pt=param0->GetSigned1Pt();
  if( !PropagateToPoint(param0,param1,mass)) return 0;
  FillHistos(param0,param1,tr1Pt);
  return 1;
}

//________________________________________________________________________
Int_t AliTrackComparison::AddTracks(AliExternalTrackParam *param0,  AliTrackPoint *point1, Double_t mass, Double_t energy, Double_t *vxyz)
{
  //Test propagation from param0 to point1
  //This function is usually called in real data analysis when only the
  //position of the track point is known. In this case, the angles 
  //in the track param are not the angles of the track momentum, but position.
  //Only the first two and last two THnSparse are meaninglful and should be 
  //filled. Set this via SetFillAll(kFALSE)

  Double_t gPos[3]= {point1->GetX(),point1->GetY(),point1->GetZ()};

  Double_t pos[3], pxyz[3];
  for(Int_t i=0; i<3; i++) 
    pos[i] = gPos[i]-vxyz[i];
  Double_t R = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
  for(Int_t i=0; i<3; i++) pxyz[i]= energy*pos[i]/R;

  Double_t cv[21];
  for (Int_t i=0; i<21;i++) cv[i]=0;
  AliExternalTrackParam * param1 = new AliExternalTrackParam(gPos,pxyz,cv,param0->Charge());

  if(!param1) return 0;
  Double_t tr1Pt = param0->GetSigned1Pt();
  if(!PropagateToPoint(param0,param1,mass)) return 0;

  FillHistos(param0,param1,tr1Pt);
  return 1;
}

//________________________________________________________________________
Bool_t AliTrackComparison::PropagateToPoint(AliExternalTrackParam *param0, AliExternalTrackParam *param1,  Double_t mass)
{
  //
  //Extrapolate is performed 
  //
  Double_t radius = param1->GetX();
  param0->Rotate(param1->GetAlpha());
  AliTracker::PropagateTrackToBxByBz(param0, radius+fStep, mass, fStep, kFALSE,0.99,-1);
  Bool_t isOK = param0->PropagateTo(radius,AliTracker::GetBz());
  return isOK;
}

//________________________________________________________________________
AliExternalTrackParam * AliTrackComparison::MakeTrack(const AliTrackReference* ref, Int_t charge)
{
  //
  // Make track out of the track ref
  // the covariance matrix - equal 0 - starting from ideal MC position
  Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
  Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
  Double_t cv[21];
  for (Int_t i=0; i<21;i++) cv[i]=0;
  AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
  return param;
}

//________________________________________________________________________
Long64_t AliTrackComparison::Merge(TCollection *const li) {
  //
  //Merge the comparison objects
  //   
  TIterator* iter = li->MakeIterator();
  AliTrackComparison* comp = 0;
  TString strName(GetName());
  while ((comp = (AliTrackComparison*)iter->Next())) {
    if (!comp->InheritsFrom(AliTrackComparison::Class())) {
      return -1;
    }
    if (strName.CompareTo(comp->GetName())!=0) return -1;
    // add histograms here...
    Add(comp);
  }
  return 0;  
}

//________________________________________________________________________
void  AliTrackComparison::Add(AliTrackComparison *const comp){
  //
  // Add THnSparse
  //
  for (Int_t i=0;i<6;i++){
    if (!fResolHisto[i]) 
      continue;
    THnSparse * h0 = (THnSparse*)fResolHisto[i];
    THnSparse * h1 = (THnSparse*)comp->GetHnSparse(i);
    if (h0&&h1) h0->Add(h1);
  }
}


//________________________________________________________________________
void AliTrackComparison::MakeHistos()
{
  //
  //Called in Init() to initialize histograms
  // 
  Double_t xminTrack[5], xmaxTrack[5];
  Int_t   binsTrack[5];
  TString axisName[5];
  TString axisTitle[5];
  TString hisNames[6]={"DeltaY","DeltaZ","DeltaSnp","DeltaTheta","Delta1Pt","1PtLoss"};
  Double_t lowBins[6]={fLowBinDY, fLowBinDZ, fLowBinDSnp, fLowBinDTheta, fLowBin1Pt, fLowBin1PtLoss};
  Double_t upBins[6]={fUpBinDY, fUpBinDZ, fUpBinDSnp, fUpBinDTheta, fUpBin1Pt, fUpBin1PtLoss};
  Int_t nBins[6] = {fNBinsDY, fNBinsDZ, fNBinsDSnp, fNBinsDTheta, fNBins1Pt, fNBins1PtLoss};
  //
  axisName[0]   ="#Delta";
  axisTitle[0]  ="#Delta";
  //
  binsTrack[1] =40;
  xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
  axisName[1]  ="tanTheta";
  axisTitle[1]  ="tan(#Theta)";
  //
  binsTrack[2] =180;
  xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi(); 
  axisName[2]  ="phi";
  axisTitle[2]  ="#phi";
  //
  binsTrack[3] =120;
  xminTrack[3] =-3.; xmaxTrack[3]=3.;   // 0.33 GeV cut 
  axisName[3]  ="1pt";
  axisTitle[3]  ="1/pt";
  //
  binsTrack[4] =22;
  xminTrack[4] =0; xmaxTrack[4]=22;  
  axisName[4]  ="LayerID";
  axisTitle[4] ="LayerID";

  for (Int_t ihis=0; ihis<6; ihis++){
    // modify ranges for bin 0
    //
    binsTrack[0]=nBins[ihis];
    xminTrack[0]=lowBins[ihis];
    xmaxTrack[0]=upBins[ihis];
    fResolHisto[ihis] = new THnSparseF(hisNames[ihis],hisNames[ihis],  5, binsTrack,xminTrack, xmaxTrack);
    for (Int_t ivar=0;ivar<5;ivar++){
      fResolHisto[ihis]->GetAxis(ivar)->SetName(axisName[ivar].Data());
      fResolHisto[ihis]->GetAxis(ivar)->SetTitle(axisTitle[ivar].Data());
    }
  }
}

//________________________________________________________________________
void AliTrackComparison::FillHistos(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t tr1Pt)
{
  //Fill the THnSparse. 
  //In case of not filling all, only the 4 out of 6 THnSparse are filed
  //
  Double_t dY     = param1->GetY()-param0->GetY();        //Residual in Y
  Double_t dZ     = param1->GetZ()-param0->GetZ();        //Residual in Z
  Double_t dSnp   = param1->GetSnp()-param0->GetSnp();    //Residual in sin(phi)
  Double_t dTheta = param1->Theta()-param0->Theta();      //Residual in Theta
  Double_t d1Pt   = param1->GetSigned1Pt()-param0->GetSigned1Pt(); //Residual in 1/pT
  Double_t d1PtLoss = param0->GetSigned1Pt()-tr1Pt;       //Corrected energy loss

  Double_t globalPos[3];
  param1->GetXYZ(globalPos);
  //printf("param1: Atan(y,x)=%5.3f | phi=%5.3f\n",TMath::ATan2(globalPos[1],globalPos[0]),param1->Phi());
  //printf("trkP=%5.3f | param0=%5.3f | param1=%5.3f\n",trP,param0->GetP(),param1->GetP());

  Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
  Double_t argu[6][5];
  for(Int_t j=0; j<6; j++)
    {
      if(!fFillAll && j<4 && j>1) continue;
      argu[j][0] = residual[j];
      argu[j][1] = param1->GetTgl();
      argu[j][2] = TMath::ATan2(globalPos[1],globalPos[0]);
      argu[j][3] = param1->GetSigned1Pt();
      argu[j][4] = fLayerID;
      fResolHisto[j]->Fill(argu[j]);
    } 
}

//________________________________________________________________________
void   AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
  //
  // make a distortion map out of the residual histogram
  // Results are written to the debug streamer - pcstream
  // Parameters:
  //   pcstream   - file to write the tree
  //   run        - run number
  //   refX       - track matching reference X
  //   type       - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
  // THnSparse axes:
  // OBJ: TAxis     #Delta  #Delta
  // OBJ: TAxis     tanTheta        tan(#Theta)
  // OBJ: TAxis     phi     #phi
  // OBJ: TAxis     1/pt    1/pt
  // OBJ: TAxis     LayerID LayerID  
  // marian.ivanov@cern.ch
  
  //Double_t refX=365;   //temporary
  Int_t    run=0;   //temporary
  TString hname = Form("%s_%d",GetName(),type);
  THnSparse * his0= GetHnSparse(type);
  
  TString dsName;
  dsName=GetName();
  dsName+=Form("Debug_%d.root",type);
  printf(" Create debug streamer \n");
  dsName.ReplaceAll(" ","");
  TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());


  const Int_t kMinEntries=10;
  Double_t bz=AliTrackerBase::GetBz();
  Int_t idim[4]={0,1,2,3};
  //
  //
  //
  Int_t nbins3=his0->GetAxis(3)->GetNbins();
  Int_t first3=his0->GetAxis(3)->GetFirst();
  Int_t last3 =his0->GetAxis(3)->GetLast();
  //
  for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){   // axis 3 - 1/pt
    his0->GetAxis(3)->SetRange(TMath::Max(ibin3-fNCombineBin,1),TMath::Min(ibin3+fNCombineBin,nbins3));
    Double_t      x3= his0->GetAxis(3)->GetBinCenter(ibin3);
    THnSparse * his3= his0->Projection(3,idim);         //projected histogram according selection 3
    //
    Int_t nbins2    = his3->GetAxis(2)->GetNbins();
    Int_t first2    = his3->GetAxis(2)->GetFirst();
    Int_t last2     = his3->GetAxis(2)->GetLast();
    //
    for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){   // axis 2 - phi
      his3->GetAxis(2)->SetRange(TMath::Max(ibin2-fNCombineBin,1),TMath::Min(ibin2+fNCombineBin,nbins2));
      Double_t      x2= his3->GetAxis(2)->GetBinCenter(ibin2);
      THnSparse * his2= his3->Projection(2,idim);         //projected histogram according selection 2
      Int_t nbins1     = his2->GetAxis(1)->GetNbins();
      Int_t first1     = his2->GetAxis(1)->GetFirst();
      Int_t last1      = his2->GetAxis(1)->GetLast();
      for (Int_t ibin1=first1; ibin1<last1; ibin1++){   //axis 1 - tan(theta)
	//
	Double_t       x1= his2->GetAxis(1)->GetBinCenter(ibin1);
	Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);

	his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1+fNCombineBin,nbins1));
	if (TMath::Abs(x1)<(fNCombineBin+0.5)*binWidth){
	  if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1,nbins1));
	  if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+fNCombineBin,nbins1));
	}
	TH1 * hisDelta = his2->Projection(0);
	//
	Double_t entries = hisDelta->GetEntries();
	Double_t mean=0, rms=0;
	if (entries>kMinEntries){
	  mean    = hisDelta->GetMean(); 
	  rms = hisDelta->GetRMS(); 
	}
	Double_t sector = 9.*x2/TMath::Pi();
	if (sector<0) sector+=18;
	Double_t dsec = sector-Int_t(sector)-0.5;
	Double_t quantiles[5]={0};
	Double_t mean75=0, rms75=0;
	Double_t prob[5]={0, 0.25,0.5,0.75,1};
	if (entries>kMinEntries){
	  hisDelta->GetQuantiles(5,quantiles,prob);
	  if(x3>0)hisDelta->GetXaxis()->SetRangeUser(quantiles[0], quantiles[3]);
	  else hisDelta->GetXaxis()->SetRangeUser(quantiles[1], quantiles[4]);
	  rms75=hisDelta->GetRMS();
	  mean75=hisDelta->GetMean();
	}
	
	Double_t pt = TMath::Abs(x3);
	Double_t z=refX*x1;
	(*pcstream)<<"distortion"<<
	  "run="<<run<<
	  "bz="<<bz<<
	  "theta="<<x1<<          // tan(theta)
	  "phi="<<x2<<            // phi
	  "z="<<z<<               // dummy z
	  "mpt="<<x3<<            // signed 1/pt
	  "1Pt="<<pt<<
	  "entries="<<entries<<   // entries
	  "mean="<<mean<<         // normal mean
	  //
	  "q25="<<quantiles[1]<<  // quantiles of distribution - importnat for asymetric distibutions
	  "q50="<<quantiles[2]<<  // quantiles of distribution - importnat for asymetric distibutions
	  "q75="<<quantiles[3]<<  // quantiles of distribution - importnat for asymetric distibutions
	  "mean75="<<mean75<<     // mean of the truncated distribution 0-75%
	  "rms75="<<rms75<<       // rms of the truncated distibution up to 0-75%
	  //
	  "rms="<<rms<<           // normal rms
	  "refX="<<refX<<         // track matching refernce plane
	  "type="<<type<<         // tpye of residuals
	  "sector="<<sector<<     // sector number according to TPC standard
	  "dsec="<<dsec<<
	  "\n";
	delete hisDelta;
	//printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
      }
      delete his2;
    }
    delete his3;
  }

  printf(" Finished! Close debug streamer \n");
  delete pcstream;
}

/*
  .x ~/rootlogon.C
   //TFile f("outFile.root");
   //TList * list = f.Get("jthaeder_OutterTracking");
  TFile f("/u/alice-st/marr/EnergyCorrection/TrainAnalysis/trunk/marr_TrackingTest.root");
  TList * list = f.Get("marr_TrackingTest");

  AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
  AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
  AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
  AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
  Double_t refX=438;

  compEMCAL-> MakeDistortionMap(refX,0);
  compEMCAL-> MakeDistortionMap(refX,1);
  compEMCAL-> MakeDistortionMap(refX,4);
  compEMCALEl-> MakeDistortionMap(refX,0);
  compEMCALEl-> MakeDistortionMap(refX,1);
  compEMCALEl-> MakeDistortionMap(refX,4);
  
  compTOF->SetNCombineBin(3)
  compTOF->MakeDistortionMap(365,4);
  compTOF->MakeDistortionMap(365,0);
  
  TFile f("emcalDistortion.root");
  

  // ideal case - no mislaignment dy only due energy loss correction
  TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
  //
  TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
  // delta 1/pt
  TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
  // (delta pt)/pt
  TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10");  // pions  - bias +-1 %
  //
  TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");  // electrons  - bias +-20 %
  //


*/
Commit	Line	Data
bdf94dc3	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 *
	8	* documentation strictly for non-commercial 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	**************************************************************************/
	15
	16	//-------------------------------------------------------------------
	17	// Base class to test the extrapolation performance from TPC to outer
	18	// detectors. Several member functions AddTracks() with different
	19	// arguments can be called to execute extrapolation and several
	20	// THnSparse are filled with residuls and basic track information
	21	//
	22	// Anthor: R.Ma, M.Ivanov 02/04/2011
	23	//--------------------------------------------------------------------
	24	/* $Id:$ */
	25
	26	#include "AliTrackComparison.h"
	27	#include "AliExternalTrackParam.h"
	28	#include "AliTrackReference.h"
	29	#include "AliTracker.h"
	30	#include "AliTrackPointArray.h"
	31	#include "THnSparse.h"
	32	#include "TParticle.h"
	33	#include "TMatrix.h"
	34	#include "TParticlePDG.h"
	35	#include "TParticle.h"
	36	#include "TTreeStream.h"
	37	#include "TAxis.h"
	38	#include "TH1D.h"
	39	#include "TTreeStream.h"
	40
	41	ClassImp(AliTrackComparison)
	42
	43	//________________________________________________________________________
	44	AliTrackComparison::AliTrackComparison()
	45	:TNamed()
	46	, fStep(1)
	47	, fLowBinDY(-10)
	48	, fUpBinDY(10)
	49	, fLowBinDZ(-10)
	50	, fUpBinDZ(10)
	51	, fLowBinDSnp(-0.5)
	52	, fUpBinDSnp(0.5)
	53	, fLowBinDTheta(-0.5)
	54	, fUpBinDTheta(0.5)
	55	, fLowBin1Pt(-3)
	56	, fUpBin1Pt(3)
	57	, fLowBin1PtLoss(-2)
	58	, fUpBin1PtLoss(2)
	59	, fNBinsDY(50)
	60	, fNBinsDZ(50)
	61	, fNBinsDSnp(50)
	62	, fNBinsDTheta(50)
	63	, fNBins1Pt(50)
	64	, fNBins1PtLoss(100)
65	, fLayerID(-1)
66	, fFillAll(kTRUE)
67	, fNCombineBin(1)
68	{
69	//
70	// Default constructor
71	//
72	for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
73	}
74
75
76	//________________________________________________________________________
77	AliTrackComparison::AliTrackComparison(const Text_t name, const Text_t title)
78	:TNamed(name,title)
79	, fStep(1)
80	, fLowBinDY(-10)
81	, fUpBinDY(10)
82	, fLowBinDZ(-10)
83	, fUpBinDZ(10)
84	, fLowBinDSnp(-0.5)
85	, fUpBinDSnp(0.5)
86	, fLowBinDTheta(-0.5)
87	, fUpBinDTheta(0.5)
88	, fLowBin1Pt(-3)
89	, fUpBin1Pt(3)
90	, fLowBin1PtLoss(-2)
91	, fUpBin1PtLoss(2)
92	, fNBinsDY(50)
93	, fNBinsDZ(50)
94	, fNBinsDSnp(50)
95	, fNBinsDTheta(50)
96	, fNBins1Pt(50)
97	, fNBins1PtLoss(100)
98	, fLayerID(-1)
99	, fFillAll(kTRUE)
100	, fNCombineBin(1)
101	{
102	//
103	// Non default cosntructor
104	//
105	for (Int_t i=0;i<6; i++) fResolHisto[i]=0;
106	}
107
108	//________________________________________________________________________
109	AliTrackComparison::AliTrackComparison(const AliTrackComparison& comp)
110	:TNamed(comp)
111	, fStep(comp.fStep)
112	, fLowBinDY(comp.fLowBinDY)
113	, fUpBinDY(comp.fUpBinDY)
114	, fLowBinDZ(comp.fLowBinDZ)
115	, fUpBinDZ(comp.fUpBinDZ)
116	, fLowBinDSnp(comp.fLowBinDSnp)
117	, fUpBinDSnp(comp.fUpBinDSnp)
118	, fLowBinDTheta(comp.fLowBinDTheta)
119	, fUpBinDTheta(comp.fUpBinDTheta)
120	, fLowBin1Pt(comp.fLowBin1Pt)
121	, fUpBin1Pt(comp.fUpBin1Pt)
122	, fLowBin1PtLoss(comp.fLowBin1PtLoss)
123	, fUpBin1PtLoss(comp.fUpBin1PtLoss)
124	, fNBinsDY(comp.fNBinsDY)
125	, fNBinsDZ(comp.fNBinsDZ)
126	, fNBinsDSnp(comp.fNBinsDSnp)
127	, fNBinsDTheta(comp.fNBinsDTheta)
128	, fNBins1Pt(comp.fNBins1Pt)
129	, fNBins1PtLoss(comp.fNBins1PtLoss)
130	, fLayerID(comp.fLayerID)
131	, fFillAll(comp.fFillAll)
132	, fNCombineBin(comp.fNCombineBin)
133	{
134	//
135	// copy constructor
136	//
137
138	for (Int_t i=0;i<6; i++) fResolHisto[i]=comp.fResolHisto[i];
139	}
140
141	//________________________________________________________________________
142	AliTrackComparison& AliTrackComparison::operator=(const AliTrackComparison& comp)
143	{
144	//
145	//
146	//
147	if(this != &comp) {
148	TNamed::operator=(comp);
149	}
150	return *this;
151	}
152
153	//________________________________________________________________________
154	void AliTrackComparison::Init(){
155	//
156	//Initilized the output histograms
157	//
158	MakeHistos();
159
160	}
161
162	//________________________________________________________________________
163	AliTrackComparison::~AliTrackComparison(){
164	//
165	//
166	//
167	}
168
169	//________________________________________________________________________
170	void AliTrackComparison::Analyze() {
171	//
172	//
173	//
174	}
175
176	//________________________________________________________________________
177	Int_t AliTrackComparison::AddTracks(AliTrackReference ref0, AliTrackReference ref1, Double_t mass, Int_t charge)
178	{
179	// Make track param out of track reference
180	// Test propagation from ref0 to ref1
181	// Fill the THnSparse with ref0 and ref1
182
183	AliExternalTrackParam *param0 = 0;
184	AliExternalTrackParam *param1 = 0;
185
186	param0=MakeTrack(ref0,charge);
187	param1=MakeTrack(ref1,charge);
188
189	if (!param0 \|\| !param1) return 0;
190
191	Double_t tr1Pt = param0->GetSigned1Pt();
192
193	if(!PropagateToPoint(param0,param1,mass)) return 0;
194
195	FillHistos(param0,param1,tr1Pt);
196	return 1;
197	}
198
199	//________________________________________________________________________
200	Int_t AliTrackComparison::AddTracks(AliExternalTrackParam param0, AliExternalTrackParam param1, Double_t mass)
201	{
202	//Test propagation from param0 to param1
203	//Fill the THnSparse with param0 and param1
204	//
205	Double_t tr1Pt=param0->GetSigned1Pt();
206	if( !PropagateToPoint(param0,param1,mass)) return 0;
207	FillHistos(param0,param1,tr1Pt);
208	return 1;
209	}
210
211	//________________________________________________________________________
212	Int_t AliTrackComparison::AddTracks(AliExternalTrackParam param0, AliTrackPoint point1, Double_t mass, Double_t energy, Double_t *vxyz)
213	{
214	//Test propagation from param0 to point1
215	//This function is usually called in real data analysis when only the
216	//position of the track point is known. In this case, the angles
217	//in the track param are not the angles of the track momentum, but position.
218	//Only the first two and last two THnSparse are meaninglful and should be
219	//filled. Set this via SetFillAll(kFALSE)
220
221	Double_t gPos[3]= {point1->GetX(),point1->GetY(),point1->GetZ()};
222
223	Double_t pos[3], pxyz[3];
224	for(Int_t i=0; i<3; i++)
225	pos[i] = gPos[i]-vxyz[i];
f6139574	226	Double_t R = TMath::Sqrt(pos[0]pos[0]+pos[1]pos[1]+pos[2]*pos[2]);
bdf94dc3	227	for(Int_t i=0; i<3; i++) pxyz[i]= energy*pos[i]/R;
	228
	229	Double_t cv[21];
	230	for (Int_t i=0; i<21;i++) cv[i]=0;
	231	AliExternalTrackParam * param1 = new AliExternalTrackParam(gPos,pxyz,cv,param0->Charge());
	232
	233	if(!param1) return 0;
	234	Double_t tr1Pt = param0->GetSigned1Pt();
	235	if(!PropagateToPoint(param0,param1,mass)) return 0;
	236
	237	FillHistos(param0,param1,tr1Pt);
	238	return 1;
	239	}
	240
	241	//________________________________________________________________________
	242	Bool_t AliTrackComparison::PropagateToPoint(AliExternalTrackParam param0, AliExternalTrackParam param1, Double_t mass)
	243	{
	244	//
	245	//Extrapolate is performed
	246	//
	247	Double_t radius = param1->GetX();
	248	param0->Rotate(param1->GetAlpha());
	249	AliTracker::PropagateTrackToBxByBz(param0, radius+fStep, mass, fStep, kFALSE,0.99,-1);
	250	Bool_t isOK = param0->PropagateTo(radius,AliTracker::GetBz());
	251	return isOK;
	252	}
	253
	254	//________________________________________________________________________
	255	AliExternalTrackParam * AliTrackComparison::MakeTrack(const AliTrackReference* ref, Int_t charge)
	256	{
	257	//
	258	// Make track out of the track ref
	259	// the covariance matrix - equal 0 - starting from ideal MC position
	260	Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
	261	Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
	262	Double_t cv[21];
	263	for (Int_t i=0; i<21;i++) cv[i]=0;
	264	AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
	265	return param;
	266	}
	267
	268	//________________________________________________________________________
	269	Long64_t AliTrackComparison::Merge(TCollection *const li) {
	270	//
	271	//Merge the comparison objects
	272	//
	273	TIterator* iter = li->MakeIterator();
	274	AliTrackComparison* comp = 0;
	275	TString strName(GetName());
	276	while ((comp = (AliTrackComparison*)iter->Next())) {
	277	if (!comp->InheritsFrom(AliTrackComparison::Class())) {
	278	return -1;
	279	}
	280	if (strName.CompareTo(comp->GetName())!=0) return -1;
	281	// add histograms here...
	282	Add(comp);
	283	}
	284	return 0;
	285	}
	286
	287	//________________________________________________________________________
	288	void AliTrackComparison::Add(AliTrackComparison *const comp){
	289	//
	290	// Add THnSparse
291	//
bdf94dc3	292	for (Int_t i=0;i<6;i++){
f6139574	293	if (!fResolHisto[i])
f6139574	294	continue;
bdf94dc3	295	THnSparse * h0 = (THnSparse*)fResolHisto[i];
	296	THnSparse * h1 = (THnSparse*)comp->GetHnSparse(i);
	297	if (h0&&h1) h0->Add(h1);
	298	}
	299	}
	300
	301
	302	//________________________________________________________________________
	303	void AliTrackComparison::MakeHistos()
	304	{
	305	//
	306	//Called in Init() to initialize histograms
	307	//
	308	Double_t xminTrack[5], xmaxTrack[5];
	309	Int_t binsTrack[5];
	310	TString axisName[5];
	311	TString axisTitle[5];
	312	TString hisNames[6]={"DeltaY","DeltaZ","DeltaSnp","DeltaTheta","Delta1Pt","1PtLoss"};
	313	Double_t lowBins[6]={fLowBinDY, fLowBinDZ, fLowBinDSnp, fLowBinDTheta, fLowBin1Pt, fLowBin1PtLoss};
	314	Double_t upBins[6]={fUpBinDY, fUpBinDZ, fUpBinDSnp, fUpBinDTheta, fUpBin1Pt, fUpBin1PtLoss};
	315	Int_t nBins[6] = {fNBinsDY, fNBinsDZ, fNBinsDSnp, fNBinsDTheta, fNBins1Pt, fNBins1PtLoss};
	316	//
	317	axisName[0] ="#Delta";
	318	axisTitle[0] ="#Delta";
	319	//
	320	binsTrack[1] =40;
	321	xminTrack[1] =-1.0; xmaxTrack[1]=1.0;
	322	axisName[1] ="tanTheta";
	323	axisTitle[1] ="tan(#Theta)";
	324	//
	325	binsTrack[2] =180;
	326	xminTrack[2] =-TMath::Pi(); xmaxTrack[2]=TMath::Pi();
	327	axisName[2] ="phi";
	328	axisTitle[2] ="#phi";
	329	//
	330	binsTrack[3] =120;
	331	xminTrack[3] =-3.; xmaxTrack[3]=3.; // 0.33 GeV cut
	332	axisName[3] ="1pt";
	333	axisTitle[3] ="1/pt";
	334	//
	335	binsTrack[4] =22;
	336	xminTrack[4] =0; xmaxTrack[4]=22;
	337	axisName[4] ="LayerID";
	338	axisTitle[4] ="LayerID";
	339
	340	for (Int_t ihis=0; ihis<6; ihis++){
	341	// modify ranges for bin 0
	342	//
	343	binsTrack[0]=nBins[ihis];
	344	xminTrack[0]=lowBins[ihis];
	345	xmaxTrack[0]=upBins[ihis];
	346	fResolHisto[ihis] = new THnSparseF(hisNames[ihis],hisNames[ihis], 5, binsTrack,xminTrack, xmaxTrack);
	347	for (Int_t ivar=0;ivar<5;ivar++){
	348	fResolHisto[ihis]->GetAxis(ivar)->SetName(axisName[ivar].Data());
	349	fResolHisto[ihis]->GetAxis(ivar)->SetTitle(axisTitle[ivar].Data());
	350	}
	351	}
	352	}
	353
	354	//________________________________________________________________________
	355	void AliTrackComparison::FillHistos(AliExternalTrackParam param0, AliExternalTrackParam param1, Double_t tr1Pt)
	356	{
	357	//Fill the THnSparse.
	358	//In case of not filling all, only the 4 out of 6 THnSparse are filed
359	//
360	Double_t dY = param1->GetY()-param0->GetY(); //Residual in Y
361	Double_t dZ = param1->GetZ()-param0->GetZ(); //Residual in Z
362	Double_t dSnp = param1->GetSnp()-param0->GetSnp(); //Residual in sin(phi)
363	Double_t dTheta = param1->Theta()-param0->Theta(); //Residual in Theta
364	Double_t d1Pt = param1->GetSigned1Pt()-param0->GetSigned1Pt(); //Residual in 1/pT
365	Double_t d1PtLoss = param0->GetSigned1Pt()-tr1Pt; //Corrected energy loss
366
367	Double_t globalPos[3];
368	param1->GetXYZ(globalPos);
369	//printf("param1: Atan(y,x)=%5.3f \| phi=%5.3f\n",TMath::ATan2(globalPos[1],globalPos[0]),param1->Phi());
370	//printf("trkP=%5.3f \| param0=%5.3f \| param1=%5.3f\n",trP,param0->GetP(),param1->GetP());
371
372	Double_t residual[6] = {dY,dZ,dSnp,dTheta,d1Pt,d1PtLoss};
373	Double_t argu[6][5];
374	for(Int_t j=0; j<6; j++)
375	{
376	if(!fFillAll && j<4 && j>1) continue;
377	argu[j][0] = residual[j];
378	argu[j][1] = param1->GetTgl();
379	argu[j][2] = TMath::ATan2(globalPos[1],globalPos[0]);
380	argu[j][3] = param1->GetSigned1Pt();
381	argu[j][4] = fLayerID;
382	fResolHisto[j]->Fill(argu[j]);
383	}
384	}
385
386	//________________________________________________________________________
387	void AliTrackComparison::MakeDistortionMap(Double_t refX, Int_t type){
388	//
389	// make a distortion map out of the residual histogram
390	// Results are written to the debug streamer - pcstream
391	// Parameters:
392	// pcstream - file to write the tree
393	// run - run number
394	// refX - track matching reference X
395	// type - 0- y 1-z,2 -snp, 3-theta, 4-1/pt, 5-corrected 1/pT
396	// THnSparse axes:
397	// OBJ: TAxis #Delta #Delta
398	// OBJ: TAxis tanTheta tan(#Theta)
399	// OBJ: TAxis phi #phi
400	// OBJ: TAxis 1/pt 1/pt
401	// OBJ: TAxis LayerID LayerID
402	// marian.ivanov@cern.ch
403
404	//Double_t refX=365; //temporary
405	Int_t run=0; //temporary
406	TString hname = Form("%s_%d",GetName(),type);
407	THnSparse * his0= GetHnSparse(type);
408
409	TString dsName;
410	dsName=GetName();
411	dsName+=Form("Debug_%d.root",type);
412	printf(" Create debug streamer \n");
413	dsName.ReplaceAll(" ","");
414	TTreeSRedirector *pcstream = new TTreeSRedirector(dsName.Data());
415
416
417	const Int_t kMinEntries=10;
418	Double_t bz=AliTrackerBase::GetBz();
419	Int_t idim[4]={0,1,2,3};
420	//
421	//
422	//
423	Int_t nbins3=his0->GetAxis(3)->GetNbins();
424	Int_t first3=his0->GetAxis(3)->GetFirst();
425	Int_t last3 =his0->GetAxis(3)->GetLast();
426	//
427	for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - 1/pt
428	his0->GetAxis(3)->SetRange(TMath::Max(ibin3-fNCombineBin,1),TMath::Min(ibin3+fNCombineBin,nbins3));
429	Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3);
430	THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3
431	//
432	Int_t nbins2 = his3->GetAxis(2)->GetNbins();
433	Int_t first2 = his3->GetAxis(2)->GetFirst();
434	Int_t last2 = his3->GetAxis(2)->GetLast();
435	//
436	for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi
437	his3->GetAxis(2)->SetRange(TMath::Max(ibin2-fNCombineBin,1),TMath::Min(ibin2+fNCombineBin,nbins2));
438	Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2);
439	THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2
440	Int_t nbins1 = his2->GetAxis(1)->GetNbins();
441	Int_t first1 = his2->GetAxis(1)->GetFirst();
442	Int_t last1 = his2->GetAxis(1)->GetLast();
443	for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - tan(theta)
444	//
445	Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1);
446	Double_t binWidth= his2->GetAxis(1)->GetBinWidth(ibin1);
447
448	his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1+fNCombineBin,nbins1));
449	if (TMath::Abs(x1)<(fNCombineBin+0.5)*binWidth){
450	if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-fNCombineBin,1),TMath::Min(ibin1,nbins1));
451	if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+fNCombineBin,nbins1));
452	}
453	TH1 * hisDelta = his2->Projection(0);
454	//
455	Double_t entries = hisDelta->GetEntries();
456	Double_t mean=0, rms=0;
457	if (entries>kMinEntries){
458	mean = hisDelta->GetMean();
459	rms = hisDelta->GetRMS();
460	}
461	Double_t sector = 9.*x2/TMath::Pi();
462	if (sector<0) sector+=18;
463	Double_t dsec = sector-Int_t(sector)-0.5;
464	Double_t quantiles[5]={0};
465	Double_t mean75=0, rms75=0;
466	Double_t prob[5]={0, 0.25,0.5,0.75,1};
467	if (entries>kMinEntries){
468	hisDelta->GetQuantiles(5,quantiles,prob);
469	if(x3>0)hisDelta->GetXaxis()->SetRangeUser(quantiles[0], quantiles[3]);
470	else hisDelta->GetXaxis()->SetRangeUser(quantiles[1], quantiles[4]);
471	rms75=hisDelta->GetRMS();
472	mean75=hisDelta->GetMean();
473	}
474
475	Double_t pt = TMath::Abs(x3);
476	Double_t z=refX*x1;
477	(*pcstream)<<"distortion"<<
478	"run="<<run<<
479	"bz="<<bz<<
480	"theta="<<x1<< // tan(theta)
481	"phi="<<x2<< // phi
482	"z="<<z<< // dummy z
483	"mpt="<<x3<< // signed 1/pt
484	"1Pt="<<pt<<
485	"entries="<<entries<< // entries
486	"mean="<<mean<< // normal mean
487	//
488	"q25="<<quantiles[1]<< // quantiles of distribution - importnat for asymetric distibutions
489	"q50="<<quantiles[2]<< // quantiles of distribution - importnat for asymetric distibutions
490	"q75="<<quantiles[3]<< // quantiles of distribution - importnat for asymetric distibutions
491	"mean75="<<mean75<< // mean of the truncated distribution 0-75%
492	"rms75="<<rms75<< // rms of the truncated distibution up to 0-75%
493	//
494	"rms="<<rms<< // normal rms
495	"refX="<<refX<< // track matching refernce plane
496	"type="<<type<< // tpye of residuals
497	"sector="<<sector<< // sector number according to TPC standard
498	"dsec="<<dsec<<
499	"\n";
500	delete hisDelta;
501	//printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean);
502	}
503	delete his2;
504	}
505	delete his3;
506	}
507
508	printf(" Finished! Close debug streamer \n");
509	delete pcstream;
510	}
511
512	/*
513	.x ~/rootlogon.C
514	//TFile f("outFile.root");
515	//TList * list = f.Get("jthaeder_OutterTracking");
516	TFile f("/u/alice-st/marr/EnergyCorrection/TrainAnalysis/trunk/marr_TrackingTest.root");
517	TList * list = f.Get("marr_TrackingTest");
518
519	AliTrackComparison * compTOF = list->FindObject("TPCOutToTOFIn");
520	AliTrackComparison * compEMCAL = list->FindObject("TPCOutToEMCalInPion");
521	AliTrackComparison * compEMCALEl = list->FindObject("TPCOutToEMCalInElec");
522	AliTrackComparison * compHMPID = list->FindObject("TPCOutToHMPIDIn");
523	Double_t refX=438;
524
525	compEMCAL-> MakeDistortionMap(refX,0);
526	compEMCAL-> MakeDistortionMap(refX,1);
527	compEMCAL-> MakeDistortionMap(refX,4);
528	compEMCALEl-> MakeDistortionMap(refX,0);
529	compEMCALEl-> MakeDistortionMap(refX,1);
530	compEMCALEl-> MakeDistortionMap(refX,4);
531
532	compTOF->SetNCombineBin(3)
533	compTOF->MakeDistortionMap(365,4);
534	compTOF->MakeDistortionMap(365,0);
535
536	TFile f("emcalDistortion.root");
537
538
539	// ideal case - no mislaignment dy only due energy loss correction
540	TPCOutToEMCalInElec_0->Draw("mean:mpt","entries>10");
541	//
542	TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10");
543	// delta 1/pt
544	TPCOutToEMCalInPion_4->Draw("mean:mpt","entries>10");
545	// (delta pt)/pt
546	TPCOutToEMCalInPion_4->Draw("mean/abs(mpt):mpt","entries>10"); // pions - bias +-1 %
547	//
548	TPCOutToEMCalInElec_4->Draw("mean/abs(mpt):mpt","entries>10"); // electrons - bias +-20 %
549	//
550
551
552
553
554	*/