[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDinfoGen.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.                  *
**************************************************************************/

/* $Id: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//
//  Tender wagon for TRD performance/calibration train
//
// In this wagon the information from
//   - ESD
//   - Friends [if available]
//   - MC [if available]
// are grouped into AliTRDtrackInfo objects and fed to worker tasks
//
//  Authors:
//    Markus Fasel <M.Fasel@gsi.de>
//    Alexandru Bercuci <A.Bercuci@gsi.de>
//
////////////////////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TObjArray.h>
#include <TObject.h>
#include <TString.h>
#include <TH1F.h>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <TChain.h>
#include <TParticle.h>

#include "AliLog.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDInputHandler.h"
#include "AliMCEventHandler.h"

#include "AliESDfriend.h"
#include "AliESDfriendTrack.h"
#include "AliESDHeader.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliMCParticle.h"
#include "AliPID.h"
#include "AliStack.h"
#include "AliTRDtrackV1.h"
#include "AliTrackReference.h"
#include "AliTRDgeometry.h"
#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "TTreeStream.h"

#include <cstdio>
#include <climits>
#include <cstring>
#include <iostream>

#include "AliTRDinfoGen.h"
#include "info/AliTRDtrackInfo.h"
#include "info/AliTRDeventInfo.h"
#include "info/AliTRDv0Info.h"
#include "info/AliTRDeventCuts.h"
#include "macros/AliTRDperformanceTrain.h"

ClassImp(AliTRDinfoGen)

const Float_t AliTRDinfoGen::fgkTPC = 290.;
const Float_t AliTRDinfoGen::fgkTOF = 365.;

const Float_t AliTRDinfoGen::fgkEvVertexZ = 15.;
const Int_t   AliTRDinfoGen::fgkEvVertexN = 1;

const Float_t AliTRDinfoGen::fgkTrkDCAxy  = 40.;
const Float_t AliTRDinfoGen::fgkTrkDCAz   = 15.;
const Int_t   AliTRDinfoGen::fgkNclTPC    = 10;
const Float_t AliTRDinfoGen::fgkPt        = 0.2;
const Float_t AliTRDinfoGen::fgkEta       = 0.9;

//____________________________________________________________________
AliTRDinfoGen::AliTRDinfoGen()
  :AliAnalysisTaskSE()
  ,fEvTrigger(NULL)
  ,fESDev(NULL)
  ,fMCev(NULL)
  ,fEventCut(NULL)
  ,fTrackCut(NULL)
  ,fTrackInfo(NULL)
  ,fEventInfo(NULL)
  ,fV0Info(NULL)
  ,fTracksBarrel(NULL)
  ,fTracksSA(NULL)
  ,fTracksKink(NULL)
  ,fV0List(NULL)
  ,fDebugStream(NULL)
{
  //
  // Default constructor
  //
  SetNameTitle("infoGen", "MC-REC TRD-track list generator");
}

//____________________________________________________________________
AliTRDinfoGen::AliTRDinfoGen(char* name)
  :AliAnalysisTaskSE(name)
  ,fEvTrigger(NULL)
  ,fESDev(NULL)
  ,fMCev(NULL)
  ,fEventCut(NULL)
  ,fTrackCut(NULL)
  ,fTrackInfo(NULL)
  ,fEventInfo(NULL)
  ,fV0Info(NULL)
  ,fTracksBarrel(NULL)
  ,fTracksSA(NULL)
  ,fTracksKink(NULL)
  ,fV0List(NULL)
  ,fDebugStream(NULL)
{
  //
  // Default constructor
  //
  SetTitle("MC-REC TRD-track list generator");
  DefineOutput(kTracksBarrel, TObjArray::Class());
  DefineOutput(kTracksSA, TObjArray::Class());
  DefineOutput(kTracksKink, TObjArray::Class());
  DefineOutput(kEventInfo, AliTRDeventInfo::Class());
  DefineOutput(kV0List, TObjArray::Class());
}

//____________________________________________________________________
AliTRDinfoGen::~AliTRDinfoGen()
{
// Destructor
  
  if(fEvTrigger) delete fEvTrigger;
  if(fTrackCut) delete fTrackCut;
  if(fEventCut) delete fEventCut;
  if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;
  if(fEventInfo) delete fEventInfo; fEventInfo = NULL;
  if(fV0Info) delete fV0Info; fV0Info = NULL;
  if(fTracksBarrel){ 
    fTracksBarrel->Delete(); delete fTracksBarrel;
    fTracksBarrel = NULL;
  }
  if(fTracksSA){ 
    fTracksSA->Delete(); delete fTracksSA;
    fTracksSA = NULL;
  }
  if(fTracksKink){ 
    fTracksKink->Delete(); delete fTracksKink;
    fTracksKink = NULL;
  }
  if(fV0List){ 
    fV0List->Delete(); delete fV0List;
    fV0List = NULL;
  }
}

//____________________________________________________________________
void AliTRDinfoGen::UserCreateOutputObjects()
{	
  //
  // Create Output Containers (TObjectArray containing 1D histograms)
  //
 
  fTrackInfo = new AliTRDtrackInfo();
  fEventInfo = new AliTRDeventInfo();
  fV0Info    = new AliTRDv0Info();
  fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);
  fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);
  fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);
  fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);
}

//____________________________________________________________________
void AliTRDinfoGen::UserExec(Option_t *){
  //
  // Run the Analysis
  //

  fESDev = dynamic_cast<AliESDEvent*>(InputEvent());
  fMCev = MCEvent();

  if(!fESDev){
    AliError("Failed retrieving ESD event");
    return;
  }

  // event selection : trigger cut
  if(UseLocalEvSelection() && fEvTrigger){ 
    Bool_t kTRIGGERED(kFALSE);
    const TObjArray *trig = fEvTrigger->Tokenize(" ");
    for(Int_t itrig=trig->GetEntriesFast(); itrig--;){
      const Char_t *trigClass(((TObjString*)(*trig)[itrig])->GetName());
      if(fESDev->IsTriggerClassFired(trigClass)) {
        AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));
        kTRIGGERED = kTRUE;
        break; 
      }
    }
    if(!kTRIGGERED) return;

    // select only physical events
    if(fESDev->GetEventType() != 7) return;
  }

  // if the required trigger is a collision trigger then apply event vertex cut
  if(UseLocalEvSelection() && IsCollision()){
    const AliESDVertex *vertex = fESDev->GetPrimaryVertex();
    if(TMath::Abs(vertex->GetZv())<1.e-10 || 
       TMath::Abs(vertex->GetZv())>fgkEvVertexZ || 
       vertex->GetNContributors()<fgkEvVertexN) {
      //PostData(0, fOutStorage);
      return;
    }
  }

  if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;

  if(!fESDfriend){
    AliError("Failed retrieving ESD friend event");
    return;
  }
  if(HasMCdata() && !fMCev){
    AliError("Failed retrieving MC event");
    return;
  }

  fTracksBarrel->Delete();
  fTracksSA->Delete();
  fTracksKink->Delete();
  fV0List->Delete();
  fEventInfo->Delete("");
  new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
  
  Bool_t *trackMap(NULL);
  AliStack * mStack(NULL);
  Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
  if(HasMCdata()){
    mStack = fMCev->Stack();
    if(!mStack){
      AliError("Failed retrieving MC Stack");
      return;
    }
    trackMap = new Bool_t[nTracksMC];
    memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
  }
  
  Double32_t dedx[100]; Int_t nSlices(0);
  Int_t nTRDout(0), nTRDin(0), nTPC(0), nclsTrklt;
  AliDebug(2, Form("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC));
  AliESDtrack *esdTrack = NULL;
  AliESDfriendTrack *esdFriendTrack = NULL;
  TObject *calObject = NULL;
  AliTRDtrackV1 *track = NULL;
  AliTRDseedV1 *tracklet = NULL;
  AliTRDcluster *cl = NULL;
  for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
    esdTrack = fESDev->GetTrack(itrk);
    AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));

    if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
    if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
    if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;

    // look at external track param
    const AliExternalTrackParam *op = esdTrack->GetOuterParam();
    Double_t xyz[3];
    if(op){
      op->GetXYZ(xyz);
      op->Global2LocalPosition(xyz, op->GetAlpha());
      AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));
    }

    // read MC info
    Int_t fPdg = -1;
    Int_t label = -1; UInt_t alab=UINT_MAX;
    if(HasMCdata()){
      label = esdTrack->GetLabel(); 
      alab = TMath::Abs(label);
      // register the track
      if(alab < UInt_t(nTracksMC)){ 
        trackMap[alab] = kTRUE; 
      } else { 
        AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
        continue; 
      }
      AliMCParticle *mcParticle = NULL; 
      if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
        AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
        continue;
      }
      fPdg = mcParticle->Particle()->GetPdgCode();
      Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
      Int_t iref = 0; AliTrackReference *ref = NULL; 
      while(iref<nRefs){
        ref = mcParticle->GetTrackReference(iref);
        if(ref->LocalX() > fgkTPC) break;
        iref++;
      }

      new(fTrackInfo) AliTRDtrackInfo();
      fTrackInfo->SetPDG(fPdg);
      fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
      Int_t jref = iref;//, kref = 0;
      while(jref<nRefs){
        ref = mcParticle->GetTrackReference(jref);
        if(ref->LocalX() > fgkTOF) break;
        AliDebug(4, Form("  trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));
        fTrackInfo->AddTrackRef(ref);
        jref++;
      }
      AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));
    } else {
      new (fTrackInfo) AliTRDtrackInfo();
      fTrackInfo->SetPDG(fPdg);
    }

    // copy some relevant info to TRD track info
    fTrackInfo->SetStatus(esdTrack->GetStatus());
    fTrackInfo->SetTrackId(esdTrack->GetID());
    Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
    fTrackInfo->SetESDpid(p);
    fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
    if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();
    memset(dedx, 0, 100*sizeof(Double32_t));
    Int_t in(0);
    for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)
      for(Int_t is=0; is<nSlices; is++) 
        dedx[in++]=esdTrack->GetTRDslice(il, is);
    for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);
    fTrackInfo->SetSlices(in, dedx);
    fTrackInfo->SetLabel(label);
    fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
    // some other Informations which we may wish to store in order to find problematic cases
    fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
    fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
    nclsTrklt = 0;
  

    // read REC info
    esdFriendTrack = fESDfriend->GetTrack(itrk);
    if(esdFriendTrack){
      Int_t icalib = 0;
      while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
        if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
        if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
        AliDebug(4, Form("TRD track OK"));
        // Set the clusters to unused
        for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
          if(!(tracklet = track->GetTracklet(ipl))) continue;
          tracklet->ResetClusterIter();
          while((cl = tracklet->NextCluster())) cl->Use(0);
        }
        fTrackInfo->SetTrack(track);
        break;
      }
      AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));
    } else AliDebug(3, "No ESD friends");
    if(op) fTrackInfo->SetOuterParam(op);

    if(DebugLevel() >= 1){
      AliTRDtrackInfo info(*fTrackInfo);
      (*DebugStream()) << "trackInfo"
      << "TrackInfo.=" << &info
      << "\n";
      info.Delete("");
    }

    ULong_t status(esdTrack->GetStatus());
    if((status&AliESDtrack::kTPCout)){
      if(!esdTrack->GetKinkIndex(0)){ // Barrel  Track Selection
        Bool_t selected(kTRUE);
        if(UseLocalTrkSelection()){
          if(esdTrack->Pt() < fgkPt) selected = kFALSE;
          if(TMath::Abs(esdTrack->Eta()) > fgkEta) selected = kFALSE;
          if(esdTrack->GetTPCNcls() < fgkNclTPC) selected = kFALSE;
          Float_t par[2], cov[3];
          esdTrack->GetImpactParameters(par, cov);
          if(IsCollision()){ // cuts on DCA
            if(TMath::Abs(par[0]) > fgkTrkDCAxy) selected = kFALSE;
            if(TMath::Abs(par[1]) > fgkTrkDCAz) selected = kFALSE;
          }
        }
        if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;
        if(selected) fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
      } else fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
    } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)) fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
    fTrackInfo->Delete("");
  }

  AliDebug(2, Form("%3d Tracks: TPCout[%d] TRDin[%d] TRDout[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRDin, nTRDout));

//   AliESDv0 *v0 = NULL;
//   for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
//     if(!(v0 = fESD->GetV0(iv0))) continue;
//     fV0List->Add(new AliTRDv0Info(v0));
//   }

  // Insert also MC tracks which are passing TRD where the track is not reconstructed
  if(HasMCdata()){
    AliDebug(10, "Output of the MC track map:");
    for(Int_t itk = 0; itk < nTracksMC;  itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));
  
    for(Int_t itk = 0; itk < nTracksMC; itk++){
      if(trackMap[itk]) continue;
      AliMCParticle *mcParticle =  (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));
      Int_t fPdg = mcParticle->Particle()->GetPdgCode();
      Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
      Int_t iref = 0; AliTrackReference *ref = NULL; 
      Int_t nRefsTRD = 0;
      new(fTrackInfo) AliTRDtrackInfo();
      fTrackInfo->SetPDG(fPdg);
      while(iref<nRefs){
        Bool_t kIN(kFALSE);
        ref = mcParticle->GetTrackReference(iref);
        if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTOF){
          fTrackInfo->AddTrackRef(ref);
          nRefsTRD++;kIN=kTRUE;
        }
        AliDebug(4, Form("  trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
        iref++;
      }
      if(!nRefsTRD){
        // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the 
        // analysis job
        fTrackInfo->Delete("");
        continue;
      }
      fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
      fTrackInfo->SetLabel(itk);
      if(DebugLevel() >= 1){
        AliTRDtrackInfo info(*fTrackInfo);
        (*DebugStream()) << "trackInfo"
        << "TrackInfo.=" << &info
        << "\n";
        info.Delete("");
      }
      AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
      fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
      fTrackInfo->Delete("");
    }
    delete[] trackMap;
  }
  PostData(kTracksBarrel, fTracksBarrel);
  PostData(kTracksSA, fTracksSA);
  PostData(kTracksKink, fTracksKink);
  PostData(kEventInfo, fEventInfo);
  PostData(kV0List, fV0List);
}

//____________________________________________________________________
void AliTRDinfoGen::SetTrigger(const Char_t *trigger)
{
  if(!fEvTrigger) fEvTrigger = new TString(trigger);
  else (*fEvTrigger) = trigger;
}

//____________________________________________________________________
TTreeSRedirector* AliTRDinfoGen::DebugStream()
{
  if(!fDebugStream){
    TDirectory *savedir = gDirectory;
    fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");
    savedir->cd();
  }
  return fDebugStream;
}
Commit	Line	Data
1ee39b3a	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	/* $Id: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	//
	20	// Tender wagon for TRD performance/calibration train
	21	//
	22	// In this wagon the information from
	23	// - ESD
	24	// - Friends [if available]
	25	// - MC [if available]
	26	// are grouped into AliTRDtrackInfo objects and fed to worker tasks
	27	//
	28	// Authors:
	29	// Markus Fasel <M.Fasel@gsi.de>
	30	// Alexandru Bercuci <A.Bercuci@gsi.de>
	31	//
	32	////////////////////////////////////////////////////////////////////////////
	33
	34	#include <TClonesArray.h>
	35	#include <TObjArray.h>
	36	#include <TObject.h>
705f8b0a	37	#include <TString.h>
1ee39b3a	38	#include <TH1F.h>
	39	#include <TFile.h>
	40	#include <TTree.h>
	41	#include <TROOT.h>
	42	#include <TChain.h>
	43	#include <TParticle.h>
	44
	45	#include "AliLog.h"
	46	#include "AliAnalysisManager.h"
	47	#include "AliESDEvent.h"
	48	#include "AliMCEvent.h"
	49	#include "AliESDInputHandler.h"
	50	#include "AliMCEventHandler.h"
	51
	52	#include "AliESDfriend.h"
	53	#include "AliESDfriendTrack.h"
	54	#include "AliESDHeader.h"
	55	#include "AliESDtrack.h"
07ece31c	56	#include "AliESDtrackCuts.h"
1ee39b3a	57	#include "AliMCParticle.h"
	58	#include "AliPID.h"
	59	#include "AliStack.h"
	60	#include "AliTRDtrackV1.h"
	61	#include "AliTrackReference.h"
	62	#include "AliTRDgeometry.h"
	63	#include "AliTRDcluster.h"
	64	#include "AliTRDseedV1.h"
	65	#include "TTreeStream.h"
	66
	67	#include <cstdio>
	68	#include <climits>
	69	#include <cstring>
	70	#include <iostream>
	71
	72	#include "AliTRDinfoGen.h"
	73	#include "info/AliTRDtrackInfo.h"
	74	#include "info/AliTRDeventInfo.h"
	75	#include "info/AliTRDv0Info.h"
07ece31c	76	#include "info/AliTRDeventCuts.h"
705f8b0a	77	#include "macros/AliTRDperformanceTrain.h"
1ee39b3a	78
	79	ClassImp(AliTRDinfoGen)
	80
	81	const Float_t AliTRDinfoGen::fgkTPC = 290.;
	82	const Float_t AliTRDinfoGen::fgkTOF = 365.;
	83
3f19c2cd	84	const Float_t AliTRDinfoGen::fgkEvVertexZ = 15.;
	85	const Int_t AliTRDinfoGen::fgkEvVertexN = 1;
	86
	87	const Float_t AliTRDinfoGen::fgkTrkDCAxy = 40.;
	88	const Float_t AliTRDinfoGen::fgkTrkDCAz = 15.;
	89	const Int_t AliTRDinfoGen::fgkNclTPC = 10;
	90	const Float_t AliTRDinfoGen::fgkPt = 0.2;
	91	const Float_t AliTRDinfoGen::fgkEta = 0.9;
	92
1ee39b3a	93	//____________________________________________________________________
db99a57a	94	AliTRDinfoGen::AliTRDinfoGen()
705f8b0a	95	:AliAnalysisTaskSE()
705f8b0a	96	,fEvTrigger(NULL)
770382d9	97	,fESDev(NULL)
770382d9	98	,fMCev(NULL)
705f8b0a	99	,fEventCut(NULL)
705f8b0a	100	,fTrackCut(NULL)
770382d9	101	,fTrackInfo(NULL)
770382d9	102	,fEventInfo(NULL)
770382d9	103	,fV0Info(NULL)
705f8b0a	104	,fTracksBarrel(NULL)
	105	,fTracksSA(NULL)
	106	,fTracksKink(NULL)
	107	,fV0List(NULL)
	108	,fDebugStream(NULL)
1ee39b3a	109	{
	110	//
	111	// Default constructor
	112	//
f2e89a4c	113	SetNameTitle("infoGen", "MC-REC TRD-track list generator");
f8f46e4d	114	}
1ee39b3a	115
07ece31c	116	//____________________________________________________________________
db99a57a	117	AliTRDinfoGen::AliTRDinfoGen(char* name)
705f8b0a	118	:AliAnalysisTaskSE(name)
705f8b0a	119	,fEvTrigger(NULL)
f8f46e4d	120	,fESDev(NULL)
f8f46e4d	121	,fMCev(NULL)
705f8b0a	122	,fEventCut(NULL)
705f8b0a	123	,fTrackCut(NULL)
f8f46e4d	124	,fTrackInfo(NULL)
f8f46e4d	125	,fEventInfo(NULL)
f8f46e4d	126	,fV0Info(NULL)
705f8b0a	127	,fTracksBarrel(NULL)
	128	,fTracksSA(NULL)
	129	,fTracksKink(NULL)
	130	,fV0List(NULL)
	131	,fDebugStream(NULL)
f8f46e4d	132	{
	133	//
	134	// Default constructor
	135	//
705f8b0a	136	SetTitle("MC-REC TRD-track list generator");
	137	DefineOutput(kTracksBarrel, TObjArray::Class());
	138	DefineOutput(kTracksSA, TObjArray::Class());
	139	DefineOutput(kTracksKink, TObjArray::Class());
	140	DefineOutput(kEventInfo, AliTRDeventInfo::Class());
	141	DefineOutput(kV0List, TObjArray::Class());
1ee39b3a	142	}
	143
	144	//____________________________________________________________________
	145	AliTRDinfoGen::~AliTRDinfoGen()
	146	{
	147	// Destructor
705f8b0a	148
705f8b0a	149	if(fEvTrigger) delete fEvTrigger;
07ece31c	150	if(fTrackCut) delete fTrackCut;
07ece31c	151	if(fEventCut) delete fEventCut;
770382d9	152	if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;
770382d9	153	if(fEventInfo) delete fEventInfo; fEventInfo = NULL;
770382d9	154	if(fV0Info) delete fV0Info; fV0Info = NULL;
705f8b0a	155	if(fTracksBarrel){
	156	fTracksBarrel->Delete(); delete fTracksBarrel;
	157	fTracksBarrel = NULL;
	158	}
	159	if(fTracksSA){
	160	fTracksSA->Delete(); delete fTracksSA;
	161	fTracksSA = NULL;
1ee39b3a	162	}
705f8b0a	163	if(fTracksKink){
	164	fTracksKink->Delete(); delete fTracksKink;
	165	fTracksKink = NULL;
	166	}
	167	if(fV0List){
	168	fV0List->Delete(); delete fV0List;
	169	fV0List = NULL;
1ee39b3a	170	}
	171	}
	172
	173	//____________________________________________________________________
f8f46e4d	174	void AliTRDinfoGen::UserCreateOutputObjects()
1ee39b3a	175	{
	176	//
	177	// Create Output Containers (TObjectArray containing 1D histograms)
	178	//
b4414720	179
1ee39b3a	180	fTrackInfo = new AliTRDtrackInfo();
	181	fEventInfo = new AliTRDeventInfo();
	182	fV0Info = new AliTRDv0Info();
705f8b0a	183	fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);
	184	fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);
	185	fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);
	186	fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);
1ee39b3a	187	}
	188
	189	//____________________________________________________________________
f8f46e4d	190	void AliTRDinfoGen::UserExec(Option_t *){
1ee39b3a	191	//
	192	// Run the Analysis
	193	//
f2e89a4c	194
f8f46e4d	195	fESDev = dynamic_cast<AliESDEvent*>(InputEvent());
	196	fMCev = MCEvent();
	197
1ee39b3a	198	if(!fESDev){
	199	AliError("Failed retrieving ESD event");
	200	return;
	201	}
3f19c2cd	202
3f19c2cd	203	// event selection : trigger cut
705f8b0a	204	if(UseLocalEvSelection() && fEvTrigger){
3f19c2cd	205	Bool_t kTRIGGERED(kFALSE);
705f8b0a	206	const TObjArray *trig = fEvTrigger->Tokenize(" ");
3f19c2cd	207	for(Int_t itrig=trig->GetEntriesFast(); itrig--;){
	208	const Char_t trigClass(((TObjString)(*trig)[itrig])->GetName());
	209	if(fESDev->IsTriggerClassFired(trigClass)) {
	210	AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));
	211	kTRIGGERED = kTRUE;
	212	break;
	213	}
	214	}
	215	if(!kTRIGGERED) return;
	216
	217	// select only physical events
	218	if(fESDev->GetEventType() != 7) return;
	219	}
	220
	221	// if the required trigger is a collision trigger then apply event vertex cut
	222	if(UseLocalEvSelection() && IsCollision()){
	223	const AliESDVertex *vertex = fESDev->GetPrimaryVertex();
	224	if(TMath::Abs(vertex->GetZv())<1.e-10 \|\|
	225	TMath::Abs(vertex->GetZv())>fgkEvVertexZ \|\|
	226	vertex->GetNContributors()<fgkEvVertexN) {
	227	//PostData(0, fOutStorage);
	228	return;
	229	}
	230	}
	231
07ece31c	232	if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;
07ece31c	233
1ee39b3a	234	if(!fESDfriend){
	235	AliError("Failed retrieving ESD friend event");
	236	return;
	237	}
	238	if(HasMCdata() && !fMCev){
	239	AliError("Failed retrieving MC event");
	240	return;
	241	}
f8f46e4d	242
705f8b0a	243	fTracksBarrel->Delete();
	244	fTracksSA->Delete();
	245	fTracksKink->Delete();
	246	fV0List->Delete();
1ee39b3a	247	fEventInfo->Delete("");
1ee39b3a	248	new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
1ee39b3a	249
770382d9	250	Bool_t *trackMap(NULL);
770382d9	251	AliStack * mStack(NULL);
1ee39b3a	252	Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
	253	if(HasMCdata()){
	254	mStack = fMCev->Stack();
	255	if(!mStack){
	256	AliError("Failed retrieving MC Stack");
	257	return;
	258	}
	259	trackMap = new Bool_t[nTracksMC];
	260	memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
	261	}
	262
f232df0d	263	Double32_t dedx[100]; Int_t nSlices(0);
770382d9	264	Int_t nTRDout(0), nTRDin(0), nTPC(0), nclsTrklt;
1ee39b3a	265	AliDebug(2, Form("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC));
770382d9	266	AliESDtrack *esdTrack = NULL;
	267	AliESDfriendTrack *esdFriendTrack = NULL;
	268	TObject *calObject = NULL;
	269	AliTRDtrackV1 *track = NULL;
	270	AliTRDseedV1 *tracklet = NULL;
	271	AliTRDcluster *cl = NULL;
1ee39b3a	272	for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
	273	esdTrack = fESDev->GetTrack(itrk);
	274	AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));
3f19c2cd	275
770382d9	276	if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
	277	if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
	278	if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;
1ee39b3a	279
	280	// look at external track param
	281	const AliExternalTrackParam *op = esdTrack->GetOuterParam();
	282	Double_t xyz[3];
	283	if(op){
	284	op->GetXYZ(xyz);
	285	op->Global2LocalPosition(xyz, op->GetAlpha());
	286	AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));
	287	}
	288
	289	// read MC info
	290	Int_t fPdg = -1;
	291	Int_t label = -1; UInt_t alab=UINT_MAX;
	292	if(HasMCdata()){
	293	label = esdTrack->GetLabel();
	294	alab = TMath::Abs(label);
	295	// register the track
	296	if(alab < UInt_t(nTracksMC)){
	297	trackMap[alab] = kTRUE;
	298	} else {
	299	AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
	300	continue;
	301	}
770382d9	302	AliMCParticle *mcParticle = NULL;
1ee39b3a	303	if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
	304	AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
	305	continue;
	306	}
	307	fPdg = mcParticle->Particle()->GetPdgCode();
	308	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
770382d9	309	Int_t iref = 0; AliTrackReference *ref = NULL;
1ee39b3a	310	while(iref<nRefs){
	311	ref = mcParticle->GetTrackReference(iref);
	312	if(ref->LocalX() > fgkTPC) break;
	313	iref++;
	314	}
	315
	316	new(fTrackInfo) AliTRDtrackInfo();
	317	fTrackInfo->SetPDG(fPdg);
	318	fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
	319	Int_t jref = iref;//, kref = 0;
	320	while(jref<nRefs){
	321	ref = mcParticle->GetTrackReference(jref);
	322	if(ref->LocalX() > fgkTOF) break;
	323	AliDebug(4, Form(" trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));
	324	fTrackInfo->AddTrackRef(ref);
	325	jref++;
	326	}
	327	AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));
	328	} else {
	329	new (fTrackInfo) AliTRDtrackInfo();
	330	fTrackInfo->SetPDG(fPdg);
	331	}
	332
	333	// copy some relevant info to TRD track info
	334	fTrackInfo->SetStatus(esdTrack->GetStatus());
	335	fTrackInfo->SetTrackId(esdTrack->GetID());
	336	Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
	337	fTrackInfo->SetESDpid(p);
	338	fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
f232df0d	339	if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();
f232df0d	340	memset(dedx, 0, 100*sizeof(Double32_t));
92c40e64	341	Int_t in(0);
	342	for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)
	343	for(Int_t is=0; is<nSlices; is++)
	344	dedx[in++]=esdTrack->GetTRDslice(il, is);
	345	for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);
	346	fTrackInfo->SetSlices(in, dedx);
1ee39b3a	347	fTrackInfo->SetLabel(label);
	348	fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
	349	// some other Informations which we may wish to store in order to find problematic cases
	350	fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
	351	fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
	352	nclsTrklt = 0;
	353
	354
	355	// read REC info
	356	esdFriendTrack = fESDfriend->GetTrack(itrk);
	357	if(esdFriendTrack){
	358	Int_t icalib = 0;
	359	while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
	360	if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
	361	if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
1ee39b3a	362	AliDebug(4, Form("TRD track OK"));
	363	// Set the clusters to unused
	364	for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
	365	if(!(tracklet = track->GetTracklet(ipl))) continue;
	366	tracklet->ResetClusterIter();
	367	while((cl = tracklet->NextCluster())) cl->Use(0);
	368	}
	369	fTrackInfo->SetTrack(track);
	370	break;
	371	}
	372	AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));
	373	} else AliDebug(3, "No ESD friends");
	374	if(op) fTrackInfo->SetOuterParam(op);
	375
	376	if(DebugLevel() >= 1){
	377	AliTRDtrackInfo info(*fTrackInfo);
	378	(*DebugStream()) << "trackInfo"
	379	<< "TrackInfo.=" << &info
	380	<< "\n";
	381	info.Delete("");
	382	}
705f8b0a	383
	384	ULong_t status(esdTrack->GetStatus());
	385	if((status&AliESDtrack::kTPCout)){
	386	if(!esdTrack->GetKinkIndex(0)){ // Barrel Track Selection
	387	Bool_t selected(kTRUE);
	388	if(UseLocalTrkSelection()){
	389	if(esdTrack->Pt() < fgkPt) selected = kFALSE;
	390	if(TMath::Abs(esdTrack->Eta()) > fgkEta) selected = kFALSE;
	391	if(esdTrack->GetTPCNcls() < fgkNclTPC) selected = kFALSE;
	392	Float_t par[2], cov[3];
	393	esdTrack->GetImpactParameters(par, cov);
	394	if(IsCollision()){ // cuts on DCA
	395	if(TMath::Abs(par[0]) > fgkTrkDCAxy) selected = kFALSE;
	396	if(TMath::Abs(par[1]) > fgkTrkDCAz) selected = kFALSE;
	397	}
	398	}
	399	if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;
	400	if(selected) fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
	401	} else fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
	402	} else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)) fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
1ee39b3a	403	fTrackInfo->Delete("");
1ee39b3a	404	}
b4414720	405
770382d9	406	AliDebug(2, Form("%3d Tracks: TPCout[%d] TRDin[%d] TRDout[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTPC, nTRDin, nTRDout));
1ee39b3a	407
770382d9	408	// AliESDv0 *v0 = NULL;
1ee39b3a	409	// for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
1ee39b3a	410	// if(!(v0 = fESD->GetV0(iv0))) continue;
705f8b0a	411	// fV0List->Add(new AliTRDv0Info(v0));
1ee39b3a	412	// }
	413
	414	// Insert also MC tracks which are passing TRD where the track is not reconstructed
	415	if(HasMCdata()){
	416	AliDebug(10, "Output of the MC track map:");
	417	for(Int_t itk = 0; itk < nTracksMC; itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));
	418
	419	for(Int_t itk = 0; itk < nTracksMC; itk++){
	420	if(trackMap[itk]) continue;
	421	AliMCParticle mcParticle = (AliMCParticle) fMCev->GetTrack(TMath::Abs(itk));
	422	Int_t fPdg = mcParticle->Particle()->GetPdgCode();
	423	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
770382d9	424	Int_t iref = 0; AliTrackReference *ref = NULL;
1ee39b3a	425	Int_t nRefsTRD = 0;
	426	new(fTrackInfo) AliTRDtrackInfo();
	427	fTrackInfo->SetPDG(fPdg);
	428	while(iref<nRefs){
770382d9	429	Bool_t kIN(kFALSE);
1ee39b3a	430	ref = mcParticle->GetTrackReference(iref);
770382d9	431	if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTOF){
1ee39b3a	432	fTrackInfo->AddTrackRef(ref);
770382d9	433	nRefsTRD++;kIN=kTRUE;
1ee39b3a	434	}
770382d9	435	AliDebug(4, Form(" trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
1ee39b3a	436	iref++;
	437	}
	438	if(!nRefsTRD){
	439	// In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
	440	// analysis job
	441	fTrackInfo->Delete("");
	442	continue;
	443	}
	444	fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
	445	fTrackInfo->SetLabel(itk);
	446	if(DebugLevel() >= 1){
	447	AliTRDtrackInfo info(*fTrackInfo);
	448	(*DebugStream()) << "trackInfo"
	449	<< "TrackInfo.=" << &info
	450	<< "\n";
	451	info.Delete("");
	452	}
770382d9	453	AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
705f8b0a	454	fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
1ee39b3a	455	fTrackInfo->Delete("");
	456	}
	457	delete[] trackMap;
	458	}
705f8b0a	459	PostData(kTracksBarrel, fTracksBarrel);
	460	PostData(kTracksSA, fTracksSA);
	461	PostData(kTracksKink, fTracksKink);
	462	PostData(kEventInfo, fEventInfo);
	463	PostData(kV0List, fV0List);
	464	}
	465
	466	//____________________________________________________________________
	467	void AliTRDinfoGen::SetTrigger(const Char_t *trigger)
	468	{
	469	if(!fEvTrigger) fEvTrigger = new TString(trigger);
	470	else (*fEvTrigger) = trigger;
	471	}
	472
	473	//____________________________________________________________________
	474	TTreeSRedirector* AliTRDinfoGen::DebugStream()
	475	{
	476	if(!fDebugStream){
	477	TDirectory *savedir = gDirectory;
	478	fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");
	479	savedir->cd();
	480	}
	481	return fDebugStream;
1ee39b3a	482	}
1ee39b3a	483
705f8b0a	484