[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEreducedEventCreatorAOD.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 <TString.h>
#include <TTree.h>

#include "AliAnalysisManager.h"
#include "AliAnalysisUtils.h"
#include "AliAODEvent.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODPid.h"
#include "AliAODTrack.h"
#include "AliAODTracklets.h"
#include "AliAODVertex.h"
#include "AliCentrality.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 "AliPIDResponse.h"
#include "AliVEvent.h"
#include "AliVTrack.h"
#include "AliVCluster.h"
#include "AliVVZERO.h"
#include "AliVZDC.h"
#include "TTreeStream.h"

#include "AliHFEreducedEventCreatorAOD.h"

ClassImp(AliHFEreducedEventCreatorAOD)

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

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

AliHFEreducedEventCreatorAOD::~AliHFEreducedEventCreatorAOD(){
  //
  // Default destructor
  //
    if(fTPCpid) delete fTPCpid;
    if(fAnalysisUtils) delete fAnalysisUtils;
    if(fHFEevent) delete fHFEevent;
}

void AliHFEreducedEventCreatorAOD::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();
  fTrackCuts->SetAOD();
  // 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 AliHFEreducedEventCreatorAOD::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;
  }

  AliAODTrack copyTrack;

  // MC info
  Bool_t mcthere = kTRUE;
  AliAODEvent *aodE = dynamic_cast<AliAODEvent *>(fInputEvent);
  if(!aodE){
  //  printf("testd\n");
    AliError("No AOD Event");
    return;
  }
  fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
  if(!fAODMCHeader){ 
      mcthere = kFALSE;
  //   return;
  }
  fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
  if(!fAODArrayMCInfo){ 
      mcthere = kFALSE;
  //  return;
  }
  else {
    fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
    fTrackCuts->SetMCEvent(aodE);
  }

  // Set for track cuts
  if(mcthere){
    // Check Monte-Carlo events for AODs
    AliDebug(1, Form("Monte-Carlo Event: %p\n", fMCEvent));
    if(fMCEvent){
      AliDebug(1, Form("Available! Number of particles: %d\n", fMCEvent->GetNumberOfTracks()));
    }
  }
  fTrackCuts->SetRecEvent(fInputEvent);

  if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
    AliDebug(1, "Event rejected by the event cuts\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();

  // 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(fInputEvent->GetPrimaryVertex()->GetNContributors());
  fHFEevent->SetNContribVertex(ncontrib);

  // 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
  AliAODTracklets *tls = aodE->GetTracklets();
  if(tls) fHFEevent->SetSPDMultiplicity(tls->GetNumberOfTracklets());

  // Look for kink mother
  Int_t numberofvertices = aodE->GetNumberOfVertices();
  Double_t listofmotherkink[numberofvertices];
  Int_t numberofmotherkink = 0;
  //printf("Number of vertices %d\n",numberofvertices);
  for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
    AliAODVertex *aodvertex = aodE->GetVertex(ivertex);
    if(!aodvertex) continue;
    //printf("Type %d\n",aodvertex->GetType());
    if(aodvertex->GetType()==AliAODVertex::kKink) {
      //printf("Find one kink\n");
      AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
      if(!mother) continue;
      Int_t idmother = mother->GetID();
      listofmotherkink[numberofmotherkink] = idmother;
      //printf("ID %d\n",idmother);
      numberofmotherkink++;
    }
  }
  //printf("Number of kink mother in the events %d\n",numberofmotherkink);

  //
  // Loop on MC tracks only
  //
  AliAODMCParticle *mctrack(NULL);
  // Monte-Carlo info
  Int_t source(5);
  if(mcthere){
    for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
      mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(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->GetPdgCode());
      
      // 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 < fAODArrayMCInfo->GetEntriesFast()){
        AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
        if(mother) hfemcpart.SetMotherPdg(mother->GetPdgCode());
      }
      
      // 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->GetPdgCode()) == 11) source = 4;
      else source = 5;
      hfemcpart.SetSource(source);

      fHFEevent->AddMCParticle(&hfemcpart);
    }
  }
  
  //
  // Loop on reconstructed tracks
  //
  TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
  Int_t counterdc=0;
  
  AliAODTrack *track = 0x0;
  for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
    // Run track loop
    track = dynamic_cast<AliAODTrack *>(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->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P());

    // 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();
    // Filter
    for(Int_t k=0; k<20; k++) {
      Int_t u = 1<<k;     
      if((track->TestFilterBit(u))) {
	      hfetrack.SetFilterBit(k);
      } 
    }

    if(mcthere){
      // Fill Monte-Carlo Information
      Int_t label = TMath::Abs(track->GetLabel());
      if(label && label < fAODArrayMCInfo->GetEntriesFast())
        mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(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->GetPdgCode());
      
      // 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 < fAODArrayMCInfo->GetEntriesFast()){
        AliAODMCParticle *mother = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(motherlabel));
        if(mother) hfetrack.SetMCMotherPdg(mother->GetPdgCode());
      }
      
      // 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->GetPdgCode()) == 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->GetITSNcls()),
          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(track->GetTPCNcls());
    hfetrack.SetTRDnclusters(track->GetTRDncls());
    hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
    hfetrack.SetTPCcrossedRows(track->GetTPCNCrossedRows());
    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)
    hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
    hfetrack.SetTRDnslices(track->GetDetPid()->GetTRDnSlices());
    hfetrack.SetTRDchi2(track->GetTRDchi2());
    AliAODPid* aodpid= track->GetDetPid();
    Double_t* arraytrdsignals;
    arraytrdsignals=aodpid->GetTRDslices();
    Int_t nslicetemp=0;
    for(Int_t iplane = 0; iplane < 6; iplane++){
	    nslicetemp=0;
	    for(Int_t b=(iplane*8);b<((iplane*8)+8);b++){
	      if(track->GetTRDntrackletsPID()>0){
		      if(arraytrdsignals[b]>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;
    for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
      if(track->GetID() == listofmotherkink[kinkmother]) {
	      kinkmotherpass = 1;
	      continue;
      }
    }
    
    // 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->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.);
    hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
    hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
    hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
    // Eta correction
    copyTrack.~AliAODTrack();
    new(&copyTrack) AliAODTrack(*track);
    if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
    if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
    hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
    hfetrack.SetTPCdEdxCorrected(copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.);
    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 AliHFEreducedEventCreatorAOD::Terminate(Option_t *){
  //
  // Terminate
  //
  AliInfo("terminating...\n");

}

Bool_t AliHFEreducedEventCreatorAOD::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	#include <TArrayI.h>
	23	#include <TBits.h>
	24	#include <TFile.h>
	25	#include <TString.h>
	26	#include <TTree.h>
	27
	28	#include "AliAnalysisManager.h"
	29	#include "AliAnalysisUtils.h"
	30	#include "AliAODEvent.h"
	31	#include "AliAODMCHeader.h"
	32	#include "AliAODMCParticle.h"
	33	#include "AliAODPid.h"
	34	#include "AliAODTrack.h"
	35	#include "AliAODTracklets.h"
	36	#include "AliAODVertex.h"
	37	#include "AliCentrality.h"
	38	#include "AliInputEventHandler.h"
	39	#include "AliHFEcuts.h"
	40	#include "AliHFEextraCuts.h"
	41	#include "AliHFEpidTPC.h"
	42	#include "AliHFEreducedEvent.h"
	43	#include "AliHFEreducedTrack.h"
	44	#include "AliHFEreducedMCParticle.h"
	45	#include "AliHFEsignalCuts.h"
	46	#include "AliLog.h"
	47	#include "AliMCEvent.h"
	48	#include "AliPIDResponse.h"
	49	#include "AliVEvent.h"
	50	#include "AliVTrack.h"
	51	#include "AliVCluster.h"
	52	#include "AliVVZERO.h"
	53	#include "AliVZDC.h"
	54	#include "TTreeStream.h"
	55
	56	#include "AliHFEreducedEventCreatorAOD.h"
	57
	58	ClassImp(AliHFEreducedEventCreatorAOD)
	59
	60	AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD():
	61	AliAnalysisTaskSE(),
	62	fAODMCHeader(NULL),
	63	fAODArrayMCInfo(NULL),
	64	fHFEtree(NULL),
65	fAnalysisUtils(NULL),
66	fHFEevent(NULL),
67	fTrackCuts(NULL),
68	fExtraCuts(NULL),
69	fSignalCuts(NULL),
70	fTPCpid(NULL),
71	fEventNumber(0),
72	fNclustersTPC(70),
73	fNclustersTPCPID(0),
74	fNclustersITS(2),
75	fRemoveFirstEvent(kFALSE)
76	{
77	//
78	// Default constructor
79	//
80	}
81
82	AliHFEreducedEventCreatorAOD::AliHFEreducedEventCreatorAOD(const char *name):
83	AliAnalysisTaskSE(name),
84	fAODMCHeader(NULL),
85	fAODArrayMCInfo(NULL),
86	fHFEtree(NULL),
87	fAnalysisUtils(NULL),
88	fHFEevent(NULL),
89	fTrackCuts(NULL),
90	fExtraCuts(NULL),
91	fSignalCuts(NULL),
92	fTPCpid(NULL),
93	fEventNumber(0),
94	fNclustersTPC(70),
95	fNclustersTPCPID(0),
96	fNclustersITS(2),
97	fRemoveFirstEvent(kFALSE)
98	{
99	//
100	// Default constructor
101	//
102	fTPCpid = new AliHFEpidTPC("QAtpcPID");
103	fAnalysisUtils = new AliAnalysisUtils;
104	DefineOutput(1, TTree::Class());
105	}
106
107	AliHFEreducedEventCreatorAOD::~AliHFEreducedEventCreatorAOD(){
108	//
109	// Default destructor
110	//
111	if(fTPCpid) delete fTPCpid;
112	if(fAnalysisUtils) delete fAnalysisUtils;
113	if(fHFEevent) delete fHFEevent;
114	}
115
116	void AliHFEreducedEventCreatorAOD::UserCreateOutputObjects(){
117	//
118	// Create debug tree, signal cuts and track cuts
119	//
120
121	//printf("test\n");
122
123	// printf("testa\n");
124	fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
125	//printf("testb\n");
126
127	fTrackCuts = new AliHFEcuts("fTrackCuts", "Basic HFE track cuts");
128	fTrackCuts->CreateStandardCuts();
129	fTrackCuts->SetAOD();
130	// Track cuts
131	fTrackCuts->SetMinNClustersTPC(fNclustersTPC);
132	fTrackCuts->SetMinRatioTPCclusters(0);
133	fTrackCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable);
134	fTrackCuts->SetMinNClustersTPCPID(fNclustersTPCPID);
135	fTrackCuts->SetMinNClustersITS(fNclustersITS);
136	// Event cuts
137	fTrackCuts->SetUseMixedVertex(kTRUE);
138	fTrackCuts->SetVertexRange(10.);
139	//printf("testa\n");
140	fTrackCuts->Initialize();
141	//printf("testb\n");
142
143	fExtraCuts = new AliHFEextraCuts("hfeExtraCuts","HFE Extra Cuts");
144
145	fHFEevent = new AliHFEreducedEvent;
146	OpenFile(1);
147	fHFEtree = new TTree("HFEtree", "HFE event tree");
148	fHFEtree->Branch("HFEevent", "AliHFEreducedEvent", fHFEevent,128000,0);
149	PostData(1, fHFEtree);
150	}
151
152	void AliHFEreducedEventCreatorAOD::UserExec(Option_t *){
153	//
154	// User Exec: Fill debug Tree
155	//
156
157	// Get PID response
158	AliPIDResponse *pid = NULL;
159	AliInputEventHandler handler = dynamic_cast<AliInputEventHandler >(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
160	if(handler){
161	// printf("testb\n");
162	pid = handler->GetPIDResponse();
163	} else {
164	AliError("No Handler");
165	}
166	if(!pid){
167	// printf("testc\n");
168	AliError("No PID response");
169	return;
170	}
171	if(!fInputEvent) {
172	AliError("No Input event");
173	return;
174	}
175
176	if(fRemoveFirstEvent){
177	if(fAnalysisUtils->IsFirstEventInChunk(fInputEvent)) return;
178	}
179
180	AliAODTrack copyTrack;
181
182	// MC info
183	Bool_t mcthere = kTRUE;
184	AliAODEvent aodE = dynamic_cast<AliAODEvent >(fInputEvent);
185	if(!aodE){
186	// printf("testd\n");
187	AliError("No AOD Event");
188	return;
189	}
190	fAODMCHeader = dynamic_cast<AliAODMCHeader *>(fInputEvent->FindListObject(AliAODMCHeader::StdBranchName()));
191	if(!fAODMCHeader){
192	mcthere = kFALSE;
193	// return;
194	}
195	fAODArrayMCInfo = dynamic_cast<TClonesArray *>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
196	if(!fAODArrayMCInfo){
197	mcthere = kFALSE;
198	// return;
199	}
200	else {
201	fSignalCuts->SetMCAODInfo(fAODArrayMCInfo);
202	fTrackCuts->SetMCEvent(aodE);
203	}
204
205	// Set for track cuts
206	if(mcthere){
207	// Check Monte-Carlo events for AODs
208	AliDebug(1, Form("Monte-Carlo Event: %p\n", fMCEvent));
209	if(fMCEvent){
210	AliDebug(1, Form("Available! Number of particles: %d\n", fMCEvent->GetNumberOfTracks()));
211	}
212	}
213	fTrackCuts->SetRecEvent(fInputEvent);
214
215	if(!fTrackCuts->CheckEventCuts("fCutsEvRec", fInputEvent)){
216	AliDebug(1, "Event rejected by the event cuts\n");
217	return;
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();
236	if(trigger & AliVEvent::kCentral) fHFEevent->SetCentralTrigger();
237	if(trigger & AliVEvent::kSemiCentral) fHFEevent->SetCentralTrigger();
238	if(trigger & AliVEvent::kEMCEJE) fHFEevent->SetEMCALTrigger();
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]);
247	Int_t ncontrib(fInputEvent->GetPrimaryVertex()->GetNContributors());
248	fHFEevent->SetNContribVertex(ncontrib);
249
250	// Get centrality
251	AliCentrality *hicent = fInputEvent->GetCentrality();
252	fHFEevent->SetCentrality(
253	hicent->GetCentralityPercentile("V0M"),
254	hicent->GetCentralityPercentile("V0A"),
255	hicent->GetCentralityPercentile("V0C"),
256	hicent->GetCentralityPercentile("TKL"),
257	hicent->GetCentralityPercentile("TRK"),
258	hicent->GetCentralityPercentile("ZNA")
259	);
260
261	// Get VZERO Information
262	AliVVZERO *vzeroinfo = fInputEvent->GetVZEROData();
263	if(vzeroinfo) fHFEevent->SetV0Multiplicity(vzeroinfo->GetMTotV0A(), vzeroinfo->GetMTotV0C());
264
265	// Get ZDC Information
266	AliVZDC *zdcinfo = fInputEvent->GetZDCData();
267	if(zdcinfo) fHFEevent->SetZDCEnergy(zdcinfo->GetZNAEnergy(), zdcinfo->GetZNCEnergy(), zdcinfo->GetZPAEnergy(), zdcinfo->GetZPCEnergy());
268
269	// Set SPD multiplicity
270	AliAODTracklets *tls = aodE->GetTracklets();
271	if(tls) fHFEevent->SetSPDMultiplicity(tls->GetNumberOfTracklets());
272
273	// Look for kink mother
274	Int_t numberofvertices = aodE->GetNumberOfVertices();
275	Double_t listofmotherkink[numberofvertices];
276	Int_t numberofmotherkink = 0;
277	//printf("Number of vertices %d\n",numberofvertices);
278	for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
279	AliAODVertex *aodvertex = aodE->GetVertex(ivertex);
280	if(!aodvertex) continue;
281	//printf("Type %d\n",aodvertex->GetType());
282	if(aodvertex->GetType()==AliAODVertex::kKink) {
283	//printf("Find one kink\n");
284	AliAODTrack mother = (AliAODTrack ) aodvertex->GetParent();
285	if(!mother) continue;
286	Int_t idmother = mother->GetID();
287	listofmotherkink[numberofmotherkink] = idmother;
288	//printf("ID %d\n",idmother);
289	numberofmotherkink++;
290	}
291	}
292	//printf("Number of kink mother in the events %d\n",numberofmotherkink);
293
294	//
295	// Loop on MC tracks only
296	//
297	AliAODMCParticle *mctrack(NULL);
298	// Monte-Carlo info
299	Int_t source(5);
300	if(mcthere){
301	for(Int_t itrack = 0; itrack < fAODArrayMCInfo->GetEntriesFast(); itrack++) {
302	mctrack = (AliAODMCParticle *)(fAODArrayMCInfo->At(itrack));
303	if(!mctrack) continue;
304	AliHFEreducedMCParticle hfemcpart;
305	if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfemcpart.SetSignal();
306	// Kinematics
307	hfemcpart.SetSignedPt(mctrack->Pt(), mctrack->Charge() > 0.);
308	hfemcpart.SetP(mctrack->P());
309	hfemcpart.SetEta(mctrack->Eta());
310	hfemcpart.SetPhi(mctrack->Phi());
311	hfemcpart.SetPdg(mctrack->GetPdgCode());
312
313	// Get Production Vertex in radial direction
314	hfemcpart.SetProductionVertex(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
315
316	// Get Mother PDG code of the particle
317	Int_t motherlabel = TMath::Abs(mctrack->GetMother());
318	if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
319	AliAODMCParticle mother = dynamic_cast<AliAODMCParticle >(fAODArrayMCInfo->At(motherlabel));
320	if(mother) hfemcpart.SetMotherPdg(mother->GetPdgCode());
321	}
322
323	// derive source
324	source = 5;
325	if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
326	else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
327	else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
328	else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
329	else if(TMath::Abs(mctrack->GetPdgCode()) == 11) source = 4;
330	else source = 5;
331	hfemcpart.SetSource(source);
332
333	fHFEevent->AddMCParticle(&hfemcpart);
334	}
335	}
336
337	//
338	// Loop on reconstructed tracks
339	//
340	TArrayI arraytrack(fInputEvent->GetNumberOfTracks());
341	Int_t counterdc=0;
342
343	AliAODTrack *track = 0x0;
344	for(Int_t itrack = 0; itrack < fInputEvent->GetNumberOfTracks(); itrack++){
345	// Run track loop
346	track = dynamic_cast<AliAODTrack *>(fInputEvent->GetTrack(itrack));
347	if(!track) continue;
348	// Cut track (Only basic track cuts)
349	// printf("testv\n");
350	if(!fTrackCuts->CheckParticleCuts(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepRecKineITSTPC, track)) continue;
351	//
352	//printf("testu\n");
353
354	// Kinematics
355	AliHFEreducedTrack hfetrack;
356	hfetrack.SetSignedPt(track->Pt(), track->Charge() > 0);
357	hfetrack.SetP(track->P());
358	hfetrack.SetEta(track->Eta());
359	hfetrack.SetPhi(track->Phi());
360	hfetrack.SetTPCmomentum(track->GetDetPid() ? track->GetDetPid()->GetTPCmomentum() : track->P());
361
362	// Track ID
363	hfetrack.SetTrackID(track->GetID());
364
365	// status
366	ULong_t status = track->GetStatus();
367	if((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) hfetrack.SetITSrefit();
368	if((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) hfetrack.SetTPCrefit();
369	if((status & AliVTrack::kTOFpid) == AliVTrack::kTOFpid) hfetrack.SetTOFpid();
370	//if((status & AliVTrack::kTOFmismatch) == AliVTrack::kTOFmismatch) hfetrack.SetTOFmismatch();
371	if(IsTOFmismatch(track, pid)) hfetrack.SetTOFmismatch(); // New version suggested by Pietro Antonioli
372	if(track->IsEMCAL()) hfetrack.SetEMCALpid();
373	// Filter
374	for(Int_t k=0; k<20; k++) {
375	Int_t u = 1<<k;
376	if((track->TestFilterBit(u))) {
377	hfetrack.SetFilterBit(k);
378	}
379	}
380
381	if(mcthere){
382	// Fill Monte-Carlo Information
383	Int_t label = TMath::Abs(track->GetLabel());
384	if(label && label < fAODArrayMCInfo->GetEntriesFast())
385	mctrack = dynamic_cast<AliAODMCParticle *>(fAODArrayMCInfo->At(label));
386	if(!mctrack) continue;
387	if(fTrackCuts->CheckParticleCuts(static_cast<UInt_t>(AliHFEcuts::kStepMCGenerated), mctrack)) hfetrack.SetMCSignal();
388	// Kinematics
389	hfetrack.SetMCSignedPt(mctrack->Pt(),mctrack->Charge() > 0.);
390	hfetrack.SetMCP(mctrack->P());
391	hfetrack.SetMCEta(mctrack->Eta());
392	hfetrack.SetMCPhi(mctrack->Phi());
393	hfetrack.SetMCPDG(mctrack->GetPdgCode());
394
395	// Get Production Vertex in radial direction
396	hfetrack.SetMCProdVtx(mctrack->Xv(),mctrack->Yv(),mctrack->Zv());
397
398	// Get Mother PDG code of the particle
399	Int_t motherlabel = TMath::Abs(mctrack->GetMother());
400	if(motherlabel >= 0 && motherlabel < fAODArrayMCInfo->GetEntriesFast()){
401	AliAODMCParticle mother = dynamic_cast<AliAODMCParticle >(fAODArrayMCInfo->At(motherlabel));
402	if(mother) hfetrack.SetMCMotherPdg(mother->GetPdgCode());
403	}
404
405	// derive source
406	source = 5;
407	if(fSignalCuts->IsCharmElectron(mctrack)) source = 0;
408	else if(fSignalCuts->IsBeautyElectron(mctrack)) source = 1;
409	else if(fSignalCuts->IsGammaElectron(mctrack)) source = 2;
410	else if(fSignalCuts->IsNonHFElectron(mctrack)) source = 3;
411	else if(TMath::Abs(mctrack->GetPdgCode()) == 11) source = 4;
412	else source = 5;
413	hfetrack.SetMCSource(source);
414	}
415
416	// HFE DCA
417	Float_t dcaxy = -999.,
418	dcaz = -999.;
419	fExtraCuts->GetImpactParameters((AliVTrack *)track,dcaxy,dcaz);
420	hfetrack.SetDCA(dcaxy, dcaz);
421
422	// Different number of clusters definitions
423	Int_t nclustersITS(track->GetITSNcls()),
424	nclustersTPC(track->GetTPCNcls()),
425	nclustersTPCall(track->GetTPCClusterMap().CountBits()),
426	nclustersTPCshared(0);
427	UChar_t nfindableTPC = track->GetTPCNclsF();
428	const TBits &sharedTPC = track->GetTPCSharedMap();
429	for(Int_t ibit = 0; ibit < 160; ibit++) if(sharedTPC.TestBitNumber(ibit)) nclustersTPCshared++;
430	hfetrack.SetITSnclusters(nclustersITS);
431	hfetrack.SetTPCnclusters(track->GetTPCNcls());
432	hfetrack.SetTRDnclusters(track->GetTRDncls());
433	hfetrack.SetTPCnclustersPID(track->GetTPCsignalN());
434	hfetrack.SetTPCcrossedRows(track->GetTPCNCrossedRows());
435	hfetrack.SetTPCnclustersAll(nclustersTPCall);
436	hfetrack.SetTPCsharedClusters(nclustersTPCshared);
437	hfetrack.SetTPCclusterRatio(nfindableTPC ? static_cast<Float_t>(nclustersTPC)/static_cast<Float_t>(nfindableTPC) : 0);
438	hfetrack.SetTPCclusterRatioAll(nfindableTPC ? static_cast<Float_t>(nclustersTPCall)/static_cast<Float_t>(nfindableTPC) : 0);
439	UChar_t itsPixel = track->GetITSClusterMap();
440	for(int ily = 0; ily < 6; ily++)
441	if(TESTBIT(itsPixel, ily)) hfetrack.SetITScluster(ily);
442
443	// TRD related quantities (Yvonne)
444	hfetrack.SetTRDntrackletsPID(track->GetTRDntrackletsPID());
445	hfetrack.SetTRDnslices(track->GetDetPid()->GetTRDnSlices());
446	hfetrack.SetTRDchi2(track->GetTRDchi2());
447	AliAODPid* aodpid= track->GetDetPid();
448	Double_t* arraytrdsignals;
449	arraytrdsignals=aodpid->GetTRDslices();
450	Int_t nslicetemp=0;
451	for(Int_t iplane = 0; iplane < 6; iplane++){
452	nslicetemp=0;
453	for(Int_t b=(iplane8);b<((iplane8)+8);b++){
454	if(track->GetTRDntrackletsPID()>0){
455	if(arraytrdsignals[b]>0.001) nslicetemp++;
456	}
457	}
458	if(nslicetemp > 0) hfetrack.SetTRDstatus(iplane);
459	}
460
461
462	// Kink
463	Int_t kink = 0;
464	if(fExtraCuts->IsKinkDaughter(track)) kink = 1;
465
466	// kink mother
467	Int_t kinkmotherpass = 0;
468	for(Int_t kinkmother = 0; kinkmother < numberofmotherkink; kinkmother++) {
469	if(track->GetID() == listofmotherkink[kinkmother]) {
470	kinkmotherpass = 1;
471	continue;
472	}
473	}
474
475	// Double counted
476	Int_t id(track->GetID());
477	for(Int_t l=0; l < counterdc; l++){
478	Int_t iTrack2 = arraytrack.At(l);
479	if(iTrack2==id){
480	hfetrack.SetDoubleCounted();
481	break;
482	}
483	}
484	// Add the id at this place
485	arraytrack.AddAt(id,counterdc);
486	counterdc++;
487
488	// PID
489	hfetrack.SetTPCdEdx(track->GetDetPid() ? track->GetDetPid()->GetTPCsignal() : 0.);
490	hfetrack.SetTPCsigmaEl(pid->NumberOfSigmasTPC(track, AliPID::kElectron));
491	hfetrack.SetTOFsigmaEl(pid->NumberOfSigmasTOF(track, AliPID::kElectron));
492	hfetrack.SetTOFmismatchProbability(pid->GetTOFMismatchProbability(track));
493	// Eta correction
494	copyTrack.~AliAODTrack();
495	new(&copyTrack) AliAODTrack(*track);
496	if(fTPCpid->HasCentralityCorrection()) fTPCpid->ApplyCentralityCorrection(&copyTrack, static_cast<Double_t>(ncontrib),AliHFEpidObject::kAODanalysis);
497	if(fTPCpid->HasEtaCorrection()) fTPCpid->ApplyEtaCorrection(&copyTrack, AliHFEpidObject::kAODanalysis);
498	hfetrack.SetTPCsigmaElCorrected(pid->NumberOfSigmasTPC(&copyTrack, AliPID::kElectron));
499	hfetrack.SetTPCdEdxCorrected(copyTrack.GetDetPid() ? copyTrack.GetDetPid()->GetTPCsignal() : 0.);
500	if(track->IsEMCAL()){
501	// EMCAL cluster
502	Double_t emcalEnergyOverP = -1.,
503	showershape[4] = {0.,0.,0.,0.};
504	hfetrack.SetEMCALSigmaEl(pid->NumberOfSigmasEMCAL(track, AliPID::kElectron, emcalEnergyOverP, &showershape[0]));
505	hfetrack.SetEMCALEoverP(emcalEnergyOverP);
506	hfetrack.SetEMCALShowerShape(showershape);
507	}
508
509	// Track finished, add NOW to the Event
510	fHFEevent->AddTrack(&hfetrack);
511	//printf("after\n");
512	}
513
514	// Fill the debug tree
515	AliInfo(Form("Number of tracks: %d\n", fHFEevent->GetNumberOfTracks()));
516	AliInfo(Form("Number of MC particles: %d\n", fHFEevent->GetNumberOfMCParticles()));
517	fHFEtree->Fill();
518
519	fEventNumber++;
520	PostData(1, fHFEtree);
521	}
522
523	void AliHFEreducedEventCreatorAOD::Terminate(Option_t *){
524	//
525	// Terminate
526	//
527	AliInfo("terminating...\n");
528
529	}
530
531	Bool_t AliHFEreducedEventCreatorAOD::IsTOFmismatch(const AliVTrack const track, const AliPIDResponse const pid) const {
532	//
533	// is TOF mismatch
534	//
535	Double_t probs[AliPID::kSPECIESC];
536	AliPIDResponse::EDetPidStatus status = pid->ComputeTOFProbability(track, AliPID::kSPECIESC, probs);
537	return status == AliPIDResponse::kDetMismatch;
538	}
539