[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEreducedEventCreatorESD.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.                  *
**************************************************************************/
//
// Debug tree task
// the tree is represented as reduced events
// 
// Authors:
//   M.Fasel <M.Fasel@gsi.de>
//
//
#include <TArrayI.h>
#include <TBits.h>
#include <TFile.h>
#include <TParticle.h>
#include <TString.h>
#include <TTree.h>

#include "AliAnalysisManager.h"
#include "AliAODMCHeader.h"
#include "AliAnalysisUtils.h"
#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliInputEventHandler.h"
#include "AliHFEcuts.h"
#include "AliHFEextraCuts.h"
#include "AliHFEpidTPC.h"
#include "AliHFEreducedEvent.h"
#include "AliHFEreducedTrack.h"
#include "AliHFEreducedMCParticle.h"
#include "AliHFEsignalCuts.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
#include "AliMCEventHandler.h"
#include "AliMultiplicity.h"
#include "AliPIDResponse.h"
#include "AliVEvent.h"
#include "AliVTrack.h"
#include "AliVCluster.h"
#include "AliVVZERO.h"
#include "AliVZDC.h"
#include "AliTRDTriggerAnalysis.h"
#include "TTreeStream.h"

#include "AliHFEreducedEventCreatorESD.h"

ClassImp(AliHFEreducedEventCreatorESD)

AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD():
  AliAnalysisTaskSE(),
  fHFEtree(NULL),
  fAnalysisUtils(NULL),
  fHFEevent(NULL),
  fTrackCuts(NULL),
  fExtraCuts(NULL),
  fSignalCuts(NULL),
  fTPCpid(NULL),
  fTRDTriggerAnalysis(NULL),
  fEventNumber(0),
  fNclustersTPC(70),
  fNclustersTPCPID(0),
  fNclustersITS(2),
  fRemoveFirstEvent(kFALSE)
{
  //
  // Default constructor
  //
}

AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD(const char *name):
  AliAnalysisTaskSE(name),
  fHFEtree(NULL),
  fAnalysisUtils(NULL),
  fHFEevent(NULL),
  fTrackCuts(NULL),
  fExtraCuts(NULL),
  fSignalCuts(NULL),
  fTPCpid(NULL),
  fTRDTriggerAnalysis(NULL),
  fEventNumber(0),
  fNclustersTPC(70),
  fNclustersTPCPID(0),
  fNclustersITS(2),
  fRemoveFirstEvent(kFALSE)
{
  //
  // Default constructor
  //
  fTPCpid = new AliHFEpidTPC("QAtpcPID");
  fAnalysisUtils = new AliAnalysisUtils;
  fTRDTriggerAnalysis = new AliTRDTriggerAnalysis();
  fTRDTriggerAnalysis->SetRequireMatchElectron(kTRUE); 
  DefineOutput(1, TTree::Class());
}

AliHFEreducedEventCreatorESD::~AliHFEreducedEventCreatorESD(){
  //
  // Default destructor
  //
  if(fAnalysisUtils) delete fAnalysisUtils;
  if(fTPCpid) delete fTPCpid;
  if(fTRDTriggerAnalysis) delete fTRDTriggerAnalysis;
  if(fHFEevent) delete fHFEevent;
  if(fSignalCuts) delete fSignalCuts;
  if(fTrackCuts) delete fTrackCuts;

}

void AliHFEreducedEventCreatorESD::UserCreateOutputObjects(){
  //
  // Create debug tree, signal cuts and track cuts
  //

  //printf("test\n");

 // printf("testa\n");
  fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
  //printf("testb\n");
  
  fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
  fTrackCuts->CreateStandardCuts();
  // Track cuts
  fTrackCuts->SetMaxChi2perClusterTPC(1000.);   // switch off this cut (for the moment);
  fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
  fTrackCuts->SetMinRatioTPCclusters(0);
  fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable); 
  fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
  fTrackCuts->SetMinNClustersITS(fNclustersITS);
  // Event cuts
  fTrackCuts->SetUseMixedVertex(kTRUE);
  fTrackCuts->SetVertexRange(10.);
  //printf("testa\n");
  fTrackCuts->Initialize();
  //printf("testb\n");

  fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");

  fHFEevent = new AliHFEreducedEvent;
  OpenFile(1);
  fHFEtree = new TTree("HFEtree", "HFE event tree");
  fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
  PostData(1, fHFEtree);
}

void AliHFEreducedEventCreatorESD::UserExec(Option_t *){
  //
  // User Exec: Fill debug Tree
  // 

  // Get PID response
  AliPIDResponse *pid = NULL;
  AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if(handler){
//    printf("testb\n");
    pid = handler->GetPIDResponse();
  } else {
    AliError("No Handler");
  }
  if(!pid){
 //   printf("testc\n");
    AliError("No PID response");
    return;
  }
  if(!fInputEvent) {
    AliError("No Input event");
    return;
  }

  if(fRemoveFirstEvent){
    if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
  }

  AliDebug(1, "Event Selected");

  AliESDtrack copyTrack;

  // MC info
  Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
  if(mcthere){ 
    fTrackCuts->SetMCEvent(fMCEvent);
    fSignalCuts->SetMCEvent(fMCEvent);
  }  

  fTrackCuts->SetRecEvent(fInputEvent);

  if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
    AliDebug(1, "Event rejected by the event cuts\n");
    return;
  }
  
  // reject pile up in case of pp
  AliESDEvent *event = dynamic_cast<AliESDEvent *>(fInputEvent);
  if(event) {
    TString beamtype = event->GetBeamType();
    //printf("beamtype %s\n",(const char*)beamtype);
    if (strstr(beamtype,"p-p")) {
      //printf("Reject\n");
      if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
	      AliDebug(1, "Event flagged as pileup\n");
	      return;
      }
    }
  }

  if(!fExtraCuts){
    fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
  }
  fExtraCuts->SetRecEventInfo(fInputEvent);


  // Make Reduced Event 
  //AliHFEreducedEvent hfeevent;
  fHFEevent->~AliHFEreducedEvent();
  new(fHFEevent)AliHFEreducedEvent();

  // Get run number
  fHFEevent->SetRunNumber(fInputEvent->GetRunNumber());

  // Derive trigger 
  UInt_t trigger = fInputHandler->IsEventSelected();
  if(trigger & AliVEvent::kMB) fHFEevent->SetMBTrigger();
  if((trigger & AliVEvent::kINT7)||(trigger & AliVEvent::kINT8)) fHFEevent->SetINTTrigger();
  if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
  if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
  if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();

  fTRDTriggerAnalysis->CalcTriggers(event);
  if(fTRDTriggerAnalysis->HasTriggeredConfirmed(AliTRDTriggerAnalysis::kHSE)) fHFEevent->SetTRDSETrigger();
  if(fTRDTriggerAnalysis->HasTriggeredConfirmed(AliTRDTriggerAnalysis::kHQU)) fHFEevent->SetTRDDQTrigger();

  // Get Primary Vertex
  const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
  Double_t vtx[3];
  vertex->GetXYZ(vtx);
  fHFEevent->SetVX(vtx[0]);
  fHFEevent->SetVY(vtx[1]);
  fHFEevent->SetVZ(vtx[2]);
  Int_t ncontrib(vertex->GetNContributors());
  fHFEevent->SetNContribVertex(ncontrib);
  Double_t vcov[6];
  vertex->GetCovarianceMatrix(vcov);
  fHFEevent->SetVertexResolution(TMath::Sqrt(vcov[5]));
  // Get Primary Vertex from SPD
  const AliVVertex *vertexSPD = event->GetPrimaryVertexSPD();
  if(vertexSPD){
    memset(vtx, 0, sizeof(Double_t) *3);
    vertexSPD->GetXYZ(vtx);
    fHFEevent->SetVXSPD(vtx[0]);
    fHFEevent->SetVYSPD(vtx[1]);
    fHFEevent->SetVZSPD(vtx[2]);
    fHFEevent->SetNContribVertexSPD(vertexSPD->GetNContributors());
    memset(vcov, 0, sizeof(Double_t)*6);
    vertex->GetCovarianceMatrix(vcov);
    fHFEevent->SetVertexResolutionSPD(TMath::Sqrt(vcov[5]));
  }

  // Get centrality
  AliCentrality *hicent = fInputEvent->GetCentrality();
  fHFEevent->SetCentrality(
    hicent->GetCentralityPercentile("V0M"),
    hicent->GetCentralityPercentile("V0A"),
    hicent->GetCentralityPercentile("V0C"),
    hicent->GetCentralityPercentile("TKL"),
    hicent->GetCentralityPercentile("TRK"),
    hicent->GetCentralityPercentile("ZNA"),
    hicent->GetCentralityPercentile("ZNC"),
    hicent->GetCentralityPercentile("CL0"),
    hicent->GetCentralityPercentile("CL1"),
    hicent->GetCentralityPercentile("CND")
  );
  
  // Get VZERO Information
  AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
  if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());

  // Get ZDC Information
  AliVZDC *zdcinfo = fInputEvent->GetZDCData();
  if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy()); 

  // Set SPD multiplicity
  const AliMultiplicity *mult = event->GetMultiplicity();
  if(mult) fHFEevent->SetSPDMultiplicity(mult->GetNumberOfTracklets());

  //
  // Loop on MC tracks only
  //
  AliMCParticle *mctrack(NULL);
  // Monte-Carlo info
  Int_t source(5);
  if(mcthere){
    for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++) {
      mctrack = (AliMCParticle *)(fMCEvent->GetTrack(itrack));
      if(!mctrack) continue;
      AliHFEreducedMCParticle hfemcpart;
      if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack))  continue;        
      hfemcpart.SetSignal();
      // Kinematics
      hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
      hfemcpart.SetP(mctrack->P());
      hfemcpart.SetEta(mctrack->Eta());
      hfemcpart.SetPhi(mctrack->Phi());
      hfemcpart.SetPdg(mctrack->PdgCode());
      
      // Get Production Vertex in radial direction
      hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());

      // Get Mother PDG code of the particle
      Int_t motherlabel = TMath::Abs(mctrack->GetMother());
      if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
        AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
        if(mother) hfemcpart.SetMotherPdg(mother->PdgCode());
      }
      
      // derive source
      source = 5;
      if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
      else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
      else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
      else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
      else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
      else source = 5;
      hfemcpart.SetSource(source);

      fHFEevent->AddMCParticle(&hfemcpart);
    }
  }
  
  //
  // Loop on reconstructed tracks
  //
  TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
  Int_t counterdc=0;
  
  AliESDtrack *track = 0x0;
  for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
    // Run track loop
    track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
    if(!track) continue;
    // Cut track (Only basic track cuts)
    // printf("testv\n");
    if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
    //
    //printf("testu\n");

    // Kinematics
    AliHFEreducedTrack hfetrack;
    hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
    hfetrack.SetP(track->P());
    hfetrack.SetEta(track->Eta());
    hfetrack.SetPhi(track->Phi());
    hfetrack.SetTPCmomentum(track->GetTPCmomentum());

    // Track ID
    hfetrack.SetTrackID(track->GetID());

    // status
    ULong_t status = track->GetStatus();
    if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
    if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
    if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
    //if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
    if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
    Bool_t isEMCAL(kFALSE);
    Int_t fClsId = track->GetEMCALcluster();
    if(fClsId >= 0) isEMCAL = kTRUE;
    AliDebug(2, Form("cluster ID: %d, EMCAL: %s", fClsId, isEMCAL ? "yes" : "no"));
    if(isEMCAL) hfetrack.SetEMCALpid();
    // no filter bits available for ESDs

    // V0 information
    AliHFEreducedTrack::EV0PID_t v0pid = AliHFEreducedTrack::kV0undef;
    if(track->TestBit(BIT(14))) v0pid = AliHFEreducedTrack::kV0electron;
    else if(track->TestBit(BIT(15))) v0pid = AliHFEreducedTrack::kV0pion;
    else if(track->TestBit(BIT(16))) v0pid = AliHFEreducedTrack::kV0proton;
    hfetrack.SetV0PID(v0pid);

    if(mcthere){
      // Fill Monte-Carlo Information
      Int_t label = TMath::Abs(track->GetLabel());
      if(label < fMCEvent->GetNumberOfTracks())
        mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
      if(mctrack){
        AliDebug(2, "Associated MC particle found");
        if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
        // Kinematics
        hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
        hfetrack.SetMCP(mctrack->P());
        hfetrack.SetMCEta(mctrack->Eta());
        hfetrack.SetMCPhi(mctrack->Phi());
        hfetrack.SetMCPDG(mctrack->PdgCode());
      
        // Get Production Vertex in radial direction
        hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
      
        // Get Mother PDG code of the particle
        Int_t motherlabel = TMath::Abs(mctrack->GetMother());
        if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
          AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherlabel));
          if(mother) hfetrack.SetMCMotherPdg(mother->PdgCode());
        }
        
        // derive source
        source = 5;
        if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
        else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
        else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
        else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
        else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
        else source = 5;
        hfetrack.SetMCSource(source); 
      } else {
        AliDebug(2, "Associated MC particle not found");
      }
    }

    // HFE DCA
    Float_t dcaxy = -999.,
            dcaz = -999.;
    fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
    hfetrack.SetDCA(dcaxy, dcaz);
    Double_t hfeImpactParam(-999.), hfeImpactParamResol(-999.);
    fExtraCuts->GetHFEImpactParameters((AliVTrack *)track,hfeImpactParam,hfeImpactParamResol);
    hfetrack.SetHFEImpactParam(hfeImpactParam,hfeImpactParamResol);

    // Different number of clusters definitions
    Int_t nclustersITS(track->GetITSclusters(NULL)),
          nclustersTPC(track->GetTPCNcls()),
          nclustersTPCall(track->GetTPCClusterMap().CountBits()),
          nclustersTPCshared(0);
    UChar_t nfindableTPC = track->GetTPCNclsF();
    const TBits &sharedTPC = track->GetTPCSharedMap();
    for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
    hfetrack.SetChi2PerTPCcluster(track->GetTPCchi2()/Double_t(nclustersTPC));
    hfetrack.SetITSnclusters(nclustersITS);
    hfetrack.SetTPCnclusters(nclustersTPC);
    hfetrack.SetTRDnclusters(track->GetTRDncls());
    hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
    hfetrack.SetTPCcrossedRows(track->GetTPCCrossedRows());
    hfetrack.SetTPCnclustersAll(nclustersTPCall);
    hfetrack.SetTPCsharedClusters(nclustersTPCshared);
    hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
    hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
    UChar_t itsPixel = track->GetITSClusterMap();
    for(int ily = 0; ily < 6; ily++) 
            if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
   
    // TRD related quantities (Yvonne)
    Int_t nslices = track->GetNumberOfTRDslices();
    hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
    hfetrack.SetTRDnslices(nslices);
    hfetrack.SetTRDchi2(track->GetTRDchi2());
    Int_t nslicetemp=0;
    for(Int_t iplane = 0; iplane < 6; iplane++){
	    nslicetemp=0;
	    for(Int_t isl = 0; isl < nslices; isl++){
	      if(track->GetTRDntrackletsPID()>0){
		      if(track->GetTRDslice(iplane, isl)>0.001) nslicetemp++;
	      }
	    }
	    if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
    }


    //test for kink tracks
    if(fExtraCuts->IsKinkMother(track)) hfetrack.SetIsKinkMother();
    else if(fExtraCuts->IsKinkDaughter(track)) hfetrack.SetIsKinkDaughter();
    
    // Double counted
    Int_t id(track->GetID());
    for(Int_t l=0; l < counterdc; l++){
      Int_t iTrack2 = arraytrack.At(l);
      if(iTrack2==id){
         hfetrack.SetDoubleCounted();
         break;
      }
    }
    // Add the id at this place
    arraytrack.AddAt(id,counterdc);
    counterdc++;

    // PID
    hfetrack.SetTPCdEdx(track->GetTPCsignal());
    hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
    hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
    hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
    hfetrack.SetITSsigmaEl(pid->NumberOfSigmasITS(track, AliPID::kElectron));
    // Eta correction
    copyTrack.~AliESDtrack();
    new(&copyTrack) AliESDtrack(*track);
    if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kESDanalysis);
    if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kESDanalysis);
    hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
    hfetrack.SetTPCdEdxCorrected(copyTrack.GetTPCsignal());
    if(isEMCAL){
      AliDebug(2, "Adding EMCAL PID information");
      // EMCAL cluster
      Double_t emcalEnergyOverP = -1.,
               showershape[4] = {0.,0.,0.,0.};
      hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
      hfetrack.SetEMCALEoverP(emcalEnergyOverP);
      hfetrack.SetEMCALShowerShape(showershape);
    }

    // Track finished, add NOW to the Event
    fHFEevent->AddTrack(&hfetrack);
    //printf("after\n");
  }
  
  // Fill the debug tree
  //AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
  //AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
  fHFEtree->Fill();

  fEventNumber++;
  PostData(1, fHFEtree);
}

void AliHFEreducedEventCreatorESD::Terminate(Option_t *){
  //
  // Terminate
  //
  AliInfo("terminating...\n");

}

Bool_t AliHFEreducedEventCreatorESD::IsTOFmismatch(const AliVTrack *const track, const AliPIDResponse *const pid) const {
  //
  // Is TOF mismatch
  //
  Double_t probs[AliPID::kSPECIESC];
  AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
  return status == AliPIDResponse::kDetMismatch;
}
Commit	Line	Data
3513afb7	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	// Debug tree task
	17	// the tree is represented as reduced events
	18	//
	19	// Authors:
	20	// M.Fasel <M.Fasel@gsi.de>
	21	//
	22	//
	23	#include <TArrayI.h>
	24	#include <TBits.h>
	25	#include <TFile.h>
	26	#include <TParticle.h>
	27	#include <TString.h>
	28	#include <TTree.h>
	29
	30	#include "AliAnalysisManager.h"
	31	#include "AliAODMCHeader.h"
	32	#include "AliAnalysisUtils.h"
	33	#include "AliCentrality.h"
	34	#include "AliESDEvent.h"
	35	#include "AliESDtrack.h"
	36	#include "AliESDVertex.h"
	37	#include "AliInputEventHandler.h"
	38	#include "AliHFEcuts.h"
	39	#include "AliHFEextraCuts.h"
	40	#include "AliHFEpidTPC.h"
	41	#include "AliHFEreducedEvent.h"
	42	#include "AliHFEreducedTrack.h"
	43	#include "AliHFEreducedMCParticle.h"
	44	#include "AliHFEsignalCuts.h"
	45	#include "AliLog.h"
	46	#include "AliMCEvent.h"
	47	#include "AliMCParticle.h"
	48	#include "AliMCEventHandler.h"
	49	#include "AliMultiplicity.h"
	50	#include "AliPIDResponse.h"
	51	#include "AliVEvent.h"
	52	#include "AliVTrack.h"
	53	#include "AliVCluster.h"
	54	#include "AliVVZERO.h"
	55	#include "AliVZDC.h"
618038b6	56	#include "AliTRDTriggerAnalysis.h"
3513afb7	57	#include "TTreeStream.h"
	58
	59	#include "AliHFEreducedEventCreatorESD.h"
	60
	61	ClassImp(AliHFEreducedEventCreatorESD)
	62
	63	AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD():
	64	AliAnalysisTaskSE(),
	65	fHFEtree(NULL),
	66	fAnalysisUtils(NULL),
	67	fHFEevent(NULL),
	68	fTrackCuts(NULL),
	69	fExtraCuts(NULL),
	70	fSignalCuts(NULL),
	71	fTPCpid(NULL),
618038b6	72	fTRDTriggerAnalysis(NULL),
3513afb7	73	fEventNumber(0),
	74	fNclustersTPC(70),
	75	fNclustersTPCPID(0),
	76	fNclustersITS(2),
	77	fRemoveFirstEvent(kFALSE)
	78	{
	79	//
	80	// Default constructor
	81	//
	82	}
	83
	84	AliHFEreducedEventCreatorESD::AliHFEreducedEventCreatorESD(const char *name):
	85	AliAnalysisTaskSE(name),
	86	fHFEtree(NULL),
	87	fAnalysisUtils(NULL),
	88	fHFEevent(NULL),
	89	fTrackCuts(NULL),
	90	fExtraCuts(NULL),
	91	fSignalCuts(NULL),
	92	fTPCpid(NULL),
618038b6	93	fTRDTriggerAnalysis(NULL),
3513afb7	94	fEventNumber(0),
	95	fNclustersTPC(70),
	96	fNclustersTPCPID(0),
	97	fNclustersITS(2),
	98	fRemoveFirstEvent(kFALSE)
	99	{
	100	//
	101	// Default constructor
	102	//
	103	fTPCpid = new AliHFEpidTPC("QAtpcPID");
	104	fAnalysisUtils = new AliAnalysisUtils;
4437a0d2	105	fTRDTriggerAnalysis = new AliTRDTriggerAnalysis();
58d0bf42	106	fTRDTriggerAnalysis->SetRequireMatchElectron(kTRUE);
3513afb7	107	DefineOutput(1, TTree::Class());
	108	}
	109
	110	AliHFEreducedEventCreatorESD::~AliHFEreducedEventCreatorESD(){
	111	//
	112	// Default destructor
	113	//
	114	if(fAnalysisUtils) delete fAnalysisUtils;
	115	if(fTPCpid) delete fTPCpid;
618038b6	116	if(fTRDTriggerAnalysis) delete fTRDTriggerAnalysis;
3513afb7	117	if(fHFEevent) delete fHFEevent;
	118	if(fSignalCuts) delete fSignalCuts;
	119	if(fTrackCuts) delete fTrackCuts;
4437a0d2	120
3513afb7	121	}
	122
	123	void AliHFEreducedEventCreatorESD::UserCreateOutputObjects(){
	124	//
	125	// Create debug tree, signal cuts and track cuts
	126	//
	127
	128	//printf("test\n");
	129
	130	// printf("testa\n");
	131	fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
	132	//printf("testb\n");
	133
	134	fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
	135	fTrackCuts->CreateStandardCuts();
	136	// Track cuts
4437a0d2	137	fTrackCuts->SetMaxChi2perClusterTPC(1000.); // switch off this cut (for the moment);
3513afb7	138	fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
	139	fTrackCuts->SetMinRatioTPCclusters(0);
	140	fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
	141	fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
	142	fTrackCuts->SetMinNClustersITS(fNclustersITS);
	143	// Event cuts
	144	fTrackCuts->SetUseMixedVertex(kTRUE);
	145	fTrackCuts->SetVertexRange(10.);
	146	//printf("testa\n");
	147	fTrackCuts->Initialize();
	148	//printf("testb\n");
	149
	150	fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
	151
	152	fHFEevent = new AliHFEreducedEvent;
	153	OpenFile(1);
	154	fHFEtree = new TTree("HFEtree", "HFE event tree");
	155	fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
	156	PostData(1, fHFEtree);
	157	}
	158
	159	void AliHFEreducedEventCreatorESD::UserExec(Option_t *){
	160	//
	161	// User Exec: Fill debug Tree
	162	//
	163
	164	// Get PID response
	165	AliPIDResponse *pid = NULL;
	166	AliInputEventHandler handler = dynamic_cast<AliInputEventHandler >(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	167	if(handler){
	168	// printf("testb\n");
	169	pid = handler->GetPIDResponse();
	170	} else {
	171	AliError("No Handler");
	172	}
	173	if(!pid){
	174	// printf("testc\n");
	175	AliError("No PID response");
	176	return;
	177	}
	178	if(!fInputEvent) {
	179	AliError("No Input event");
	180	return;
	181	}
	182
	183	if(fRemoveFirstEvent){
	184	if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
	185	}
	186
	187	AliDebug(1, "Event Selected");
	188
	189	AliESDtrack copyTrack;
	190
	191	// MC info
	192	Bool_t mcthere = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) != NULL;
	193	if(mcthere){
	194	fTrackCuts->SetMCEvent(fMCEvent);
	195	fSignalCuts->SetMCEvent(fMCEvent);
	196	}
	197
	198	fTrackCuts->SetRecEvent(fInputEvent);
	199
	200	if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
	201	AliDebug(1, "Event rejected by the event cuts\n");
202	return;
203	}
204
205	// reject pile up in case of pp
206	AliESDEvent event = dynamic_cast<AliESDEvent >(fInputEvent);
207	if(event) {
208	TString beamtype = event->GetBeamType();
209	//printf("beamtype %s\n",(const char*)beamtype);
210	if (strstr(beamtype,"p-p")) {
211	//printf("Reject\n");
212	if(fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5)){
7e695443	213	AliDebug(1, "Event flagged as pileup\n");
7e695443	214	return;
3513afb7	215	}
	216	}
	217	}
	218
	219	if(!fExtraCuts){
	220	fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
	221	}
	222	fExtraCuts->SetRecEventInfo(fInputEvent);
	223
	224
	225	// Make Reduced Event
	226	//AliHFEreducedEvent hfeevent;
	227	fHFEevent->~AliHFEreducedEvent();
	228	new(fHFEevent)AliHFEreducedEvent();
	229
	230	// Get run number
	231	fHFEevent->SetRunNumber(fInputEvent->GetRunNumber());
	232
	233	// Derive trigger
	234	UInt_t trigger = fInputHandler->IsEventSelected();
	235	if(trigger & AliVEvent::kMB) fHFEevent->SetMBTrigger();
420bd2ea	236	if((trigger & AliVEvent::kINT7)\|\|(trigger & AliVEvent::kINT8)) fHFEevent->SetINTTrigger();
3513afb7	237	if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
	238	if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
	239	if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();
4437a0d2	240
618038b6	241	fTRDTriggerAnalysis->CalcTriggers(event);
58d0bf42	242	if(fTRDTriggerAnalysis->HasTriggeredConfirmed(AliTRDTriggerAnalysis::kHSE)) fHFEevent->SetTRDSETrigger();
58d0bf42	243	if(fTRDTriggerAnalysis->HasTriggeredConfirmed(AliTRDTriggerAnalysis::kHQU)) fHFEevent->SetTRDDQTrigger();
3513afb7	244
	245	// Get Primary Vertex
	246	const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
	247	Double_t vtx[3];
	248	vertex->GetXYZ(vtx);
	249	fHFEevent->SetVX(vtx[0]);
	250	fHFEevent->SetVY(vtx[1]);
	251	fHFEevent->SetVZ(vtx[2]);
7e695443	252	Int_t ncontrib(vertex->GetNContributors());
3513afb7	253	fHFEevent->SetNContribVertex(ncontrib);
7e695443	254	Double_t vcov[6];
	255	vertex->GetCovarianceMatrix(vcov);
	256	fHFEevent->SetVertexResolution(TMath::Sqrt(vcov[5]));
	257	// Get Primary Vertex from SPD
	258	const AliVVertex *vertexSPD = event->GetPrimaryVertexSPD();
	259	if(vertexSPD){
	260	memset(vtx, 0, sizeof(Double_t) *3);
	261	vertexSPD->GetXYZ(vtx);
	262	fHFEevent->SetVXSPD(vtx[0]);
	263	fHFEevent->SetVYSPD(vtx[1]);
	264	fHFEevent->SetVZSPD(vtx[2]);
	265	fHFEevent->SetNContribVertexSPD(vertexSPD->GetNContributors());
	266	memset(vcov, 0, sizeof(Double_t)*6);
	267	vertex->GetCovarianceMatrix(vcov);
	268	fHFEevent->SetVertexResolutionSPD(TMath::Sqrt(vcov[5]));
	269	}
3513afb7	270
	271	// Get centrality
	272	AliCentrality *hicent = fInputEvent->GetCentrality();
	273	fHFEevent->SetCentrality(
	274	hicent->GetCentralityPercentile("V0M"),
	275	hicent->GetCentralityPercentile("V0A"),
	276	hicent->GetCentralityPercentile("V0C"),
	277	hicent->GetCentralityPercentile("TKL"),
	278	hicent->GetCentralityPercentile("TRK"),
4437a0d2	279	hicent->GetCentralityPercentile("ZNA"),
	280	hicent->GetCentralityPercentile("ZNC"),
	281	hicent->GetCentralityPercentile("CL0"),
	282	hicent->GetCentralityPercentile("CL1"),
	283	hicent->GetCentralityPercentile("CND")
3513afb7	284	);
	285
	286	// Get VZERO Information
	287	AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
	288	if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());
	289
	290	// Get ZDC Information
	291	AliVZDC *zdcinfo = fInputEvent->GetZDCData();
	292	if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy());
	293
	294	// Set SPD multiplicity
	295	const AliMultiplicity *mult = event->GetMultiplicity();
	296	if(mult) fHFEevent->SetSPDMultiplicity(mult->GetNumberOfTracklets());
	297
	298	//
	299	// Loop on MC tracks only
	300	//
	301	AliMCParticle *mctrack(NULL);
	302	// Monte-Carlo info
	303	Int_t source(5);
	304	if(mcthere){
	305	for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++) {
	306	mctrack = (AliMCParticle *)(fMCEvent->GetTrack(itrack));
	307	if(!mctrack) continue;
	308	AliHFEreducedMCParticle hfemcpart;
4437a0d2	309	if(!fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) continue;
4437a0d2	310	hfemcpart.SetSignal();
3513afb7	311	// Kinematics
	312	hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
	313	hfemcpart.SetP(mctrack->P());
	314	hfemcpart.SetEta(mctrack->Eta());
	315	hfemcpart.SetPhi(mctrack->Phi());
	316	hfemcpart.SetPdg(mctrack->PdgCode());
	317
	318	// Get Production Vertex in radial direction
	319	hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
	320
	321	// Get Mother PDG code of the particle
	322	Int_t motherlabel = TMath::Abs(mctrack->GetMother());
	323	if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
	324	AliMCParticle mother = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(motherlabel));
	325	if(mother) hfemcpart.SetMotherPdg(mother->PdgCode());
	326	}
	327
	328	// derive source
	329	source = 5;
	330	if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
	331	else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
	332	else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
	333	else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
	334	else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
	335	else source = 5;
	336	hfemcpart.SetSource(source);
	337
	338	fHFEevent->AddMCParticle(&hfemcpart);
	339	}
	340	}
	341
	342	//
	343	// Loop on reconstructed tracks
	344	//
	345	TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
	346	Int_t counterdc=0;
	347
	348	AliESDtrack *track = 0x0;
	349	for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
	350	// Run track loop
	351	track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrack));
	352	if(!track) continue;
	353	// Cut track (Only basic track cuts)
	354	// printf("testv\n");
	355	if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
	356	//
	357	//printf("testu\n");
	358
	359	// Kinematics
	360	AliHFEreducedTrack hfetrack;
	361	hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
	362	hfetrack.SetP(track->P());
	363	hfetrack.SetEta(track->Eta());
	364	hfetrack.SetPhi(track->Phi());
	365	hfetrack.SetTPCmomentum(track->GetTPCmomentum());
	366
	367	// Track ID
	368	hfetrack.SetTrackID(track->GetID());
	369
	370	// status
	371	ULong_t status = track->GetStatus();
	372	if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
	373	if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
	374	if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
375	//if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
376	if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
4437a0d2	377	Bool_t isEMCAL(kFALSE);
	378	Int_t fClsId = track->GetEMCALcluster();
	379	if(fClsId >= 0) isEMCAL = kTRUE;
	380	AliDebug(2, Form("cluster ID: %d, EMCAL: %s", fClsId, isEMCAL ? "yes" : "no"));
	381	if(isEMCAL) hfetrack.SetEMCALpid();
3513afb7	382	// no filter bits available for ESDs
3513afb7	383
4437a0d2	384	// V0 information
	385	AliHFEreducedTrack::EV0PID_t v0pid = AliHFEreducedTrack::kV0undef;
	386	if(track->TestBit(BIT(14))) v0pid = AliHFEreducedTrack::kV0electron;
	387	else if(track->TestBit(BIT(15))) v0pid = AliHFEreducedTrack::kV0pion;
	388	else if(track->TestBit(BIT(16))) v0pid = AliHFEreducedTrack::kV0proton;
	389	hfetrack.SetV0PID(v0pid);
	390
3513afb7	391	if(mcthere){
	392	// Fill Monte-Carlo Information
	393	Int_t label = TMath::Abs(track->GetLabel());
4437a0d2	394	if(label < fMCEvent->GetNumberOfTracks())
3513afb7	395	mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
4437a0d2	396	if(mctrack){
	397	AliDebug(2, "Associated MC particle found");
	398	if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
	399	// Kinematics
	400	hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
	401	hfetrack.SetMCP(mctrack->P());
	402	hfetrack.SetMCEta(mctrack->Eta());
	403	hfetrack.SetMCPhi(mctrack->Phi());
	404	hfetrack.SetMCPDG(mctrack->PdgCode());
3513afb7	405
4437a0d2	406	// Get Production Vertex in radial direction
4437a0d2	407	hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
3513afb7	408
4437a0d2	409	// Get Mother PDG code of the particle
	410	Int_t motherlabel = TMath::Abs(mctrack->GetMother());
	411	if(motherlabel >= 0 && motherlabel < fMCEvent->GetNumberOfTracks()){
	412	AliMCParticle mother = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(motherlabel));
	413	if(mother) hfetrack.SetMCMotherPdg(mother->PdgCode());
	414	}
	415
	416	// derive source
	417	source = 5;
	418	if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
	419	else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
	420	else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
	421	else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
	422	else if(TMath::Abs(mctrack->PdgCode()) == 11) source = 4;
	423	else source = 5;
	424	hfetrack.SetMCSource(source);
	425	} else {
	426	AliDebug(2, "Associated MC particle not found");
3513afb7	427	}
3513afb7	428	}
	429
	430	// HFE DCA
	431	Float_t dcaxy = -999.,
	432	dcaz = -999.;
	433	fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
	434	hfetrack.SetDCA(dcaxy, dcaz);
4437a0d2	435	Double_t hfeImpactParam(-999.), hfeImpactParamResol(-999.);
	436	fExtraCuts->GetHFEImpactParameters((AliVTrack *)track,hfeImpactParam,hfeImpactParamResol);
	437	hfetrack.SetHFEImpactParam(hfeImpactParam,hfeImpactParamResol);
3513afb7	438
	439	// Different number of clusters definitions
	440	Int_t nclustersITS(track->GetITSclusters(NULL)),
	441	nclustersTPC(track->GetTPCNcls()),
	442	nclustersTPCall(track->GetTPCClusterMap().CountBits()),
	443	nclustersTPCshared(0);
	444	UChar_t nfindableTPC = track->GetTPCNclsF();
	445	const TBits &sharedTPC = track->GetTPCSharedMap();
	446	for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
4437a0d2	447	hfetrack.SetChi2PerTPCcluster(track->GetTPCchi2()/Double_t(nclustersTPC));
3513afb7	448	hfetrack.SetITSnclusters(nclustersITS);
	449	hfetrack.SetTPCnclusters(nclustersTPC);
	450	hfetrack.SetTRDnclusters(track->GetTRDncls());
	451	hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
	452	hfetrack.SetTPCcrossedRows(track->GetTPCCrossedRows());
	453	hfetrack.SetTPCnclustersAll(nclustersTPCall);
	454	hfetrack.SetTPCsharedClusters(nclustersTPCshared);
	455	hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
	456	hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
	457	UChar_t itsPixel = track->GetITSClusterMap();
	458	for(int ily = 0; ily < 6; ily++)
	459	if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
	460
	461	// TRD related quantities (Yvonne)
	462	Int_t nslices = track->GetNumberOfTRDslices();
	463	hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
	464	hfetrack.SetTRDnslices(nslices);
	465	hfetrack.SetTRDchi2(track->GetTRDchi2());
	466	Int_t nslicetemp=0;
	467	for(Int_t iplane = 0; iplane < 6; iplane++){
	468	nslicetemp=0;
	469	for(Int_t isl = 0; isl < nslices; isl++){
	470	if(track->GetTRDntrackletsPID()>0){
	471	if(track->GetTRDslice(iplane, isl)>0.001) nslicetemp++;
	472	}
	473	}
	474	if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
	475	}
	476
	477
4437a0d2	478	//test for kink tracks
	479	if(fExtraCuts->IsKinkMother(track)) hfetrack.SetIsKinkMother();
	480	else if(fExtraCuts->IsKinkDaughter(track)) hfetrack.SetIsKinkDaughter();
3513afb7	481
	482	// Double counted
	483	Int_t id(track->GetID());
	484	for(Int_t l=0; l < counterdc; l++){
	485	Int_t iTrack2 = arraytrack.At(l);
	486	if(iTrack2==id){
	487	hfetrack.SetDoubleCounted();
	488	break;
	489	}
	490	}
	491	// Add the id at this place
	492	arraytrack.AddAt(id,counterdc);
	493	counterdc++;
	494
	495	// PID
	496	hfetrack.SetTPCdEdx(track->GetTPCsignal());
	497	hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
	498	hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
	499	hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
4437a0d2	500	hfetrack.SetITSsigmaEl(pid->NumberOfSigmasITS(track, AliPID::kElectron));
3513afb7	501	// Eta correction
	502	copyTrack.~AliESDtrack();
	503	new(&copyTrack) AliESDtrack(*track);
	504	if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kESDanalysis);
	505	if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kESDanalysis);
	506	hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
	507	hfetrack.SetTPCdEdxCorrected(copyTrack.GetTPCsignal());
4437a0d2	508	if(isEMCAL){
4437a0d2	509	AliDebug(2, "Adding EMCAL PID information");
3513afb7	510	// EMCAL cluster
	511	Double_t emcalEnergyOverP = -1.,
	512	showershape[4] = {0.,0.,0.,0.};
	513	hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
	514	hfetrack.SetEMCALEoverP(emcalEnergyOverP);
	515	hfetrack.SetEMCALShowerShape(showershape);
	516	}
	517
	518	// Track finished, add NOW to the Event
	519	fHFEevent->AddTrack(&hfetrack);
	520	//printf("after\n");
	521	}
	522
	523	// Fill the debug tree
	524	//AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
	525	//AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
	526	fHFEtree->Fill();
	527
	528	fEventNumber++;
	529	PostData(1, fHFEtree);
	530	}
	531
	532	void AliHFEreducedEventCreatorESD::Terminate(Option_t *){
	533	//
	534	// Terminate
	535	//
	536	AliInfo("terminating...\n");
	537
	538	}
	539
	540	Bool_t AliHFEreducedEventCreatorESD::IsTOFmismatch(const AliVTrack const track, const AliPIDResponse const pid) const {
	541	//
	542	// Is TOF mismatch
	543	//
	544	Double_t probs[AliPID::kSPECIESC];
	545	AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
	546	return status == AliPIDResponse::kDetMismatch;
	547	}
	548