[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(fDebugStream) delete fDebugStream;
  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"
    "            Barrel[%d] SA[%d] Kink[%d]"
    ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry()
    , nTPC, nTRDin, nTRDout
    ,fTracksBarrel->GetEntriesFast()
    ,fTracksSA->GetEntriesFast()
    ,fTracksKink->GetEntriesFast()
  ));


//   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
0205fb04	149	if(fDebugStream) delete fDebugStream;
705f8b0a	150	if(fEvTrigger) delete fEvTrigger;
07ece31c	151	if(fTrackCut) delete fTrackCut;
07ece31c	152	if(fEventCut) delete fEventCut;
770382d9	153	if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;
770382d9	154	if(fEventInfo) delete fEventInfo; fEventInfo = NULL;
770382d9	155	if(fV0Info) delete fV0Info; fV0Info = NULL;
705f8b0a	156	if(fTracksBarrel){
	157	fTracksBarrel->Delete(); delete fTracksBarrel;
	158	fTracksBarrel = NULL;
	159	}
	160	if(fTracksSA){
	161	fTracksSA->Delete(); delete fTracksSA;
	162	fTracksSA = NULL;
1ee39b3a	163	}
705f8b0a	164	if(fTracksKink){
	165	fTracksKink->Delete(); delete fTracksKink;
	166	fTracksKink = NULL;
	167	}
	168	if(fV0List){
	169	fV0List->Delete(); delete fV0List;
	170	fV0List = NULL;
1ee39b3a	171	}
	172	}
	173
	174	//____________________________________________________________________
f8f46e4d	175	void AliTRDinfoGen::UserCreateOutputObjects()
1ee39b3a	176	{
	177	//
	178	// Create Output Containers (TObjectArray containing 1D histograms)
	179	//
b4414720	180
1ee39b3a	181	fTrackInfo = new AliTRDtrackInfo();
	182	fEventInfo = new AliTRDeventInfo();
	183	fV0Info = new AliTRDv0Info();
705f8b0a	184	fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);
	185	fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);
	186	fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);
	187	fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);
1ee39b3a	188	}
	189
	190	//____________________________________________________________________
f8f46e4d	191	void AliTRDinfoGen::UserExec(Option_t *){
1ee39b3a	192	//
	193	// Run the Analysis
	194	//
f2e89a4c	195
f8f46e4d	196	fESDev = dynamic_cast<AliESDEvent*>(InputEvent());
	197	fMCev = MCEvent();
	198
1ee39b3a	199	if(!fESDev){
	200	AliError("Failed retrieving ESD event");
	201	return;
	202	}
3f19c2cd	203
3f19c2cd	204	// event selection : trigger cut
705f8b0a	205	if(UseLocalEvSelection() && fEvTrigger){
3f19c2cd	206	Bool_t kTRIGGERED(kFALSE);
705f8b0a	207	const TObjArray *trig = fEvTrigger->Tokenize(" ");
3f19c2cd	208	for(Int_t itrig=trig->GetEntriesFast(); itrig--;){
	209	const Char_t trigClass(((TObjString)(*trig)[itrig])->GetName());
	210	if(fESDev->IsTriggerClassFired(trigClass)) {
	211	AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));
	212	kTRIGGERED = kTRUE;
	213	break;
	214	}
	215	}
	216	if(!kTRIGGERED) return;
	217
	218	// select only physical events
	219	if(fESDev->GetEventType() != 7) return;
	220	}
	221
	222	// if the required trigger is a collision trigger then apply event vertex cut
	223	if(UseLocalEvSelection() && IsCollision()){
	224	const AliESDVertex *vertex = fESDev->GetPrimaryVertex();
	225	if(TMath::Abs(vertex->GetZv())<1.e-10 \|\|
	226	TMath::Abs(vertex->GetZv())>fgkEvVertexZ \|\|
	227	vertex->GetNContributors()<fgkEvVertexN) {
	228	//PostData(0, fOutStorage);
	229	return;
	230	}
	231	}
	232
07ece31c	233	if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;
07ece31c	234
1ee39b3a	235	if(!fESDfriend){
	236	AliError("Failed retrieving ESD friend event");
	237	return;
	238	}
	239	if(HasMCdata() && !fMCev){
	240	AliError("Failed retrieving MC event");
	241	return;
	242	}
f8f46e4d	243
705f8b0a	244	fTracksBarrel->Delete();
	245	fTracksSA->Delete();
	246	fTracksKink->Delete();
	247	fV0List->Delete();
1ee39b3a	248	fEventInfo->Delete("");
1ee39b3a	249	new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
1ee39b3a	250
770382d9	251	Bool_t *trackMap(NULL);
770382d9	252	AliStack * mStack(NULL);
1ee39b3a	253	Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
	254	if(HasMCdata()){
	255	mStack = fMCev->Stack();
	256	if(!mStack){
	257	AliError("Failed retrieving MC Stack");
	258	return;
	259	}
	260	trackMap = new Bool_t[nTracksMC];
	261	memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
	262	}
	263
f232df0d	264	Double32_t dedx[100]; Int_t nSlices(0);
770382d9	265	Int_t nTRDout(0), nTRDin(0), nTPC(0), nclsTrklt;
1ee39b3a	266	AliDebug(2, Form("Entry[%3d] Tracks: ESD[%d] MC[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC));
770382d9	267	AliESDtrack *esdTrack = NULL;
	268	AliESDfriendTrack *esdFriendTrack = NULL;
	269	TObject *calObject = NULL;
	270	AliTRDtrackV1 *track = NULL;
	271	AliTRDseedV1 *tracklet = NULL;
	272	AliTRDcluster *cl = NULL;
1ee39b3a	273	for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
	274	esdTrack = fESDev->GetTrack(itrk);
	275	AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));
3f19c2cd	276
770382d9	277	if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
	278	if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
	279	if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;
1ee39b3a	280
	281	// look at external track param
	282	const AliExternalTrackParam *op = esdTrack->GetOuterParam();
	283	Double_t xyz[3];
	284	if(op){
	285	op->GetXYZ(xyz);
	286	op->Global2LocalPosition(xyz, op->GetAlpha());
	287	AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));
	288	}
	289
	290	// read MC info
	291	Int_t fPdg = -1;
	292	Int_t label = -1; UInt_t alab=UINT_MAX;
	293	if(HasMCdata()){
	294	label = esdTrack->GetLabel();
	295	alab = TMath::Abs(label);
	296	// register the track
	297	if(alab < UInt_t(nTracksMC)){
	298	trackMap[alab] = kTRUE;
	299	} else {
	300	AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
	301	continue;
	302	}
770382d9	303	AliMCParticle *mcParticle = NULL;
1ee39b3a	304	if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
	305	AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
	306	continue;
	307	}
	308	fPdg = mcParticle->Particle()->GetPdgCode();
	309	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
770382d9	310	Int_t iref = 0; AliTrackReference *ref = NULL;
1ee39b3a	311	while(iref<nRefs){
	312	ref = mcParticle->GetTrackReference(iref);
	313	if(ref->LocalX() > fgkTPC) break;
	314	iref++;
	315	}
	316
	317	new(fTrackInfo) AliTRDtrackInfo();
	318	fTrackInfo->SetPDG(fPdg);
	319	fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
	320	Int_t jref = iref;//, kref = 0;
	321	while(jref<nRefs){
	322	ref = mcParticle->GetTrackReference(jref);
	323	if(ref->LocalX() > fgkTOF) break;
	324	AliDebug(4, Form(" trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));
	325	fTrackInfo->AddTrackRef(ref);
	326	jref++;
	327	}
	328	AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));
	329	} else {
	330	new (fTrackInfo) AliTRDtrackInfo();
	331	fTrackInfo->SetPDG(fPdg);
	332	}
	333
	334	// copy some relevant info to TRD track info
	335	fTrackInfo->SetStatus(esdTrack->GetStatus());
	336	fTrackInfo->SetTrackId(esdTrack->GetID());
	337	Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
	338	fTrackInfo->SetESDpid(p);
	339	fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
f232df0d	340	if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();
f232df0d	341	memset(dedx, 0, 100*sizeof(Double32_t));
92c40e64	342	Int_t in(0);
	343	for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)
	344	for(Int_t is=0; is<nSlices; is++)
	345	dedx[in++]=esdTrack->GetTRDslice(il, is);
	346	for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);
	347	fTrackInfo->SetSlices(in, dedx);
1ee39b3a	348	fTrackInfo->SetLabel(label);
	349	fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
	350	// some other Informations which we may wish to store in order to find problematic cases
	351	fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
	352	fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
	353	nclsTrklt = 0;
	354
	355
	356	// read REC info
	357	esdFriendTrack = fESDfriend->GetTrack(itrk);
	358	if(esdFriendTrack){
	359	Int_t icalib = 0;
	360	while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
	361	if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
	362	if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
1ee39b3a	363	AliDebug(4, Form("TRD track OK"));
	364	// Set the clusters to unused
	365	for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
	366	if(!(tracklet = track->GetTracklet(ipl))) continue;
	367	tracklet->ResetClusterIter();
	368	while((cl = tracklet->NextCluster())) cl->Use(0);
	369	}
	370	fTrackInfo->SetTrack(track);
	371	break;
	372	}
	373	AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));
	374	} else AliDebug(3, "No ESD friends");
	375	if(op) fTrackInfo->SetOuterParam(op);
	376
	377	if(DebugLevel() >= 1){
	378	AliTRDtrackInfo info(*fTrackInfo);
	379	(*DebugStream()) << "trackInfo"
	380	<< "TrackInfo.=" << &info
	381	<< "\n";
	382	info.Delete("");
	383	}
705f8b0a	384
	385	ULong_t status(esdTrack->GetStatus());
	386	if((status&AliESDtrack::kTPCout)){
	387	if(!esdTrack->GetKinkIndex(0)){ // Barrel Track Selection
	388	Bool_t selected(kTRUE);
	389	if(UseLocalTrkSelection()){
	390	if(esdTrack->Pt() < fgkPt) selected = kFALSE;
	391	if(TMath::Abs(esdTrack->Eta()) > fgkEta) selected = kFALSE;
	392	if(esdTrack->GetTPCNcls() < fgkNclTPC) selected = kFALSE;
	393	Float_t par[2], cov[3];
	394	esdTrack->GetImpactParameters(par, cov);
	395	if(IsCollision()){ // cuts on DCA
	396	if(TMath::Abs(par[0]) > fgkTrkDCAxy) selected = kFALSE;
	397	if(TMath::Abs(par[1]) > fgkTrkDCAz) selected = kFALSE;
	398	}
	399	}
	400	if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;
	401	if(selected) fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
	402	} else fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
	403	} else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)) fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
1ee39b3a	404	fTrackInfo->Delete("");
1ee39b3a	405	}
afca20ef	406	AliDebug(2, Form(
	407	"%3d Tracks: TPCout[%d] TRDin[%d] TRDout[%d]\n"
	408	" Barrel[%d] SA[%d] Kink[%d]"
	409	,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry()
	410	, nTPC, nTRDin, nTRDout
	411	,fTracksBarrel->GetEntriesFast()
	412	,fTracksSA->GetEntriesFast()
	413	,fTracksKink->GetEntriesFast()
	414	));
b4414720	415
1ee39b3a	416
770382d9	417	// AliESDv0 *v0 = NULL;
1ee39b3a	418	// for(Int_t iv0=0; iv0<fESD->GetNumberOfV0s(); iv0++){
1ee39b3a	419	// if(!(v0 = fESD->GetV0(iv0))) continue;
705f8b0a	420	// fV0List->Add(new AliTRDv0Info(v0));
1ee39b3a	421	// }
	422
	423	// Insert also MC tracks which are passing TRD where the track is not reconstructed
	424	if(HasMCdata()){
	425	AliDebug(10, "Output of the MC track map:");
	426	for(Int_t itk = 0; itk < nTracksMC; itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));
	427
	428	for(Int_t itk = 0; itk < nTracksMC; itk++){
	429	if(trackMap[itk]) continue;
	430	AliMCParticle mcParticle = (AliMCParticle) fMCev->GetTrack(TMath::Abs(itk));
	431	Int_t fPdg = mcParticle->Particle()->GetPdgCode();
	432	Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
770382d9	433	Int_t iref = 0; AliTrackReference *ref = NULL;
1ee39b3a	434	Int_t nRefsTRD = 0;
	435	new(fTrackInfo) AliTRDtrackInfo();
	436	fTrackInfo->SetPDG(fPdg);
	437	while(iref<nRefs){
770382d9	438	Bool_t kIN(kFALSE);
1ee39b3a	439	ref = mcParticle->GetTrackReference(iref);
770382d9	440	if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTOF){
1ee39b3a	441	fTrackInfo->AddTrackRef(ref);
770382d9	442	nRefsTRD++;kIN=kTRUE;
1ee39b3a	443	}
770382d9	444	AliDebug(4, Form(" trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
1ee39b3a	445	iref++;
	446	}
	447	if(!nRefsTRD){
	448	// In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
	449	// analysis job
	450	fTrackInfo->Delete("");
	451	continue;
	452	}
	453	fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
	454	fTrackInfo->SetLabel(itk);
	455	if(DebugLevel() >= 1){
	456	AliTRDtrackInfo info(*fTrackInfo);
	457	(*DebugStream()) << "trackInfo"
	458	<< "TrackInfo.=" << &info
	459	<< "\n";
	460	info.Delete("");
	461	}
770382d9	462	AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
705f8b0a	463	fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
1ee39b3a	464	fTrackInfo->Delete("");
	465	}
	466	delete[] trackMap;
	467	}
705f8b0a	468	PostData(kTracksBarrel, fTracksBarrel);
	469	PostData(kTracksSA, fTracksSA);
	470	PostData(kTracksKink, fTracksKink);
	471	PostData(kEventInfo, fEventInfo);
	472	PostData(kV0List, fV0List);
	473	}
	474
	475	//____________________________________________________________________
	476	void AliTRDinfoGen::SetTrigger(const Char_t *trigger)
	477	{
	478	if(!fEvTrigger) fEvTrigger = new TString(trigger);
	479	else (*fEvTrigger) = trigger;
	480	}
	481
	482	//____________________________________________________________________
	483	TTreeSRedirector* AliTRDinfoGen::DebugStream()
	484	{
	485	if(!fDebugStream){
	486	TDirectory *savedir = gDirectory;
	487	fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");
	488	savedir->cd();
	489	}
	490	return fDebugStream;
1ee39b3a	491	}
1ee39b3a	492
705f8b0a	493