[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronHelper.cxx

/*************************************************************************
* Copyright(c) 1998-2009, 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.                  *
**************************************************************************/

//
// Dielectron helper functions wrapped in a namespace
// 
//
// Authors: 
//   Jens Wiechula <Jens.Wiechula@cern.ch> 
//   Frederick Kramer <Frederick.Kramer@cern.ch> 
//   Julian Book <Julian.Book@cern.ch>


#include <TError.h>
#include <TMath.h>
#include <TObjString.h>
#include <TObjArray.h>
#include <TVectorD.h>
#include <TF1.h>
#include <TRandom.h>
#include <TProfile.h>

#include <AliVEvent.h>
#include <AliVParticle.h>
#include <AliKFParticle.h>
#include <AliESDtrackCuts.h>
#include <AliESDEvent.h>
#include <AliMCEvent.h>
#include <AliAODEvent.h>
#include <AliAODTracklets.h>
#include <AliMultiplicity.h>
#include <AliStack.h>

#include "AliDielectronVarCuts.h"
#include "AliDielectronTrackCuts.h"
#include "AliDielectronVarManager.h"
#include "AliDielectronHelper.h"

//_____________________________________________________________________________
TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make logarithmic binning
  // the user has to delete the array afterwards!!!
  //
  
  //check limits
  if (xmin<1e-20 || xmax<1e-20){
    Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
    return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax);
  }
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t expMax=TMath::Log(last/first);
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
  }
  return binLim;
}

//_____________________________________________________________________________
TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make linear binning
  // the user has to delete the array afterwards!!!
  //
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t binWidth=(last-first)/nbinsX;
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first+binWidth*(Double_t)i;
  }
  return binLim;
}

//_____________________________________________________________________________
TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
{
  //
  // Make arbitrary binning, bins separated by a ','
  //
  TString limits(bins);
  if (limits.IsNull()){
    Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
    return 0x0;
  }
  
  TObjArray *arr=limits.Tokenize(",");
  Int_t nLimits=arr->GetEntries();
  if (nLimits<2){
    Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
    delete arr;
    return 0x0;
  }
  
  TVectorD *binLimits=new TVectorD(nLimits);
  for (Int_t iLim=0; iLim<nLimits; ++iLim){
    (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
  }
  
  delete arr;
  return binLimits;
}


//_____________________________________________________________________________
Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
  // determination of Nch
  if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;

  AliStack *stack = ((AliMCEvent*)ev)->Stack();

  if (!stack) return -1;

  Int_t nParticles = stack->GetNtrack();
  Int_t nCh = 0;

  // count..
  for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
    if (!stack->IsPhysicalPrimary(iMc)) continue;

    TParticle* particle = stack->Particle(iMc);
    if (!particle) continue;
    if (particle->GetPDG()->Charge() == 0) continue;

    Float_t eta = particle->Eta();
    if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
  }

  return nCh;
}


//_____________________________________________________________________________
Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){
  // Compute the collision multiplicity based on AOD or ESD tracklets
  // Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()

  if (!ev) return -1;

  Int_t nTracklets = 0;
  Int_t nAcc = 0;
  
  if (ev->IsA() == AliAODEvent::Class()) {
    AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
    if (!tracklets) return -1;
    nTracklets = tracklets->GetNumberOfTracklets();
    for (Int_t nn = 0; nn < nTracklets; nn++) {
      Double_t theta = tracklets->GetTheta(nn);
      Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
      if (TMath::Abs(eta) < etaRange) nAcc++;
    }
  } else if (ev->IsA() == AliESDEvent::Class()) {
    nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
    for (Int_t nn = 0; nn < nTracklets; nn++) {
      Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
      if (TMath::Abs(eta) < etaRange) nAcc++;
    }
  } else return -1;

  return nAcc;
}


//________________________________________________________________
Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
  //
  // Correct the number of accepted tracklets based on the period and event vertex
  //

  Int_t runNo = event->GetRunNumber();

  Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
  Double_t refMult = 0.0;   // reference multiplicity
  
  if(runNo>114930 && runNo<117223) period = 0;
  if(runNo>119158 && runNo<120830) period = 1;
  if(runNo>122373 && runNo<126438) period = 2;
  if(runNo>127711 && runNo<130841) period = 3;
  if(period<0 || period>3) return uncorrectedNacc;

  if(type<0 || type>8) return uncorrectedNacc;
  if(type == 0) refMult = 5.0;         // SPD tracklets in |eta|<0.5 
  if(type == 1) refMult = 9.5;         // SPD tracklets in |eta|<1.0
  if(type == 2) refMult = 13.0;        // SPD tracklets in |eta|<1.6
  if(type == 3) refMult = 6.0;         // ITSTPC+ in |eta|<0.5
  if(type == 4) refMult = 12.0;        // ITSTPC+ in |eta|<1.0
  if(type == 5) refMult = 16.0;        // ITSTPC+ in |eta|<1.6
  if(type == 6) refMult = 6.0;         // ITSSA+ in |eta|<0.5
  if(type == 7) refMult = 12.0;        // ITSSA+ in |eta|<1.0
  if(type == 8) refMult = 15.0;        // ITSSA+ in |eta|<1.6

  if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;

  TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
  if(!estimatorAvg) return uncorrectedNacc;

  Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));

  Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);

  Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));

  return correctedNacc;
}

//_____________________________________________________________________________
Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
  // put a robust Nacc definition here

  if (!ev) return -1;
  
  AliDielectronVarCuts *varCuts = new AliDielectronVarCuts("VarCuts","VarCuts");
  varCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
  varCuts->AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
  varCuts->AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
  varCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
  varCuts->AddCut(AliDielectronVarManager::kNclsTPC,     70.0, 160.0);
  varCuts->AddCut(AliDielectronVarManager::kKinkIndex0,  -0.5,   0.5);   //noKinks
    
  AliDielectronTrackCuts *trkCuts = new AliDielectronTrackCuts("TrkCuts","TrkCuts");
  trkCuts->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
  trkCuts->SetRequireITSRefit(kTRUE);
  trkCuts->SetRequireTPCRefit(kTRUE);

  Int_t nRecoTracks = ev->GetNumberOfTracks();
  Int_t nAcc = 0;
  
  for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
    AliVTrack *track        = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
    if (!track) continue;
    if (!trkCuts->IsSelected(track)) continue;
    if (!varCuts->IsSelected(track)) continue;
    nAcc++;
  }
  
  delete varCuts;
  delete trkCuts;

  return nAcc;
}

//_____________________________________________________________________________
Double_t AliDielectronHelper::GetITSTPCMatchEff(const AliVEvent *ev){
  // recalulate the its-tpc matching efficiecy

  if (!ev) return -1;

  AliDielectronVarCuts *varCutsTPC = new AliDielectronVarCuts("VarCutsTPC","VarCutsTPC");
  varCutsTPC->AddCut(AliDielectronVarManager::kImpactParXY, -1.0,   1.0);
  varCutsTPC->AddCut(AliDielectronVarManager::kImpactParZ,  -3.0,   3.0);
  varCutsTPC->AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
  varCutsTPC->AddCut(AliDielectronVarManager::kTPCchi2Cl,    0.0,   4.0);
  varCutsTPC->AddCut(AliDielectronVarManager::kNclsTPC,     50.0, 160.0);
  AliDielectronTrackCuts *trkCutsTPC = new AliDielectronTrackCuts("TrkCutsTPC","TrkCutsTPC");
  trkCutsTPC->SetRequireTPCRefit(kTRUE);

  AliDielectronVarCuts *varCutsITS = new AliDielectronVarCuts("VarCutsITS","VarCutsITS");
  varCutsITS->AddCut(AliDielectronVarManager::kEta,         -0.9,   0.9);
  AliDielectronTrackCuts *trkCutsITS = new AliDielectronTrackCuts("TrkCutsITS","TrkCutsITS");
  trkCutsITS->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
  trkCutsITS->SetRequireITSRefit(kTRUE);


  Int_t nRecoTracks = ev->GetNumberOfTracks();
  Double_t nTPC = 0, nITS = 0;

  for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
    AliVTrack *track        = static_cast<AliVTrack*>(ev->GetTrack(iTrack));
    if (!track) continue;

    if(!trkCutsITS->IsSelected(track)) continue;
    if(!varCutsITS->IsSelected(track)) continue;
    nITS+=1.;

    if(!trkCutsTPC->IsSelected(track)) continue;
    if(!varCutsTPC->IsSelected(track)) continue;
    nTPC+=1.;

  }

  delete varCutsITS;
  delete trkCutsITS;
  delete varCutsTPC;
  delete trkCutsTPC;

  printf(" tracks TPC %.3e ITS %.3e = %.5f \n",nTPC,nITS,(nITS>0. ? nTPC/nITS : -1));
  return (nITS>0. ? nTPC/nITS : -1);
}


//_____________________________________________________________________________
void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
  // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
  if (!kfParticle) return;
  Double_t dx = 0.;
  Double_t dy = 0.;
  Double_t dz = 0.;

  if (ev){
    dx = ev->GetPrimaryVertex()->GetX()-0.;
    dy = ev->GetPrimaryVertex()->GetY()-0.;
    dz = ev->GetPrimaryVertex()->GetZ()-0.;
  }
  
  kfParticle->X() = kfParticle->GetX() - dx;
  kfParticle->Y() = kfParticle->GetY() - dy;
  kfParticle->Z() = kfParticle->GetZ() - dz;
  
  
  // Rotate the kf particle
  Double_t c = cos(angle);
  Double_t s = sin(angle);
  
  Double_t mA[8][ 8];
  for( Int_t i=0; i<8; i++ ){
    for( Int_t j=0; j<8; j++){
      mA[i][j] = 0;
    }
  }
  for( int i=0; i<8; i++ ){
    mA[i][i] = 1;
  }
  mA[0][0] =  c;  mA[0][1] = s;
  mA[1][0] = -s;  mA[1][1] = c;
  mA[3][3] =  c;  mA[3][4] = s;
  mA[4][3] = -s;  mA[4][4] = c;
  
  Double_t mAC[8][8];
  Double_t mAp[8];
  
  for( Int_t i=0; i<8; i++ ){
    mAp[i] = 0;
    for( Int_t k=0; k<8; k++){
      mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
    }
  }
  
  for( Int_t i=0; i<8; i++){
    kfParticle->Parameter(i) = mAp[i];
  }
  
  for( Int_t i=0; i<8; i++ ){
    for( Int_t j=0; j<8; j++ ){
      mAC[i][j] = 0;
      for( Int_t k=0; k<8; k++ ){
        mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
      }
    }
  }
  
  for( Int_t i=0; i<8; i++ ){
    for( Int_t j=0; j<=i; j++ ){
      Double_t xx = 0;
      for( Int_t k=0; k<8; k++){
        xx+= mAC[i][k]*mA[j][k];
      }
      kfParticle->Covariance(i,j) = xx;
    }
  }
  
  kfParticle->X() = kfParticle->GetX() + dx;
  kfParticle->Y() = kfParticle->GetY() + dy;
  kfParticle->Z() = kfParticle->GetZ() + dz;
  
}

//_____________________________________________________________________________
Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
  // TODO: add AODs
  // counting number of mother particles generated in given eta range and 2 particle decay
  if (!ev || ev->IsA()!=AliMCEvent::Class()) return -1;
  
  AliStack *stack = ((AliMCEvent*)ev)->Stack();
  
  if (!stack) return -1;
  
  Int_t nParticles = stack->GetNtrack();
  Int_t nMothers   = 0;
  
  // count..
  for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
    
    TParticle* particle = stack->Particle(iMc);
    if (!particle) continue;
    if (particle->GetPdgCode() != pdgMother)               continue;
    if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;

    if (particle->GetNDaughters() != 2)                 continue;
    // 1st daugther
    if (particle->GetFirstDaughter()>=nParticles ||
	particle->GetFirstDaughter()<0             ) continue;
    
    TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
    if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter)     continue;
    if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
    
    // 2nd daughter
    if (particle->GetLastDaughter()>=nParticles ||
	particle->GetLastDaughter()<0             ) continue;

    TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
    if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter)     continue;
    if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
    
    // primary
    if (prim != -1) {
      if(particle->IsPrimary() != prim) continue;
    }
    nMothers++;
  }
  return nMothers;
}
Commit	Line	Data
ffbede40	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, 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	// Dielectron helper functions wrapped in a namespace
	18	//
	19	//
	20	// Authors:
	21	// Jens Wiechula <Jens.Wiechula@cern.ch>
ba15fdfb	22	// Frederick Kramer <Frederick.Kramer@cern.ch>
41a64bfb	23	// Julian Book <Julian.Book@cern.ch>
ffbede40	24
	25
	26
	27
	28	#include <TError.h>
	29	#include <TMath.h>
	30	#include <TObjString.h>
	31	#include <TObjArray.h>
	32	#include <TVectorD.h>
5720c765	33	#include <TF1.h>
	34	#include <TRandom.h>
	35	#include <TProfile.h>
ffbede40	36
1201a1a9	37	#include <AliVEvent.h>
ba15fdfb	38	#include <AliVParticle.h>
1201a1a9	39	#include <AliKFParticle.h>
ba15fdfb	40	#include <AliESDtrackCuts.h>
	41	#include <AliESDEvent.h>
	42	#include <AliMCEvent.h>
	43	#include <AliAODEvent.h>
	44	#include <AliAODTracklets.h>
	45	#include <AliMultiplicity.h>
	46	#include <AliStack.h>
1201a1a9	47
a823f01b	48	#include "AliDielectronVarCuts.h"
a823f01b	49	#include "AliDielectronTrackCuts.h"
5720c765	50	#include "AliDielectronVarManager.h"
ffbede40	51	#include "AliDielectronHelper.h"
	52
	53	//_____________________________________________________________________________
	54	TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
	55	{
	56	//
	57	// Make logarithmic binning
	58	// the user has to delete the array afterwards!!!
	59	//
	60
	61	//check limits
	62	if (xmin<1e-20 \|\| xmax<1e-20){
	63	Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
	64	return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax);
	65	}
	66	if (xmax<xmin){
	67	Double_t tmp=xmin;
	68	xmin=xmax;
	69	xmax=tmp;
	70	}
	71	TVectorD *binLim=new TVectorD(nbinsX+1);
	72	Double_t first=xmin;
	73	Double_t last=xmax;
	74	Double_t expMax=TMath::Log(last/first);
	75	for (Int_t i=0; i<nbinsX+1; ++i){
	76	(binLim)[i]=firstTMath::Exp(expMax/nbinsX*(Double_t)i);
	77	}
	78	return binLim;
	79	}
	80
	81	//_____________________________________________________________________________
	82	TVectorD* AliDielectronHelper::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
	83	{
	84	//
	85	// Make linear binning
	86	// the user has to delete the array afterwards!!!
	87	//
	88	if (xmax<xmin){
	89	Double_t tmp=xmin;
	90	xmin=xmax;
	91	xmax=tmp;
	92	}
	93	TVectorD *binLim=new TVectorD(nbinsX+1);
	94	Double_t first=xmin;
	95	Double_t last=xmax;
	96	Double_t binWidth=(last-first)/nbinsX;
	97	for (Int_t i=0; i<nbinsX+1; ++i){
	98	(binLim)[i]=first+binWidth(Double_t)i;
	99	}
	100	return binLim;
	101	}
	102
	103	//_____________________________________________________________________________
	104	TVectorD* AliDielectronHelper::MakeArbitraryBinning(const char* bins)
	105	{
	106	//
	107	// Make arbitrary binning, bins separated by a ','
	108	//
	109	TString limits(bins);
	110	if (limits.IsNull()){
	111	Error("AliDielectronHelper::MakeArbitraryBinning","Bin Limit string is empty, cannot add the variable");
	112	return 0x0;
	113	}
	114
115	TObjArray *arr=limits.Tokenize(",");
116	Int_t nLimits=arr->GetEntries();
117	if (nLimits<2){
118	Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable");
119	delete arr;
120	return 0x0;
121	}
122
123	TVectorD *binLimits=new TVectorD(nLimits);
124	for (Int_t iLim=0; iLim<nLimits; ++iLim){
125	(binLimits)[iLim]=(static_cast<TObjString>(arr->At(iLim)))->GetString().Atof();
126	}
127
128	delete arr;
129	return binLimits;
130	}
131
ba15fdfb	132
	133	//_____________________________________________________________________________
	134	Int_t AliDielectronHelper::GetNch(const AliMCEvent *ev, Double_t etaRange){
	135	// determination of Nch
	136	if (!ev \|\| ev->IsA()!=AliMCEvent::Class()) return -1;
	137
	138	AliStack stack = ((AliMCEvent)ev)->Stack();
	139
	140	if (!stack) return -1;
	141
	142	Int_t nParticles = stack->GetNtrack();
	143	Int_t nCh = 0;
	144
	145	// count..
	146	for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
	147	if (!stack->IsPhysicalPrimary(iMc)) continue;
	148
	149	TParticle* particle = stack->Particle(iMc);
	150	if (!particle) continue;
	151	if (particle->GetPDG()->Charge() == 0) continue;
	152
	153	Float_t eta = particle->Eta();
	154	if (TMath::Abs(eta) < TMath::Abs(etaRange)) nCh++;
	155	}
	156
	157	return nCh;
	158	}
	159
	160
	161	//_____________________________________________________________________________
5720c765	162	Int_t AliDielectronHelper::GetNaccTrcklts(const AliVEvent *ev, Double_t etaRange){
ba15fdfb	163	// Compute the collision multiplicity based on AOD or ESD tracklets
	164	// Code taken from: AliAnalysisTaskMuonCollisionMultiplicity::ComputeMultiplicity()
	165
	166	if (!ev) return -1;
	167
	168	Int_t nTracklets = 0;
	169	Int_t nAcc = 0;
ba15fdfb	170
	171	if (ev->IsA() == AliAODEvent::Class()) {
	172	AliAODTracklets tracklets = ((AliAODEvent)ev)->GetTracklets();
0c09cae4	173	if (!tracklets) return -1;
ba15fdfb	174	nTracklets = tracklets->GetNumberOfTracklets();
	175	for (Int_t nn = 0; nn < nTracklets; nn++) {
	176	Double_t theta = tracklets->GetTheta(nn);
	177	Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
	178	if (TMath::Abs(eta) < etaRange) nAcc++;
	179	}
	180	} else if (ev->IsA() == AliESDEvent::Class()) {
	181	nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
	182	for (Int_t nn = 0; nn < nTracklets; nn++) {
	183	Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
	184	if (TMath::Abs(eta) < etaRange) nAcc++;
	185	}
	186	} else return -1;
	187
	188	return nAcc;
	189	}
	190
	191
5720c765	192
	193	//________________________________________________________________
	194	Double_t AliDielectronHelper::GetNaccTrckltsCorrected(const AliVEvent *event, Double_t uncorrectedNacc, Double_t vtxZ, Int_t type) {
	195	//
	196	// Correct the number of accepted tracklets based on the period and event vertex
	197	//
	198
	199	Int_t runNo = event->GetRunNumber();
	200
	201	Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
	202	Double_t refMult = 0.0; // reference multiplicity
	203
	204	if(runNo>114930 && runNo<117223) period = 0;
	205	if(runNo>119158 && runNo<120830) period = 1;
	206	if(runNo>122373 && runNo<126438) period = 2;
	207	if(runNo>127711 && runNo<130841) period = 3;
	208	if(period<0 \|\| period>3) return uncorrectedNacc;
	209
	210	if(type<0 \|\| type>8) return uncorrectedNacc;
	211	if(type == 0) refMult = 5.0; // SPD tracklets in \|eta\|<0.5
	212	if(type == 1) refMult = 9.5; // SPD tracklets in \|eta\|<1.0
	213	if(type == 2) refMult = 13.0; // SPD tracklets in \|eta\|<1.6
	214	if(type == 3) refMult = 6.0; // ITSTPC+ in \|eta\|<0.5
	215	if(type == 4) refMult = 12.0; // ITSTPC+ in \|eta\|<1.0
	216	if(type == 5) refMult = 16.0; // ITSTPC+ in \|eta\|<1.6
	217	if(type == 6) refMult = 6.0; // ITSSA+ in \|eta\|<0.5
	218	if(type == 7) refMult = 12.0; // ITSSA+ in \|eta\|<1.0
	219	if(type == 8) refMult = 15.0; // ITSSA+ in \|eta\|<1.6
	220
	221	if(TMath::Abs(vtxZ)>10.0) return uncorrectedNacc;
	222
	223	TProfile* estimatorAvg = AliDielectronVarManager::GetEstimatorHistogram(period, type);
	224	if(!estimatorAvg) return uncorrectedNacc;
	225
	226	Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
	227
	228	Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
	229
	230	Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
	231
	232	return correctedNacc;
	233	}
	234
ba15fdfb	235	//_____________________________________________________________________________
a823f01b	236	Int_t AliDielectronHelper::GetNacc(const AliVEvent *ev){
ba15fdfb	237	// put a robust Nacc definition here
ba15fdfb	238
a823f01b	239	if (!ev) return -1;
e4339752	240
a823f01b	241	AliDielectronVarCuts *varCuts = new AliDielectronVarCuts("VarCuts","VarCuts");
	242	varCuts->AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0);
	243	varCuts->AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0);
	244	varCuts->AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
e4339752	245	varCuts->AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0);
a823f01b	246	varCuts->AddCut(AliDielectronVarManager::kNclsTPC, 70.0, 160.0);
a823f01b	247	varCuts->AddCut(AliDielectronVarManager::kKinkIndex0, -0.5, 0.5); //noKinks
	248
	249	AliDielectronTrackCuts *trkCuts = new AliDielectronTrackCuts("TrkCuts","TrkCuts");
e4339752	250	trkCuts->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
a823f01b	251	trkCuts->SetRequireITSRefit(kTRUE);
a823f01b	252	trkCuts->SetRequireTPCRefit(kTRUE);
ba15fdfb	253
	254	Int_t nRecoTracks = ev->GetNumberOfTracks();
	255	Int_t nAcc = 0;
a823f01b	256
ba15fdfb	257	for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
7c0a2a4d	258	AliVTrack track = static_cast<AliVTrack>(ev->GetTrack(iTrack));
a823f01b	259	if (!track) continue;
41a64bfb	260	if (!trkCuts->IsSelected(track)) continue;
	261	if (!varCuts->IsSelected(track)) continue;
	262	nAcc++;
ba15fdfb	263	}
a823f01b	264
	265	delete varCuts;
	266	delete trkCuts;
	267
ba15fdfb	268	return nAcc;
	269	}
	270
41a64bfb	271	//_____________________________________________________________________________
	272	Double_t AliDielectronHelper::GetITSTPCMatchEff(const AliVEvent *ev){
	273	// recalulate the its-tpc matching efficiecy
	274
	275	if (!ev) return -1;
	276
	277	AliDielectronVarCuts *varCutsTPC = new AliDielectronVarCuts("VarCutsTPC","VarCutsTPC");
	278	varCutsTPC->AddCut(AliDielectronVarManager::kImpactParXY, -1.0, 1.0);
	279	varCutsTPC->AddCut(AliDielectronVarManager::kImpactParZ, -3.0, 3.0);
	280	varCutsTPC->AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
	281	varCutsTPC->AddCut(AliDielectronVarManager::kTPCchi2Cl, 0.0, 4.0);
	282	varCutsTPC->AddCut(AliDielectronVarManager::kNclsTPC, 50.0, 160.0);
	283	AliDielectronTrackCuts *trkCutsTPC = new AliDielectronTrackCuts("TrkCutsTPC","TrkCutsTPC");
	284	trkCutsTPC->SetRequireTPCRefit(kTRUE);
	285
	286	AliDielectronVarCuts *varCutsITS = new AliDielectronVarCuts("VarCutsITS","VarCutsITS");
	287	varCutsITS->AddCut(AliDielectronVarManager::kEta, -0.9, 0.9);
	288	AliDielectronTrackCuts *trkCutsITS = new AliDielectronTrackCuts("TrkCutsITS","TrkCutsITS");
	289	trkCutsITS->SetClusterRequirementITS(AliDielectronTrackCuts::kSPD, AliDielectronTrackCuts::kAny);
	290	trkCutsITS->SetRequireITSRefit(kTRUE);
	291
	292
	293	Int_t nRecoTracks = ev->GetNumberOfTracks();
	294	Double_t nTPC = 0, nITS = 0;
	295
	296	for (Int_t iTrack = 0; iTrack < nRecoTracks; iTrack++) {
	297	AliVTrack track = static_cast<AliVTrack>(ev->GetTrack(iTrack));
	298	if (!track) continue;
	299
	300	if(!trkCutsITS->IsSelected(track)) continue;
	301	if(!varCutsITS->IsSelected(track)) continue;
	302	nITS+=1.;
	303
	304	if(!trkCutsTPC->IsSelected(track)) continue;
	305	if(!varCutsTPC->IsSelected(track)) continue;
	306	nTPC+=1.;
	307
	308	}
	309
	310	delete varCutsITS;
	311	delete trkCutsITS;
	312	delete varCutsTPC;
	313	delete trkCutsTPC;
	314
	315	printf(" tracks TPC %.3e ITS %.3e = %.5f \n",nTPC,nITS,(nITS>0. ? nTPC/nITS : -1));
	316	return (nITS>0. ? nTPC/nITS : -1);
	317	}
	318
	319
1201a1a9	320	//_____________________________________________________________________________
	321	void AliDielectronHelper::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle, const AliVEvent * const ev){
	322	// Before rotate needs to be moved to position 0,0,0, ; move back after rotation
	323	if (!kfParticle) return;
	324	Double_t dx = 0.;
	325	Double_t dy = 0.;
	326	Double_t dz = 0.;
	327
	328	if (ev){
	329	dx = ev->GetPrimaryVertex()->GetX()-0.;
	330	dy = ev->GetPrimaryVertex()->GetY()-0.;
	331	dz = ev->GetPrimaryVertex()->GetZ()-0.;
	332	}
	333
	334	kfParticle->X() = kfParticle->GetX() - dx;
	335	kfParticle->Y() = kfParticle->GetY() - dy;
	336	kfParticle->Z() = kfParticle->GetZ() - dz;
	337
	338
	339	// Rotate the kf particle
	340	Double_t c = cos(angle);
	341	Double_t s = sin(angle);
	342
	343	Double_t mA[8][ 8];
	344	for( Int_t i=0; i<8; i++ ){
	345	for( Int_t j=0; j<8; j++){
	346	mA[i][j] = 0;
	347	}
	348	}
	349	for( int i=0; i<8; i++ ){
	350	mA[i][i] = 1;
	351	}
	352	mA[0][0] = c; mA[0][1] = s;
	353	mA[1][0] = -s; mA[1][1] = c;
	354	mA[3][3] = c; mA[3][4] = s;
	355	mA[4][3] = -s; mA[4][4] = c;
	356
	357	Double_t mAC[8][8];
	358	Double_t mAp[8];
	359
	360	for( Int_t i=0; i<8; i++ ){
	361	mAp[i] = 0;
	362	for( Int_t k=0; k<8; k++){
	363	mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
	364	}
	365	}
	366
	367	for( Int_t i=0; i<8; i++){
	368	kfParticle->Parameter(i) = mAp[i];
	369	}
	370
	371	for( Int_t i=0; i<8; i++ ){
	372	for( Int_t j=0; j<8; j++ ){
	373	mAC[i][j] = 0;
	374	for( Int_t k=0; k<8; k++ ){
	375	mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
	376	}
	377	}
	378	}
	379
	380	for( Int_t i=0; i<8; i++ ){
	381	for( Int_t j=0; j<=i; j++ ){
	382	Double_t xx = 0;
	383	for( Int_t k=0; k<8; k++){
384	xx+= mAC[i][k]*mA[j][k];
385	}
386	kfParticle->Covariance(i,j) = xx;
387	}
388	}
389
390	kfParticle->X() = kfParticle->GetX() + dx;
391	kfParticle->Y() = kfParticle->GetY() + dy;
392	kfParticle->Z() = kfParticle->GetZ() + dz;
393
394	}
395
5720c765	396	//_____________________________________________________________________________
5720c765	397	Int_t AliDielectronHelper::GetNMothers(const AliMCEvent *ev, Double_t etaRange, Int_t pdgMother, Int_t pdgDaughter, Int_t prim){
40875e45	398	// TODO: add AODs
5720c765	399	// counting number of mother particles generated in given eta range and 2 particle decay
	400	if (!ev \|\| ev->IsA()!=AliMCEvent::Class()) return -1;
	401
	402	AliStack stack = ((AliMCEvent)ev)->Stack();
	403
	404	if (!stack) return -1;
	405
	406	Int_t nParticles = stack->GetNtrack();
	407	Int_t nMothers = 0;
	408
	409	// count..
	410	for (Int_t iMc = 0; iMc < nParticles; ++iMc) {
	411
	412	TParticle* particle = stack->Particle(iMc);
	413	if (!particle) continue;
	414	if (particle->GetPdgCode() != pdgMother) continue;
	415	if (TMath::Abs(particle->Eta()) > TMath::Abs(etaRange)) continue;
	416
	417	if (particle->GetNDaughters() != 2) continue;
	418	// 1st daugther
	419	if (particle->GetFirstDaughter()>=nParticles \|\|
	420	particle->GetFirstDaughter()<0 ) continue;
	421
	422	TParticle* dau1 = stack->Particle(particle->GetFirstDaughter());
	423	if (TMath::Abs(dau1->GetPdgCode()) != pdgDaughter) continue;
	424	if (TMath::Abs(dau1->Eta()) > TMath::Abs(etaRange)) continue;
	425
	426	// 2nd daughter
	427	if (particle->GetLastDaughter()>=nParticles \|\|
	428	particle->GetLastDaughter()<0 ) continue;
	429
	430	TParticle* dau2 = stack->Particle(particle->GetLastDaughter());
	431	if (TMath::Abs(dau2->GetPdgCode()) != pdgDaughter) continue;
	432	if (TMath::Abs(dau2->Eta()) > TMath::Abs(etaRange)) continue;
	433
	434	// primary
	435	if (prim != -1) {
	436	if(particle->IsPrimary() != prim) continue;
	437	}
	438	nMothers++;
	439	}
	440	return nMothers;
	441	}
	442
	443