[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->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::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")
  );
  
  // 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)) 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
    if(track->IsEMCAL()) hfetrack.SetEMCALpid();
    // no filter bits available for ESDs

    if(mcthere){
      // Fill Monte-Carlo Information
      Int_t label = TMath::Abs(track->GetLabel());
      if(label && label < fMCEvent->GetNumberOfTracks())
        mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
      if(!mctrack) continue;
      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); 
    }

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

    // 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.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);
    }


    // Kink
    Int_t kink = 0;
    if(fExtraCuts->IsKinkDaughter(track)) kink = 1;

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