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