[u/mrichter/AliRoot.git] / PWGPP / TPC / AliTPCPIDEtaTree.cxx

#include "TChain.h"
#include "TTree.h"
#include "TF1.h"
#include "TAxis.h"
#include "TH2I.h"

#include "THnSparse.h"

#include "AliMCParticle.h"
//#include "AliStack.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDInputHandler.h"
#include "AliInputEventHandler.h"

#include "AliVVertex.h"
#include "AliAnalysisFilter.h"
#include "AliPID.h"
#include "AliPIDResponse.h"
#include "AliTPCPIDResponse.h"
//#include "AliTPCParamSR.h"

#include "AliTPCPIDEtaTree.h"

/*
This task determines the eta dependence of the TPC signal.
For this purpose, only tracks fulfilling some standard quality cuts are taken into account.
The obtained data can be used to derive the functional behaviour of the eta dependence.
Such a function can be plugged into this task to correct for the eta dependence and to see
if there is then really no eta dependence left.

Class written by Benjamin Hess.
Contact: bhess@cern.ch
*/

ClassImp(AliTPCPIDEtaTree)

//________________________________________________________________________
AliTPCPIDEtaTree::AliTPCPIDEtaTree()
  : AliTPCPIDBase()
  , fStoreMultiplicity(kFALSE)
  , fStoreNumOfSubthresholdclusters(kFALSE)
  , fStoreNumClustersInActiveVolume(kFALSE)
  , fDoAdditionalQA(kFALSE)
  , fPtpcPionCut(1.0)
  , fPtpc(0)
  , fPt(0)
  , fDeDx(0)
  , fDeDxExpected(0)
  , fTanTheta(0)
  //, fSinAlpha(0)
  //, fY(0)
  , fPhiPrime(0)
  , fTPCsignalN(0)
  , fTPCsignalNsubthreshold(0)
  , fNumTPCClustersInActiveVolume(0)
  , fPIDtype(0)
  , fMultiplicity(0)
  , fCorrectdEdxEtaDependence(kFALSE)
  , fCorrectdEdxMultiplicityDependence(kFALSE)
  , fTree(0x0)
  , fTreePions(0x0)
  , fOutputContainer(0x0)
  , fhTOFqa(0x0)
  , fhMultiplicityQA(0x0)
{
  // default Constructor
}

//________________________________________________________________________
AliTPCPIDEtaTree::AliTPCPIDEtaTree(const char *name)
  : AliTPCPIDBase(name)
  , fStoreMultiplicity(kFALSE)
  , fStoreNumOfSubthresholdclusters(kFALSE)
  , fStoreNumClustersInActiveVolume(kFALSE)
  , fDoAdditionalQA(kFALSE)
  , fPtpcPionCut(1.0)
  , fPtpc(0)
  , fPt(0)
  , fDeDx(0)
  , fDeDxExpected(0)
  , fTanTheta(0)
  //, fSinAlpha(0)
  //, fY(0)
  , fPhiPrime(0)
  , fTPCsignalN(0)
  , fTPCsignalNsubthreshold(0)
  , fNumTPCClustersInActiveVolume(0)
  , fPIDtype(0)
  , fMultiplicity(0)
  , fCorrectdEdxEtaDependence(kFALSE)
  , fCorrectdEdxMultiplicityDependence(kFALSE)
  , fTree(0x0)
  , fTreePions(0x0)
  , fOutputContainer(0x0)
  , fhTOFqa(0x0)
  , fhMultiplicityQA(0x0)
{
  // Constructor

  // Define input and output slots here
  DefineInput(0, TChain::Class());
  DefineOutput(1, TTree::Class());
  DefineOutput(2, TTree::Class());
  DefineOutput(3, TObjArray::Class());
}


//________________________________________________________________________
AliTPCPIDEtaTree::~AliTPCPIDEtaTree()
{
  // dtor
  
  //delete fTree;
  //fTree = 0x0;
  
  delete fOutputContainer;
  fOutputContainer = 0x0;
}


//________________________________________________________________________
void AliTPCPIDEtaTree::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
  
  AliTPCPIDBase::UserCreateOutputObjects();
  
  if (fDoAdditionalQA) {
    OpenFile(3);
    
    fOutputContainer = new TObjArray(2);
    fOutputContainer->SetName(GetName());
    fOutputContainer->SetOwner(kTRUE);
    
    const Int_t nBins = 6;
    // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
    Int_t bins[nBins] =    {   48,  48,    9,   30, 110,    4  };
    Double_t xmin[nBins] = {  0.3, 0.3, -0.9, -5.0, 0.0,   0.  };
    Double_t xmax[nBins] = {  2.7, 2.7,  0.9,  5.0, 1.1, 3200  };
    
    
    fhTOFqa = new THnSparseI("hTOFqa","", nBins, bins, xmin, xmax);
    fhTOFqa->GetAxis(0)->SetTitle("p_{Vertex} (GeV/c)");
    fhTOFqa->GetAxis(1)->SetTitle("p_{TPC_inner} (GeV/c)");
    fhTOFqa->GetAxis(2)->SetTitle("#eta");
    fhTOFqa->GetAxis(3)->SetTitle("n #sigma_{p} TOF");
    fhTOFqa->GetAxis(4)->SetTitle("#beta");
    fhTOFqa->GetAxis(5)->SetTitle("Multiplicity");
    
    fOutputContainer->Add(fhTOFqa);
    
    
    fhMultiplicityQA = new TH2I("hMultiplicityQA", "Multiplicity check; Contributors to primary vertex per event; Total number of tracks per Event",
                                120, 0, 6000, 400, 0, 20000);
    fOutputContainer->Add(fhMultiplicityQA);
  }
  else {
    fOutputContainer = new TObjArray(1);
    fOutputContainer->SetName(GetName());
    fOutputContainer->SetOwner(kTRUE);
  }
  
  OpenFile(1);
  
  fTree = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal");
  fTree->Branch("pTPC", &fPtpc);
  //fTree->Branch("pT", &fPt);
  fTree->Branch("dEdx", &fDeDx);
  fTree->Branch("dEdxExpected", &fDeDxExpected);
  fTree->Branch("tanTheta", &fTanTheta);
  //fTree->Branch("sinAlpha", &fSinAlpha);
  //fTree->Branch("y", &fY);
  //TODO fTree->Branch("phiPrime", &fPhiPrime);
  fTree->Branch("tpcSignalN", &fTPCsignalN);
  
  if (fStoreNumOfSubthresholdclusters)
    fTree->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
  
  if (fStoreNumClustersInActiveVolume)
    fTree->Branch("numTPCClustersInActiveVolume", &fNumTPCClustersInActiveVolume);
  
  fTree->Branch("pidType", &fPIDtype);
  
  if (fStoreMultiplicity)  {
    fTree->Branch("fMultiplicity", &fMultiplicity);
  }
  
  OpenFile(2);
  
  fTreePions = new TTree("fTreePions", "Tree for analysis of #eta dependence of TPC signal - Pions");
  fTreePions->Branch("pTPC", &fPtpc);
  fTreePions->Branch("pT", &fPt);
  fTreePions->Branch("dEdx", &fDeDx);
  fTreePions->Branch("dEdxExpected", &fDeDxExpected);
  fTreePions->Branch("tanTheta", &fTanTheta);
  fTreePions->Branch("tpcSignalN", &fTPCsignalN);
  
  if (fStoreNumOfSubthresholdclusters)
    fTreePions->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
  
  if (fStoreMultiplicity)  {
    fTreePions->Branch("fMultiplicity", &fMultiplicity);
  }
  
  PostData(1, fTree);
  PostData(2, fTreePions);
  PostData(3, fOutputContainer);
}


//________________________________________________________________________
void AliTPCPIDEtaTree::UserExec(Option_t *)
{
  // Main loop
  // Called for each event
  
  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }
  
  fMC = dynamic_cast<AliMCEvent*>(MCEvent());
  
  if (!fPIDResponse || !fV0KineCuts)
    return;
  
  if (!GetVertexIsOk(fESD))
    return;

  const AliVVertex* primaryVertex = fESD->GetPrimaryVertexTracks(); 
  if (!primaryVertex)
    return;
  //fMultiplicity = primaryVertex->GetNContributors();
  
  
  if(primaryVertex->GetNContributors() <= 0) 
    return;
  
  fMultiplicity = fESD->GetNumberOfTracks(); 
    
  if (fDoAdditionalQA) {
    Int_t nTotTracks = fESD->GetNumberOfTracks();
    fhMultiplicityQA->Fill(primaryVertex->GetNContributors(), nTotTracks);
  }
  
  Double_t magField = fESD->GetMagneticField();
  
  // Fill V0 arrays with V0s
  FillV0PIDlist();
  
  //AliTPCParamSR par;
  //par.Update();
  
  // Track loop to fill a Train spectrum
  for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
    Bool_t isPr = kFALSE;
    Bool_t isPi = kFALSE;
  
    AliESDtrack* track = fESD->GetTrack(iTracks);
    if (!track) {
      Printf("ERROR: Could not receive track %d", iTracks);
      continue;
    }
    
    if (TMath::Abs(track->Eta()) > fEtaCut)
      continue;
    
    // Do not distinguish positively and negatively charged V0's
    Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
    if (v0tagAllCharges == -99) {
      AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, fESD->GetNumberOfTracks(),
                    fV0tags != 0x0));
      continue;
    }
    
    Bool_t isV0prNotMC = (v0tagAllCharges == 16) && !fMC; // Only accept V0 protons for data, not for MC
    Bool_t isV0piNotMC = (v0tagAllCharges == 15) && !fMC; // Only accept V0 pions for data, not for MC
    
    Bool_t isV0NotMC = isV0prNotMC || isV0piNotMC;
    
    // Apply track cut
    if(!isV0NotMC && fTrackFilter && !fTrackFilter->IsSelected(track))
      continue;
    
    // Note: For V0's, the cut on ncl can be done via the tree (value is stored). One should not cut on geometry
    // for V0's since they have different topology
    if (!isV0NotMC && GetUseTPCCutMIGeo()) {
      if (!TPCCutMIGeo(track, InputEvent()))
        continue;
    }
    
    
    fPtpc = track->GetTPCmomentum();

    if (fPtpc > 5.)
      continue;
    
    fPt = track->Pt();
    fPhiPrime = GetPhiPrime(track->Phi(), magField, track->Charge());
    
    Int_t label = track->GetLabel();

    AliMCParticle* mcTrack = 0x0;
    
    if (fMC) {
      if (label < 0)
        continue; // If MC is available, reject tracks with negative ESD label
      mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
      if (!mcTrack) {
        Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
        continue;
      }
      
      /*// Only accept MC primaries
      if (!fMC->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
        continue;
      }*/
    }

    if (fMC) {
      Int_t pdg = mcTrack->PdgCode();
      
      if (TMath::Abs(pdg) == 2212) // Proton
        isPr = kTRUE;
      else if ((pdg == 111 || TMath::Abs(pdg) == 211) && fPtpc <= fPtpcPionCut) // Pion below desired momentum threshold
        isPi = kTRUE;
      else
        continue; // Only take protons and pions
      
      fPIDtype = kMCid;
      /*
      if (pdg == 111 || TMath::Abs(pdg) == 211) {//Pion
        binMC = 3;
      }
      else if (TMath::Abs(pdg) == 311 || TMath::Abs(pdg) == 321) {//Kaon
        binMC = 1;
      }
      else if (TMath::Abs(pdg) == 2212) {//Proton
        binMC = 4;
      }
      else if (TMath::Abs(pdg) == 11) {//Electron
        binMC = 0;
      }
      else if (TMath::Abs(pdg) == 13) {//Muon
        binMC = 2;
      }*/
    }
    
    fDeDx = track->GetTPCsignal();
    
    UInt_t status = track->GetStatus();
    Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
    Bool_t hasTOFtime = status&AliESDtrack::kTIME;
    Bool_t hasTOF     = kFALSE;
    if (hasTOFout && hasTOFtime)
      hasTOF = kTRUE;
    Float_t length = track->GetIntegratedLength();
    // Length check only for primaries, not for V0's!
    if (length < 350 && !isV0NotMC)
      hasTOF = kFALSE;
      
    if (!fMC) {    
      // Note: Normally, the track cuts include a cut on primaries. Therefore, if the particle is a V0, it would usually
      // NOT be found via the primary cuts. This means that the PIDtype describes more or less disjoint sets
      if (isV0prNotMC) {
        // V0
        isPr = kTRUE;
        
        if (!hasTOF) {
          fPIDtype = kV0idNoTOF;
        }
        else {
          if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0)
            fPIDtype = kV0idPlusTOFaccepted;
          else 
            fPIDtype = kV0idPlusTOFrejected;
        }
      }
      else if (isV0piNotMC) {
        // V0 pion
        if (fPtpc > fPtpcPionCut || fDeDx > 140.) // Reject pions above desired momentum threshold and also reject too high dEdx
          continue;
        
        isPi = kTRUE;
        
        if (!hasTOF) {
          fPIDtype = kV0idNoTOF;
        }
        else {
          if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) < 3.0)
            fPIDtype = kV0idPlusTOFaccepted;
          else  
            fPIDtype = kV0idPlusTOFrejected;
        }
      }
      else { 
        // Non-V0
        isPr = kTRUE; // If particle is accepted, it is a proton (for pions, only V0s are used)
        
        if (fPtpc < 4.0 && //TODO was 2.7 // Do not accept non-V0s above this threshold -> High contamination!!!
            (fDeDx >= 50. / TMath::Power(fPtpc, 1.3))) {// Pattern recognition instead of TPC cut to be ~independent of old TPC expected dEdx
          if (fPtpc < 0.6) {
            fPIDtype = kTPCid;
          }
          // fPtpc >= 0.6
          else if (hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0) {
              fPIDtype = kTPCandTOFid;
          }
          else
            continue; // Reject particle
        }
        else
          continue; // Reject particle
      }
    }
    
    if (fDoAdditionalQA) {
      Double_t tofTime = track->GetTOFsignal();//in ps
      Double_t tof = tofTime * 1E-3; // ns, average T0 fill subtracted, no info from T0detector   
      Double_t length2 = track->GetIntegratedLength();
      Double_t c = TMath::C() * 1.E-9;// m/ns
      length2 = length2 * 0.01; // in meters
      tof = tof * c;
      Double_t beta = length2 / tof;
      
      Double_t nSigmaTOF = hasTOF ? fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton) : 999;
      
      // p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
      Double_t entry[6] = { track->P(), fPtpc, track->Eta(), nSigmaTOF, beta, (Double_t)fMultiplicity };
      fhTOFqa->Fill(entry);
    }
    
    // Prepare entry for tree (some quantities have already been set)
    // Turn eta correction off -> Important here is the pure spline value and the selection via cuts. The correction
    // can be re-done manually, if needed. But it cannot be undone!
    
    if (isPr)
      fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, 
                                                                       fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
    else if (isPi)
      fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, 
                                                                       fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
    else
      fDeDxExpected = -1.;
    fTanTheta = track->GetInnerParam()->GetTgl();
    //fSinAlpha = track->GetInnerParam()->GetSnp();
    //fY = track->GetInnerParam()->GetY();
    fTPCsignalN = track->GetTPCsignalN();
    
    if (fStoreNumOfSubthresholdclusters) {
      const AliTPCdEdxInfo *clsInfo = track->GetTPCdEdxInfo();
      
      if (clsInfo) {
        Double32_t signal[4] = {0}; // 0: IROC only; 1: OROC short padas; 2: OROC long pads; 3:OROC combined
        Char_t ncl[3] = {0};        // clusters above threshold; 0: IROC; 1: OROC short; 2: OROC long
        Char_t nrows[3] = {0};      // clusters above and below threshold; 0: IROC; 1: OROC short; 2: OROC long

        clsInfo->GetTPCSignalRegionInfo(signal, ncl, nrows);

        fTPCsignalNsubthreshold = nrows[0] + nrows[1] + nrows[2] - ncl[0] - ncl[1] - ncl[2]; 
      }
      else
        fTPCsignalNsubthreshold = 200;// Set to invalid value
    }
    
    
    if (fStoreNumClustersInActiveVolume) 
      fNumTPCClustersInActiveVolume = track->GetLengthInActiveZone(1, 1.8, 220, magField);
    else
      fNumTPCClustersInActiveVolume = 200.;//Set to invalid value
 
    
    /* START
    const Double_t tanTheta = track->GetInnerParam()->GetTgl();

    // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
    // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
    const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
    const Double_t curvature = magField * constant / track->Pt(); // in 1./cm
    
    //const Double_t angleIROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(0,  par.GetNRow(0)  / 2.) * curvature) * 0.5, 1.));
    //const Double_t angleOROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRow(36) / 2.) * curvature) * 0.5, 1.));
    
    Double_t averageddzdr = 0.;
    Int_t nParts = 0;

    for (Double_t r = 85; r < 245; r++) {
      Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5);
      
      // Cut on |sin(phi)| as during reco
      if (TMath::Abs(sinPhiLocal) <= 0.95) {
        const Double_t phiLocal = TMath::ASin(sinPhiLocal);
        const Double_t tanPhiLocal = TMath::Tan(phiLocal);
        
        averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal); 
        nParts++;
      }
    }
    
    if (nParts > 0)
      averageddzdr /= nParts; 
    else {
      AliError("Problems during determination of dz/dr. Skipping track!");
      continue;
    }
    

    //printf("padRow 0 / 36: %f / %f\n", par.GetPadRowRadii(0,  par.GetNRow(0) / 2.), par.GetPadRowRadii(36, par.GetNRow(36) / 2.));
    //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\nIROC/OROC/averageOld/average: %f / %f / %f / //%f\ntanThetaGlobalFromTheta/tanTheta/dzdr: %f / %f / %f\n\n",
    //        track->Pt(), constant, magField, 1./curvature, angleIROC, angleOROC, angleAverage, averageAngle, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, dzdr);
    //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
    //        track->Pt(), constant, magField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);

  
    fTanTheta = averageddzdr;
    */
    
    if (isPr)
      fTree->Fill();
    else if (isPi)
      fTreePions->Fill();
  } //track loop 

  // Post output data.
  PostData(1, fTree);
  PostData(2, fTreePions);

  if (fDoAdditionalQA)
    PostData(3, fOutputContainer);
  
  // Clear the V0 PID arrays
  ClearV0PIDlist();
}      

//________________________________________________________________________
void AliTPCPIDEtaTree::Terminate(const Option_t *)
{
  // Called once at the end of the query
}
Commit	Line	Data
88b71b9f	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TF1.h"
	4	#include "TAxis.h"
	5	#include "TH2I.h"
	6
	7	#include "THnSparse.h"
	8
	9	#include "AliMCParticle.h"
	10	//#include "AliStack.h"
	11
	12	#include "AliAnalysisTask.h"
	13	#include "AliAnalysisManager.h"
	14
	15	#include "AliESDEvent.h"
	16	#include "AliMCEvent.h"
	17	#include "AliESDInputHandler.h"
	18	#include "AliInputEventHandler.h"
	19
	20	#include "AliVVertex.h"
	21	#include "AliAnalysisFilter.h"
	22	#include "AliPID.h"
	23	#include "AliPIDResponse.h"
	24	#include "AliTPCPIDResponse.h"
	25	//#include "AliTPCParamSR.h"
	26
	27	#include "AliTPCPIDEtaTree.h"
	28
	29	/*
	30	This task determines the eta dependence of the TPC signal.
	31	For this purpose, only tracks fulfilling some standard quality cuts are taken into account.
	32	The obtained data can be used to derive the functional behaviour of the eta dependence.
	33	Such a function can be plugged into this task to correct for the eta dependence and to see
	34	if there is then really no eta dependence left.
	35
	36	Class written by Benjamin Hess.
	37	Contact: bhess@cern.ch
	38	*/
	39
	40	ClassImp(AliTPCPIDEtaTree)
	41
	42	//________________________________________________________________________
	43	AliTPCPIDEtaTree::AliTPCPIDEtaTree()
caa5a9ed	44	: AliTPCPIDBase()
88b71b9f	45	, fStoreMultiplicity(kFALSE)
	46	, fStoreNumOfSubthresholdclusters(kFALSE)
	47	, fStoreNumClustersInActiveVolume(kFALSE)
	48	, fDoAdditionalQA(kFALSE)
	49	, fPtpcPionCut(1.0)
	50	, fPtpc(0)
	51	, fPt(0)
	52	, fDeDx(0)
	53	, fDeDxExpected(0)
	54	, fTanTheta(0)
	55	//, fSinAlpha(0)
	56	//, fY(0)
	57	, fPhiPrime(0)
	58	, fTPCsignalN(0)
	59	, fTPCsignalNsubthreshold(0)
	60	, fNumTPCClustersInActiveVolume(0)
	61	, fPIDtype(0)
	62	, fMultiplicity(0)
	63	, fCorrectdEdxEtaDependence(kFALSE)
	64	, fCorrectdEdxMultiplicityDependence(kFALSE)
	65	, fTree(0x0)
	66	, fTreePions(0x0)
	67	, fOutputContainer(0x0)
	68	, fhTOFqa(0x0)
	69	, fhMultiplicityQA(0x0)
	70	{
	71	// default Constructor
	72	}
	73
	74	//________________________________________________________________________
	75	AliTPCPIDEtaTree::AliTPCPIDEtaTree(const char *name)
caa5a9ed	76	: AliTPCPIDBase(name)
88b71b9f	77	, fStoreMultiplicity(kFALSE)
	78	, fStoreNumOfSubthresholdclusters(kFALSE)
	79	, fStoreNumClustersInActiveVolume(kFALSE)
	80	, fDoAdditionalQA(kFALSE)
	81	, fPtpcPionCut(1.0)
	82	, fPtpc(0)
	83	, fPt(0)
	84	, fDeDx(0)
	85	, fDeDxExpected(0)
	86	, fTanTheta(0)
	87	//, fSinAlpha(0)
	88	//, fY(0)
	89	, fPhiPrime(0)
	90	, fTPCsignalN(0)
	91	, fTPCsignalNsubthreshold(0)
	92	, fNumTPCClustersInActiveVolume(0)
	93	, fPIDtype(0)
	94	, fMultiplicity(0)
	95	, fCorrectdEdxEtaDependence(kFALSE)
	96	, fCorrectdEdxMultiplicityDependence(kFALSE)
	97	, fTree(0x0)
	98	, fTreePions(0x0)
	99	, fOutputContainer(0x0)
	100	, fhTOFqa(0x0)
	101	, fhMultiplicityQA(0x0)
	102	{
	103	// Constructor
	104
	105	// Define input and output slots here
	106	DefineInput(0, TChain::Class());
	107	DefineOutput(1, TTree::Class());
	108	DefineOutput(2, TTree::Class());
	109	DefineOutput(3, TObjArray::Class());
	110	}
	111
	112
	113	//________________________________________________________________________
	114	AliTPCPIDEtaTree::~AliTPCPIDEtaTree()
	115	{
	116	// dtor
	117
	118	//delete fTree;
	119	//fTree = 0x0;
	120
	121	delete fOutputContainer;
	122	fOutputContainer = 0x0;
	123	}
	124
	125
	126	//________________________________________________________________________
	127	void AliTPCPIDEtaTree::UserCreateOutputObjects()
	128	{
	129	// Create histograms
	130	// Called once
	131
caa5a9ed	132	AliTPCPIDBase::UserCreateOutputObjects();
88b71b9f	133
	134	if (fDoAdditionalQA) {
	135	OpenFile(3);
	136
	137	fOutputContainer = new TObjArray(2);
	138	fOutputContainer->SetName(GetName());
	139	fOutputContainer->SetOwner(kTRUE);
	140
	141	const Int_t nBins = 6;
	142	// p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
	143	Int_t bins[nBins] = { 48, 48, 9, 30, 110, 4 };
	144	Double_t xmin[nBins] = { 0.3, 0.3, -0.9, -5.0, 0.0, 0. };
	145	Double_t xmax[nBins] = { 2.7, 2.7, 0.9, 5.0, 1.1, 3200 };
	146
	147
	148	fhTOFqa = new THnSparseI("hTOFqa","", nBins, bins, xmin, xmax);
	149	fhTOFqa->GetAxis(0)->SetTitle("p_{Vertex} (GeV/c)");
	150	fhTOFqa->GetAxis(1)->SetTitle("p_{TPC_inner} (GeV/c)");
	151	fhTOFqa->GetAxis(2)->SetTitle("#eta");
	152	fhTOFqa->GetAxis(3)->SetTitle("n #sigma_{p} TOF");
	153	fhTOFqa->GetAxis(4)->SetTitle("#beta");
	154	fhTOFqa->GetAxis(5)->SetTitle("Multiplicity");
	155
	156	fOutputContainer->Add(fhTOFqa);
	157
	158
	159
	160	fhMultiplicityQA = new TH2I("hMultiplicityQA", "Multiplicity check; Contributors to primary vertex per event; Total number of tracks per Event",
	161	120, 0, 6000, 400, 0, 20000);
	162	fOutputContainer->Add(fhMultiplicityQA);
	163	}
	164	else {
	165	fOutputContainer = new TObjArray(1);
	166	fOutputContainer->SetName(GetName());
	167	fOutputContainer->SetOwner(kTRUE);
	168	}
	169
	170	OpenFile(1);
	171
	172	fTree = new TTree("fTree", "Tree for analysis of #eta dependence of TPC signal");
	173	fTree->Branch("pTPC", &fPtpc);
	174	//fTree->Branch("pT", &fPt);
	175	fTree->Branch("dEdx", &fDeDx);
	176	fTree->Branch("dEdxExpected", &fDeDxExpected);
	177	fTree->Branch("tanTheta", &fTanTheta);
	178	//fTree->Branch("sinAlpha", &fSinAlpha);
	179	//fTree->Branch("y", &fY);
	180	//TODO fTree->Branch("phiPrime", &fPhiPrime);
	181	fTree->Branch("tpcSignalN", &fTPCsignalN);
	182
	183	if (fStoreNumOfSubthresholdclusters)
	184	fTree->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
	185
	186	if (fStoreNumClustersInActiveVolume)
	187	fTree->Branch("numTPCClustersInActiveVolume", &fNumTPCClustersInActiveVolume);
	188
	189	fTree->Branch("pidType", &fPIDtype);
	190
	191	if (fStoreMultiplicity) {
	192	fTree->Branch("fMultiplicity", &fMultiplicity);
	193	}
	194
	195	OpenFile(2);
	196
197	fTreePions = new TTree("fTreePions", "Tree for analysis of #eta dependence of TPC signal - Pions");
198	fTreePions->Branch("pTPC", &fPtpc);
199	fTreePions->Branch("pT", &fPt);
200	fTreePions->Branch("dEdx", &fDeDx);
201	fTreePions->Branch("dEdxExpected", &fDeDxExpected);
202	fTreePions->Branch("tanTheta", &fTanTheta);
203	fTreePions->Branch("tpcSignalN", &fTPCsignalN);
204
205	if (fStoreNumOfSubthresholdclusters)
206	fTreePions->Branch("tpcSignalNsubthreshold", &fTPCsignalNsubthreshold);
207
208	if (fStoreMultiplicity) {
209	fTreePions->Branch("fMultiplicity", &fMultiplicity);
210	}
211
212	PostData(1, fTree);
213	PostData(2, fTreePions);
214	PostData(3, fOutputContainer);
215	}
216
217
218	//________________________________________________________________________
219	void AliTPCPIDEtaTree::UserExec(Option_t *)
220	{
221	// Main loop
222	// Called for each event
223
224	fESD = dynamic_cast<AliESDEvent*>(InputEvent());
225	if (!fESD) {
226	Printf("ERROR: fESD not available");
227	return;
228	}
229
230	fMC = dynamic_cast<AliMCEvent*>(MCEvent());
231
232	if (!fPIDResponse \|\| !fV0KineCuts)
233	return;
234
235	if (!GetVertexIsOk(fESD))
236	return;
237
238	const AliVVertex* primaryVertex = fESD->GetPrimaryVertexTracks();
239	if (!primaryVertex)
240	return;
241	//fMultiplicity = primaryVertex->GetNContributors();
242
243
244	if(primaryVertex->GetNContributors() <= 0)
245	return;
246
247	fMultiplicity = fESD->GetNumberOfTracks();
248
249	if (fDoAdditionalQA) {
250	Int_t nTotTracks = fESD->GetNumberOfTracks();
251	fhMultiplicityQA->Fill(primaryVertex->GetNContributors(), nTotTracks);
252	}
253
254	Double_t magField = fESD->GetMagneticField();
255
256	// Fill V0 arrays with V0s
257	FillV0PIDlist();
258
259	//AliTPCParamSR par;
260	//par.Update();
261
262	// Track loop to fill a Train spectrum
263	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
264	Bool_t isPr = kFALSE;
265	Bool_t isPi = kFALSE;
266
267	AliESDtrack* track = fESD->GetTrack(iTracks);
268	if (!track) {
269	Printf("ERROR: Could not receive track %d", iTracks);
270	continue;
271	}
272
273	if (TMath::Abs(track->Eta()) > fEtaCut)
274	continue;
275
276	// Do not distinguish positively and negatively charged V0's
277	Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
278	if (v0tagAllCharges == -99) {
279	AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, fESD->GetNumberOfTracks(),
280	fV0tags != 0x0));
281	continue;
282	}
283
284	Bool_t isV0prNotMC = (v0tagAllCharges == 16) && !fMC; // Only accept V0 protons for data, not for MC
285	Bool_t isV0piNotMC = (v0tagAllCharges == 15) && !fMC; // Only accept V0 pions for data, not for MC
286
287	Bool_t isV0NotMC = isV0prNotMC \|\| isV0piNotMC;
288
289	// Apply track cut
290	if(!isV0NotMC && fTrackFilter && !fTrackFilter->IsSelected(track))
291	continue;
292
293	// Note: For V0's, the cut on ncl can be done via the tree (value is stored). One should not cut on geometry
294	// for V0's since they have different topology
295	if (!isV0NotMC && GetUseTPCCutMIGeo()) {
296	if (!TPCCutMIGeo(track, InputEvent()))
297	continue;
298	}
299
300
301	fPtpc = track->GetTPCmomentum();
302
303	if (fPtpc > 5.)
304	continue;
305
306	fPt = track->Pt();
307	fPhiPrime = GetPhiPrime(track->Phi(), magField, track->Charge());
308
309	Int_t label = track->GetLabel();
310
311	AliMCParticle* mcTrack = 0x0;
312
313	if (fMC) {
314	if (label < 0)
315	continue; // If MC is available, reject tracks with negative ESD label
316	mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
317	if (!mcTrack) {
318	Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
319	continue;
320	}
321
322	/*// Only accept MC primaries
323	if (!fMC->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
324	continue;
325	}*/
326	}
327
328	if (fMC) {
329	Int_t pdg = mcTrack->PdgCode();
330
331	if (TMath::Abs(pdg) == 2212) // Proton
332	isPr = kTRUE;
333	else if ((pdg == 111 \|\| TMath::Abs(pdg) == 211) && fPtpc <= fPtpcPionCut) // Pion below desired momentum threshold
334	isPi = kTRUE;
335	else
336	continue; // Only take protons and pions
337
338	fPIDtype = kMCid;
339	/*
340	if (pdg == 111 \|\| TMath::Abs(pdg) == 211) {//Pion
341	binMC = 3;
342	}
343	else if (TMath::Abs(pdg) == 311 \|\| TMath::Abs(pdg) == 321) {//Kaon
344	binMC = 1;
345	}
346	else if (TMath::Abs(pdg) == 2212) {//Proton
347	binMC = 4;
348	}
349	else if (TMath::Abs(pdg) == 11) {//Electron
350	binMC = 0;
351	}
352	else if (TMath::Abs(pdg) == 13) {//Muon
353	binMC = 2;
354	}*/
355	}
356
357	fDeDx = track->GetTPCsignal();
358
359	UInt_t status = track->GetStatus();
360	Bool_t hasTOFout = status&AliESDtrack::kTOFout;
361	Bool_t hasTOFtime = status&AliESDtrack::kTIME;
362	Bool_t hasTOF = kFALSE;
363	if (hasTOFout && hasTOFtime)
364	hasTOF = kTRUE;
365	Float_t length = track->GetIntegratedLength();
366	// Length check only for primaries, not for V0's!
367	if (length < 350 && !isV0NotMC)
368	hasTOF = kFALSE;
369
370	if (!fMC) {
371	// Note: Normally, the track cuts include a cut on primaries. Therefore, if the particle is a V0, it would usually
372	// NOT be found via the primary cuts. This means that the PIDtype describes more or less disjoint sets
373	if (isV0prNotMC) {
374	// V0
375	isPr = kTRUE;
376
377	if (!hasTOF) {
378	fPIDtype = kV0idNoTOF;
379	}
380	else {
381	if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0)
382	fPIDtype = kV0idPlusTOFaccepted;
383	else
384	fPIDtype = kV0idPlusTOFrejected;
385	}
386	}
387	else if (isV0piNotMC) {
388	// V0 pion
389	if (fPtpc > fPtpcPionCut \|\| fDeDx > 140.) // Reject pions above desired momentum threshold and also reject too high dEdx
390	continue;
391
392	isPi = kTRUE;
393
394	if (!hasTOF) {
395	fPIDtype = kV0idNoTOF;
396	}
397	else {
398	if (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) < 3.0)
399	fPIDtype = kV0idPlusTOFaccepted;
400	else
401	fPIDtype = kV0idPlusTOFrejected;
402	}
403	}
404	else {
405	// Non-V0
406	isPr = kTRUE; // If particle is accepted, it is a proton (for pions, only V0s are used)
407
408	if (fPtpc < 4.0 && //TODO was 2.7 // Do not accept non-V0s above this threshold -> High contamination!!!
409	(fDeDx >= 50. / TMath::Power(fPtpc, 1.3))) {// Pattern recognition instead of TPC cut to be ~independent of old TPC expected dEdx
410	if (fPtpc < 0.6) {
411	fPIDtype = kTPCid;
412	}
413	// fPtpc >= 0.6
414	else if (hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0) {
415	fPIDtype = kTPCandTOFid;
416	}
417	else
418	continue; // Reject particle
419	}
420	else
421	continue; // Reject particle
422	}
423	}
424
425	if (fDoAdditionalQA) {
426	Double_t tofTime = track->GetTOFsignal();//in ps
427	Double_t tof = tofTime * 1E-3; // ns, average T0 fill subtracted, no info from T0detector
428	Double_t length2 = track->GetIntegratedLength();
429	Double_t c = TMath::C() * 1.E-9;// m/ns
430	length2 = length2 * 0.01; // in meters
431	tof = tof * c;
432	Double_t beta = length2 / tof;
433
434	Double_t nSigmaTOF = hasTOF ? fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton) : 999;
435
436	// p_vert, p_TPC, eta, nSigmaTOF, beta, multiplicity
0e60aeb1	437	Double_t entry[6] = { track->P(), fPtpc, track->Eta(), nSigmaTOF, beta, (Double_t)fMultiplicity };
88b71b9f	438	fhTOFqa->Fill(entry);
	439	}
	440
	441	// Prepare entry for tree (some quantities have already been set)
	442	// Turn eta correction off -> Important here is the pure spline value and the selection via cuts. The correction
	443	// can be re-done manually, if needed. But it cannot be undone!
	444
	445	if (isPr)
	446	fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
	447	fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
	448	else if (isPi)
	449	fDeDxExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
	450	fCorrectdEdxEtaDependence, fCorrectdEdxMultiplicityDependence);
	451	else
	452	fDeDxExpected = -1.;
	453	fTanTheta = track->GetInnerParam()->GetTgl();
	454	//fSinAlpha = track->GetInnerParam()->GetSnp();
	455	//fY = track->GetInnerParam()->GetY();
	456	fTPCsignalN = track->GetTPCsignalN();
	457
	458	if (fStoreNumOfSubthresholdclusters) {
	459	const AliTPCdEdxInfo *clsInfo = track->GetTPCdEdxInfo();
	460
	461	if (clsInfo) {
	462	Double32_t signal[4] = {0}; // 0: IROC only; 1: OROC short padas; 2: OROC long pads; 3:OROC combined
	463	Char_t ncl[3] = {0}; // clusters above threshold; 0: IROC; 1: OROC short; 2: OROC long
	464	Char_t nrows[3] = {0}; // clusters above and below threshold; 0: IROC; 1: OROC short; 2: OROC long
	465
	466	clsInfo->GetTPCSignalRegionInfo(signal, ncl, nrows);
	467
	468	fTPCsignalNsubthreshold = nrows[0] + nrows[1] + nrows[2] - ncl[0] - ncl[1] - ncl[2];
	469	}
	470	else
	471	fTPCsignalNsubthreshold = 200;// Set to invalid value
	472	}
	473
	474
	475	if (fStoreNumClustersInActiveVolume)
	476	fNumTPCClustersInActiveVolume = track->GetLengthInActiveZone(1, 1.8, 220, magField);
	477	else
	478	fNumTPCClustersInActiveVolume = 200.;//Set to invalid value
	479
	480
	481
	482	/* START
	483	const Double_t tanTheta = track->GetInnerParam()->GetTgl();
	484
	485	// Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss),
	486	// pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm)
	487	const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01;
	488	const Double_t curvature = magField * constant / track->Pt(); // in 1./cm
	489
	490	//const Double_t angleIROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(0, par.GetNRow(0) / 2.) * curvature) * 0.5, 1.));
	491	//const Double_t angleOROC = TMath::ASin(TMath::Min(TMath::Abs(par.GetPadRowRadii(36, par.GetNRow(36) / 2.) * curvature) * 0.5, 1.));
	492
	493	Double_t averageddzdr = 0.;
	494	Int_t nParts = 0;
	495
	496	for (Double_t r = 85; r < 245; r++) {
	497	Double_t sinPhiLocal = TMath::Abs(rcurvature0.5);
	498
	499	// Cut on \|sin(phi)\| as during reco
	500	if (TMath::Abs(sinPhiLocal) <= 0.95) {
	501	const Double_t phiLocal = TMath::ASin(sinPhiLocal);
502	const Double_t tanPhiLocal = TMath::Tan(phiLocal);
503
504	averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal);
505	nParts++;
506	}
507	}
508
509	if (nParts > 0)
510	averageddzdr /= nParts;
511	else {
512	AliError("Problems during determination of dz/dr. Skipping track!");
513	continue;
514	}
515
516
517	//printf("padRow 0 / 36: %f / %f\n", par.GetPadRowRadii(0, par.GetNRow(0) / 2.), par.GetPadRowRadii(36, par.GetNRow(36) / 2.));
518	//printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\nIROC/OROC/averageOld/average: %f / %f / %f / //%f\ntanThetaGlobalFromTheta/tanTheta/dzdr: %f / %f / %f\n\n",
519	// track->Pt(), constant, magField, 1./curvature, angleIROC, angleOROC, angleAverage, averageAngle, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, dzdr);
520	//printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n",
521	// track->Pt(), constant, magField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr);
522
523
524	fTanTheta = averageddzdr;
525	*/
526
527	if (isPr)
528	fTree->Fill();
529	else if (isPi)
530	fTreePions->Fill();
531	} //track loop
532
533	// Post output data.
534	PostData(1, fTree);
535	PostData(2, fTreePions);
536
537	if (fDoAdditionalQA)
538	PostData(3, fOutputContainer);
539
540	// Clear the V0 PID arrays
541	ClearV0PIDlist();
542	}
543
544	//________________________________________________________________________
545	void AliTPCPIDEtaTree::Terminate(const Option_t *)
546	{
547	// Called once at the end of the query
548	}