[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDefficiencyMC.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: AliTRDefficiencyMC.cxx 27496 2008-07-22 08:35:45Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  Reconstruction QA                                                     //
//                                                                        //
//  Authors:                                                              //
//    Markus Fasel <M.Fasel@gsi.de>                                       //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <TObjArray.h>
#include <TClonesArray.h>
#include <TProfile.h>
#include <TMath.h>
#include "TTreeStream.h"

#include "AliMagF.h"
#include "AliPID.h"
#include "AliMathBase.h"
#include "AliTrackReference.h"
#include "AliAnalysisManager.h"

#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "Cal/AliTRDCalPID.h"
#include "info/AliTRDtrackInfo.h"
#include "AliTRDinfoGen.h"
#include "AliTRDefficiencyMC.h"

ClassImp(AliTRDefficiencyMC)

//_____________________________________________________________________________
AliTRDefficiencyMC::AliTRDefficiencyMC()
  :AliTRDrecoTask("EfficiencyMC", "Combined Tracking Efficiency")
{
  //
  // Default constructor
  //
}


//_____________________________________________________________________________
void AliTRDefficiencyMC::CreateOutputObjects(){
  //
  // Create output objects
  //

  OpenFile(0, "RECREATE");
  fContainer = Histos();

}

//_____________________________________________________________________________
void AliTRDefficiencyMC::Exec(Option_t *){
  //
  // Execute the task:
  //
  // Loop over TrackInfos
  // 1st: check if there is a trackTRD
  // 2nd: put conditions on the track:
  //      - check if we did not register it before
  //      - check if we have Track References for the track 
  // 3rd: Register track:
  //      - accepted if both conditions are fulfilled
  //      - contamination if at least one condition is not fulfilled
  // 4th: check Monte-Carlo Track wheter findable or not if there is no TRD track in this track info
  // 5th: register MC track as rejected if findable and not jet registered
  // Match the registers accepted and rejected and clean register rejected
  // Fill the histograms
  //
  const Int_t kArraySize = 10000;     // Store indices of track references in an array
  Int_t indexAccepted[kArraySize], indexRejected[kArraySize], indexContamination[kArraySize];
  memset(indexAccepted, 0, sizeof(Int_t) * kArraySize);
  memset(indexRejected, 0, sizeof(Int_t) * kArraySize);
  memset(indexContamination, 0, sizeof(Int_t) * kArraySize);
  Int_t naccepted = 0, nrejected = 0, ncontamination = 0;
  Bool_t isContamination = kFALSE;
  
  Int_t nTrackInfos = fTracks->GetEntriesFast();
  AliTRDtrackV1 *trackTRD = 0x0;
  AliTRDtrackInfo *trkInf = 0x0;
  for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
    trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
    if(!trkInf) continue;
    if(trkInf->GetTrack() || trkInf->GetNumberOfClustersRefit()){
      isContamination = IsRegistered(trkInf,indexAccepted,naccepted);
      if(!trkInf->GetNTrackRefs()){
        // We reject the track since the Monte Carlo Information is missing
        AliDebug(1, Form("Error: Track Reference missing for Track %d", trkInf->GetLabel()));
        isContamination = kTRUE;
        // Debugging
        if(trackTRD && DebugLevel() > 5) FillStreamTrackWOMC(trkInf);
      }	
      if(isContamination){
        // reject kink (we count these only once)
        if(trkInf->GetKinkIndex()) continue;
        // Register track as contamination
        indexContamination[ncontamination++]=itinf;
        continue;
      }
      // Accept track
      AliDebug(4, "Accept track");
      // Register track as accepted
      indexAccepted[naccepted++] = itinf;
    }else{
      if(IsFindable(trkInf)){
        // register track as rejected if not already registered there
        // Attention:
        // These track infos are not!!! registered as contamination
        if(!IsRegistered(trkInf, indexRejected, nrejected)) indexRejected[nrejected++] = itinf;
      }
    }
  }
  // we have to check if the rejected tracks are registered as found
  // a possible source for this:
  // first the track is not found by the barrel tracking but it is later found
  // by the stand alone tracking, then two track info objects with the same 
  // label would be created
  // Attention:
  // these tracks are not! registered as contamination
  Int_t tmprejected[kArraySize]; Int_t nrej = nrejected;
  memcpy(tmprejected, indexRejected, sizeof(Int_t) * nrejected);
  nrejected = 0;
  for(Int_t irej = 0; irej < nrej; irej++){
    trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej]));
    if(!IsRegistered(trkInf,indexAccepted,naccepted)) indexRejected[nrejected++] = tmprejected[irej];
  }
  // Fill Histograms
  FillHistograms(naccepted, &indexAccepted[0], kAccepted);
  FillHistograms(nrejected, &indexRejected[0], kRejected);
  FillHistograms(ncontamination, &indexContamination[0], kContamination);
  Int_t nall = naccepted + nrejected;
  //if(DebugLevel()>=1)
  AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, nall ? 1.E2*Float_t(naccepted)/Float_t(nall) : 0.));
  AliInfo(Form("%3d Tracks: ALL[%3d] Contamination[%3d | %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ncontamination, ncontamination, nall ? 1.E2*Float_t(ncontamination)/Float_t(nall + ncontamination) : 0.));

  PostData(0, fContainer);
}


//_____________________________________________________________________________
Bool_t AliTRDefficiencyMC::PostProcess()
{
  //
  // Post Process 
  //
  // Change the histogram style
  // For species histograms apply the colors connected with the given particle species
  //
  TH1 *histo = dynamic_cast<TH1 *>(fContainer->At(kEfficiencyHistogram));
  histo->SetMarkerStyle(22);
  histo->SetMarkerColor(kBlue);
  histo->GetXaxis()->SetTitle("p [GeV/c]");
  histo->GetXaxis()->SetMoreLogLabels();
  histo->GetYaxis()->SetTitle("Efficiency [%]");
  histo->GetYaxis()->SetRangeUser(0.99, 1.005);

  histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationHistogram));
  histo->SetMarkerStyle(22);
  histo->SetMarkerColor(kBlue);
  histo->GetXaxis()->SetTitle("p [GeV/c]");
  histo->GetXaxis()->SetMoreLogLabels();
  histo->GetYaxis()->SetTitle("Contamination [%]");
  
  // Species Efficiency Histograms
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    histo = dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram  + ispec));
    histo->SetMarkerStyle(22);
    histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
    histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
    histo->GetXaxis()->SetTitle("p [GeV/c]");
    histo->GetXaxis()->SetMoreLogLabels();
    histo->GetYaxis()->SetTitle("Efficiency [%]");
    histo->GetYaxis()->SetRangeUser(0.99, 1.005);
  }
  
  // Species Contamination Histograms
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram  + ispec));
    histo->SetMarkerStyle(22);
    histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
    histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
    histo->GetXaxis()->SetTitle("p [GeV/c]");
    histo->GetXaxis()->SetMoreLogLabels();
    histo->GetYaxis()->SetTitle("Contamination [%]");
  }
  
  fNRefFigures = 6;
  return kTRUE;
}

//_____________________________________________________________________________
Bool_t AliTRDefficiencyMC::GetRefFigure(Int_t ifig){
  //
  // Plot the histograms
  //
  if(ifig >= fNRefFigures) return kFALSE;
  if(ifig < 2){
    (dynamic_cast<TH1 *>(fContainer->At(ifig)))->Draw("e1");
    return kTRUE;
  }
  switch(ifig){
  case 2:
    (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram)))->Draw("e1");
    for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
      (dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram + ispec)))->Draw("e1same");
    break;
  case 3:
    (dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram)))->Draw("e1");
    for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
      (dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram + ispec)))->Draw("e1same");
    break;
  case 4:
    (dynamic_cast<TH1 *>(fContainer->At(kEfficiencyNoPID)))->Draw("e1");
    break;
  case 5:
    (dynamic_cast<TH1 *>(fContainer->At(kContaminationNoPID)))->Draw("e1");
    break;
  }
  return kTRUE;
}

//_____________________________________________________________________________
TObjArray *AliTRDefficiencyMC::Histos(){
  //
  // Create the histograms
  //
  const Int_t nbins = 11;

  if(fContainer) return fContainer;
  Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
  
  TString species[AliPID::kSPECIES] = {"Pions", "Muons", "Electrons", "Kaons", "Protons"};
  TString species_short[AliPID::kSPECIES] = {"Pi", "Mu", "El", "Ka", "Pr"};
  
  fContainer = new TObjArray();
  fContainer->AddAt(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyHistogram);
  fContainer->AddAt(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins), kContaminationHistogram);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
    fContainer->AddAt(new TProfile(Form("trEffComb%s", species_short[ispec].Data()), Form("Combined Tracking Efficiency %s", species[ispec].Data()), nbins, xbins), kEfficiencySpeciesHistogram + ispec);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
    fContainer->AddAt(new TProfile(Form("trContComb%s", species_short[ispec].Data()), Form("Combined Tracking Contamination %s", species[ispec].Data()), nbins, xbins), kContaminationSpeciesHistogram + ispec);
  fContainer->AddAt(new TProfile("trEffCombNoPID", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyNoPID);
  fContainer->AddAt(new TProfile("trContCombNoPID", "Combined Tracking Contamination", nbins, xbins), kContaminationNoPID);
  return fContainer;
}

//_____________________________________________________________________________
Bool_t AliTRDefficiencyMC::IsFindable(AliTRDtrackInfo * const trkInf){
  //
  // Apply Cuts on the Monte Carlo track references
  // return whether track is findable or not
  //
  const Float_t kAlpha = 0.349065850;
  
  AliDebug(10, "Analysing Track References\n");
  // Check if track is findable
  Double_t mom = 0.;
  Float_t xmin = 10000.0, xmax = 0.0; 
  Float_t ymin = 0.0, ymax = 0.0;
  Float_t zmin = 0.0, zmax = 0.0;
  Float_t lastx = 0.0, x = 0.0;
  Int_t nLayers = 0;
  Int_t sector[20];
  AliTrackReference *trackRef = 0x0;
  for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
    trackRef = trkInf->GetTrackRef(itr);
    AliDebug(10, Form("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z()));
    x = trackRef->LocalX(); 
        
    // Be Sure that we are inside TRD
    if(x < AliTRDinfoGen::GetTPCx() || x > AliTRDinfoGen::GetTOFx()) continue;	
    sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
    if(x < xmin){
      xmin = trackRef->LocalX();
      ymin = trackRef->LocalY();
      zmin = trackRef->Z();
      mom = trackRef->P();
    } else if(x > xmax){
      xmax = trackRef->LocalX();
      ymax = trackRef->LocalY();
      zmax = trackRef->Z();
    }
    if(itr > 0){
      Float_t dist = TMath::Abs(x - lastx);
      AliDebug(10, Form("x = %f, lastx = %f, dist = %f\n", x, lastx, dist));
      if(TMath::Abs(dist - 3.7) < 0.1) nLayers++; 	// ref(i+1) has to be larger than ref(i)
    }
    lastx = x;
  }

  // Apply cuts
  Bool_t findable = kTRUE;
  if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
    if(mom < 0.55) findable = kFALSE;									// momentum cut at 0.6
    Double_t yangle = (ymax -ymin)/(xmax - xmin);
    Double_t zangle = (zmax -zmin)/(xmax - xmin);
    AliDebug(10, Form("\ntrack: y-Angle = %f, z-Angle = %f\nnLayers = %d", yangle, zangle, nLayers));
    if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
    if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
    if(nLayers < 4) findable = kFALSE;
    if(!trkInf->IsPrimary()) findable = kFALSE;
    Bool_t samesec = kTRUE;
    for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
      if(sector[iref] != sector[0]) samesec = kFALSE;
    if(!samesec) findable = kFALSE;		// Discard all tracks which are passing more than one sector
    if(DebugLevel()){
      Double_t trackAngle = TMath::ATan(yangle);
      Bool_t primary = trkInf->IsPrimary();
      (*DebugStream()) << "NotFoundTrack"
        << "Momentum=" 	<< mom
        << "trackAngle="<< trackAngle
        << "NLayers="	<< nLayers
        << "Primary="	<< primary
        << "\n";
    }
  }
  else
    findable = kFALSE;
  return findable;
}

//_____________________________________________________________________________
void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, FillingMode_t mode){
  //
  // Fill Histograms in three different modes:
  // 1st tracks which are found and accepted
  // 2nd tracks which are not found and not already counted
  // 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
  //
  const Int_t pid[AliPID::kSPECIES] = {211,13,11,321,2212};
  Double_t trkmom = 0.;   // the track momentum
  Int_t trkpid = -1;      // particle species
  AliTRDtrackInfo *trkInf = 0x0;
  for(Int_t itk = 0; itk < nTracks; itk++){
    trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk]));
    AliDebug(3, Form("Accepted MC track: %d\n", trkInf->GetLabel()));
    if(trkInf->GetNTrackRefs()){
      // use Monte-Carlo Information for Momentum and PID
      trkmom = trkInf->GetTrackRef(0)->P();
      trkpid = trkInf->GetPDG();
    }else{
      // Use TPC Momentum
      trkmom = trkInf->GetTrack()->P();
    }
    switch(mode){
      case kAccepted:
        (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 1);
        (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
        break;
      case kRejected:
        (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 0);
        (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
        break;
      case kContamination:
        (dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 1);
        break;
    }
    // Fill species histogram
    Int_t partSpec = -1;
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      if(trkpid == pid[ispec]) partSpec = ispec;
    }
    if(partSpec >= 0){
      switch(mode){
        case kAccepted:
          (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 1);
          (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
          break;
        case kRejected:
          (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 0);          (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
          break;
        case kContamination:
          (dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 1);
          break;
      }
    } else {
      // The particle Type is not registered
      (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyNoPID)))->Fill(trkmom, 1);
      (dynamic_cast<TProfile *>(fContainer->At(kContaminationNoPID)))->Fill(trkmom, 1);
    }
  }
}

//_____________________________________________________________________________
void AliTRDefficiencyMC::FillStreamTrackWOMC(AliTRDtrackInfo * const trkInf){
  // fill debug stream
  // we want to know:
  //  1. The event number
  //  2. The track label
  //  3. The TRD track label
  //  4. The frequency of the TRD Label
  //  5. Momentum from TPC (NO STAND ALONE TRACK)
  //  6. TPC Phi angle
  //  7. the TRD track
  //  8. Monte Carlo PID
  //  9. We check the Labels of the TRD track according to them we search the maching Monte-Carlo track.
  //     From the matching Monte-Carlo track we store trackRefs, phi and momentum
  // 10. We may also want to keep the kink index
  Double_t mom = trkInf->GetESDinfo()->GetOuterParam()->P();
  Int_t event = (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry();
  Int_t label = trkInf->GetLabel();
  Int_t kinkIndex = trkInf->GetKinkIndex();
  Int_t pdg = trkInf->GetPDG();
  Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
  Int_t labelsTRD[180];	// Container for the cluster labels
  Int_t sortlabels[360];	// Cluster Labels sorted according their occurancy
  AliTRDseedV1 *tracklet = 0x0;
  AliTRDcluster *c = 0x0;
  Int_t nclusters = 0x0;
  AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
  for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
    tracklet = trackTRD->GetTracklet(il);
    if(!tracklet) continue;
    tracklet->ResetClusterIter();
    c = 0x0;
    while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
  }
  // Determine Label and Frequency
  AliMathBase::Freq(nclusters, const_cast<const Int_t *>(&labelsTRD[0]), &sortlabels[0], kTRUE);
  Int_t labelTRD = sortlabels[0];
  Int_t freqTRD = sortlabels[1];
  // find the track info object matching to the TRD track
  AliTRDtrackInfo *realtrack = 0;
  TObjArrayIter rtiter(fTracks);
  while((realtrack = (AliTRDtrackInfo *)rtiter())){
    if(realtrack->GetLabel() != labelTRD) continue;
    break;
  }
  TClonesArray trackRefs("AliTrackReference");
  Int_t realPdg = -1;
  Double_t realP = 0.;
  Double_t realPhi = 0.;
  if(realtrack){
    // pack the track references into the trackRefsContainer
    for(Int_t iref = 0; iref < realtrack->GetNTrackRefs(); iref++){
    new(trackRefs[iref])AliTrackReference(*(realtrack->GetTrackRef(iref)));
    }
    realPdg = realtrack->GetPDG();
    if(realtrack->GetNTrackRefs()){
      realP = realtrack->GetTrackRef(0)->P();
      realPhi = realtrack->GetTrackRef(0)->Phi();
    }
  }
  (*DebugStream()) << "TrackingEffMCfake"
    << "Event="	<< event
    << "Label=" << label
    << "labelTRD=" << labelTRD
    << "FreqTRDlabel=" << freqTRD
    << "TPCp="	<< mom
    << "phiTPC=" << phiTPC
    << "trackTRD=" << trackTRD
    << "PDG="	<< pdg
    << "TrackRefs=" << &trackRefs
    << "RealPDG=" << realPdg
    << "RealP="	<< realP
    << "RealPhi" << realPhi
    << "KinkIndex=" << kinkIndex
    << "\n";
}

//_____________________________________________________________________________
Bool_t AliTRDefficiencyMC::IsRegistered(AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
  //
  // Checks if track is registered in a given mode
  //
  Bool_t isRegistered = kFALSE;
  for(Int_t il = 0; il < nTracks; il++){
    if((dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))->GetLabel() == trkInf->GetLabel()){
      isRegistered = kTRUE;        
      break;
    }
  }
  return isRegistered;
}
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: AliTRDefficiencyMC.cxx 27496 2008-07-22 08:35:45Z cblume $ */
	17
	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Reconstruction QA //
	21	// //
	22	// Authors: //
	23	// Markus Fasel <M.Fasel@gsi.de> //
	24	// //
	25	////////////////////////////////////////////////////////////////////////////
	26
	27	#include <TObjArray.h>
	28	#include <TClonesArray.h>
	29	#include <TProfile.h>
	30	#include <TMath.h>
	31	#include "TTreeStream.h"
	32
	33	#include "AliMagF.h"
	34	#include "AliPID.h"
	35	#include "AliMathBase.h"
	36	#include "AliTrackReference.h"
	37	#include "AliAnalysisManager.h"
	38
	39	#include "AliTRDcluster.h"
	40	#include "AliTRDseedV1.h"
	41	#include "AliTRDtrackV1.h"
	42	#include "Cal/AliTRDCalPID.h"
	43	#include "info/AliTRDtrackInfo.h"
	44	#include "AliTRDinfoGen.h"
	45	#include "AliTRDefficiencyMC.h"
	46
	47	ClassImp(AliTRDefficiencyMC)
	48
	49	//_____________________________________________________________________________
	50	AliTRDefficiencyMC::AliTRDefficiencyMC()
	51	:AliTRDrecoTask("EfficiencyMC", "Combined Tracking Efficiency")
	52	{
	53	//
	54	// Default constructor
	55	//
	56	}
	57
	58
	59	//_____________________________________________________________________________
	60	void AliTRDefficiencyMC::CreateOutputObjects(){
	61	//
	62	// Create output objects
	63	//
	64
65	OpenFile(0, "RECREATE");
66	fContainer = Histos();
67
68	}
69
70	//_____________________________________________________________________________
71	void AliTRDefficiencyMC::Exec(Option_t *){
72	//
73	// Execute the task:
74	//
75	// Loop over TrackInfos
76	// 1st: check if there is a trackTRD
77	// 2nd: put conditions on the track:
78	// - check if we did not register it before
79	// - check if we have Track References for the track
80	// 3rd: Register track:
81	// - accepted if both conditions are fulfilled
82	// - contamination if at least one condition is not fulfilled
83	// 4th: check Monte-Carlo Track wheter findable or not if there is no TRD track in this track info
84	// 5th: register MC track as rejected if findable and not jet registered
85	// Match the registers accepted and rejected and clean register rejected
86	// Fill the histograms
87	//
88	const Int_t kArraySize = 10000; // Store indices of track references in an array
89	Int_t indexAccepted[kArraySize], indexRejected[kArraySize], indexContamination[kArraySize];
90	memset(indexAccepted, 0, sizeof(Int_t) * kArraySize);
91	memset(indexRejected, 0, sizeof(Int_t) * kArraySize);
92	memset(indexContamination, 0, sizeof(Int_t) * kArraySize);
93	Int_t naccepted = 0, nrejected = 0, ncontamination = 0;
94	Bool_t isContamination = kFALSE;
95
96	Int_t nTrackInfos = fTracks->GetEntriesFast();
97	AliTRDtrackV1 *trackTRD = 0x0;
98	AliTRDtrackInfo *trkInf = 0x0;
99	for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
100	trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
101	if(!trkInf) continue;
102	if(trkInf->GetTrack() \|\| trkInf->GetNumberOfClustersRefit()){
103	isContamination = IsRegistered(trkInf,indexAccepted,naccepted);
104	if(!trkInf->GetNTrackRefs()){
105	// We reject the track since the Monte Carlo Information is missing
106	AliDebug(1, Form("Error: Track Reference missing for Track %d", trkInf->GetLabel()));
107	isContamination = kTRUE;
108	// Debugging
109	if(trackTRD && DebugLevel() > 5) FillStreamTrackWOMC(trkInf);
110	}
111	if(isContamination){
112	// reject kink (we count these only once)
113	if(trkInf->GetKinkIndex()) continue;
114	// Register track as contamination
115	indexContamination[ncontamination++]=itinf;
116	continue;
117	}
118	// Accept track
119	AliDebug(4, "Accept track");
120	// Register track as accepted
121	indexAccepted[naccepted++] = itinf;
122	}else{
123	if(IsFindable(trkInf)){
124	// register track as rejected if not already registered there
125	// Attention:
126	// These track infos are not!!! registered as contamination
127	if(!IsRegistered(trkInf, indexRejected, nrejected)) indexRejected[nrejected++] = itinf;
128	}
129	}
130	}
131	// we have to check if the rejected tracks are registered as found
132	// a possible source for this:
133	// first the track is not found by the barrel tracking but it is later found
134	// by the stand alone tracking, then two track info objects with the same
135	// label would be created
136	// Attention:
137	// these tracks are not! registered as contamination
138	Int_t tmprejected[kArraySize]; Int_t nrej = nrejected;
139	memcpy(tmprejected, indexRejected, sizeof(Int_t) * nrejected);
140	nrejected = 0;
141	for(Int_t irej = 0; irej < nrej; irej++){
142	trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej]));
143	if(!IsRegistered(trkInf,indexAccepted,naccepted)) indexRejected[nrejected++] = tmprejected[irej];
144	}
145	// Fill Histograms
146	FillHistograms(naccepted, &indexAccepted[0], kAccepted);
147	FillHistograms(nrejected, &indexRejected[0], kRejected);
148	FillHistograms(ncontamination, &indexContamination[0], kContamination);
149	Int_t nall = naccepted + nrejected;
150	//if(DebugLevel()>=1)
151	AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d \| %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall, naccepted, nall ? 1.E2*Float_t(naccepted)/Float_t(nall) : 0.));
152	AliInfo(Form("%3d Tracks: ALL[%3d] Contamination[%3d \| %5.2f%%] \n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nall + ncontamination, ncontamination, nall ? 1.E2*Float_t(ncontamination)/Float_t(nall + ncontamination) : 0.));
153
154	PostData(0, fContainer);
155	}
156
157
158	//_____________________________________________________________________________
159	Bool_t AliTRDefficiencyMC::PostProcess()
160	{
161	//
162	// Post Process
163	//
164	// Change the histogram style
165	// For species histograms apply the colors connected with the given particle species
166	//
167	TH1 histo = dynamic_cast<TH1 >(fContainer->At(kEfficiencyHistogram));
168	histo->SetMarkerStyle(22);
169	histo->SetMarkerColor(kBlue);
170	histo->GetXaxis()->SetTitle("p [GeV/c]");
171	histo->GetXaxis()->SetMoreLogLabels();
172	histo->GetYaxis()->SetTitle("Efficiency [%]");
173	histo->GetYaxis()->SetRangeUser(0.99, 1.005);
174
175	histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationHistogram));
176	histo->SetMarkerStyle(22);
177	histo->SetMarkerColor(kBlue);
178	histo->GetXaxis()->SetTitle("p [GeV/c]");
179	histo->GetXaxis()->SetMoreLogLabels();
180	histo->GetYaxis()->SetTitle("Contamination [%]");
181
182	// Species Efficiency Histograms
183	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
184	histo = dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram + ispec));
185	histo->SetMarkerStyle(22);
186	histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
187	histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
188	histo->GetXaxis()->SetTitle("p [GeV/c]");
189	histo->GetXaxis()->SetMoreLogLabels();
190	histo->GetYaxis()->SetTitle("Efficiency [%]");
191	histo->GetYaxis()->SetRangeUser(0.99, 1.005);
192	}
193
194	// Species Contamination Histograms
195	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
196	histo = dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram + ispec));
197	histo->SetMarkerStyle(22);
198	histo->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
199	histo->SetMarkerColor(AliTRDCalPID::GetPartColor(ispec));
200	histo->GetXaxis()->SetTitle("p [GeV/c]");
201	histo->GetXaxis()->SetMoreLogLabels();
202	histo->GetYaxis()->SetTitle("Contamination [%]");
203	}
204
205	fNRefFigures = 6;
206	return kTRUE;
207	}
208
209	//_____________________________________________________________________________
210	Bool_t AliTRDefficiencyMC::GetRefFigure(Int_t ifig){
211	//
212	// Plot the histograms
213	//
214	if(ifig >= fNRefFigures) return kFALSE;
215	if(ifig < 2){
216	(dynamic_cast<TH1 *>(fContainer->At(ifig)))->Draw("e1");
217	return kTRUE;
218	}
219	switch(ifig){
220	case 2:
221	(dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram)))->Draw("e1");
222	for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
223	(dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram + ispec)))->Draw("e1same");
224	break;
225	case 3:
226	(dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram)))->Draw("e1");
227	for(Int_t ispec = 1; ispec < AliPID::kSPECIES; ispec++)
228	(dynamic_cast<TH1 *>(fContainer->At(kContaminationSpeciesHistogram + ispec)))->Draw("e1same");
229	break;
230	case 4:
231	(dynamic_cast<TH1 *>(fContainer->At(kEfficiencyNoPID)))->Draw("e1");
232	break;
233	case 5:
234	(dynamic_cast<TH1 *>(fContainer->At(kContaminationNoPID)))->Draw("e1");
235	break;
236	}
237	return kTRUE;
238	}
239
240	//_____________________________________________________________________________
241	TObjArray *AliTRDefficiencyMC::Histos(){
242	//
243	// Create the histograms
244	//
245	const Int_t nbins = 11;
246
247	if(fContainer) return fContainer;
248	Float_t xbins[nbins+1] = {.5, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
249
250	TString species[AliPID::kSPECIES] = {"Pions", "Muons", "Electrons", "Kaons", "Protons"};
251	TString species_short[AliPID::kSPECIES] = {"Pi", "Mu", "El", "Ka", "Pr"};
252
253	fContainer = new TObjArray();
254	fContainer->AddAt(new TProfile("trEffComb", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyHistogram);
255	fContainer->AddAt(new TProfile("trContComb", "Combined Tracking Contamination", nbins, xbins), kContaminationHistogram);
256	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
257	fContainer->AddAt(new TProfile(Form("trEffComb%s", species_short[ispec].Data()), Form("Combined Tracking Efficiency %s", species[ispec].Data()), nbins, xbins), kEfficiencySpeciesHistogram + ispec);
258	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
259	fContainer->AddAt(new TProfile(Form("trContComb%s", species_short[ispec].Data()), Form("Combined Tracking Contamination %s", species[ispec].Data()), nbins, xbins), kContaminationSpeciesHistogram + ispec);
260	fContainer->AddAt(new TProfile("trEffCombNoPID", "Combined Tracking Efficiency", nbins, xbins), kEfficiencyNoPID);
261	fContainer->AddAt(new TProfile("trContCombNoPID", "Combined Tracking Contamination", nbins, xbins), kContaminationNoPID);
262	return fContainer;
263	}
264
265	//_____________________________________________________________________________
266	Bool_t AliTRDefficiencyMC::IsFindable(AliTRDtrackInfo * const trkInf){
267	//
268	// Apply Cuts on the Monte Carlo track references
269	// return whether track is findable or not
270	//
271	const Float_t kAlpha = 0.349065850;
272
273	AliDebug(10, "Analysing Track References\n");
274	// Check if track is findable
275	Double_t mom = 0.;
276	Float_t xmin = 10000.0, xmax = 0.0;
277	Float_t ymin = 0.0, ymax = 0.0;
278	Float_t zmin = 0.0, zmax = 0.0;
279	Float_t lastx = 0.0, x = 0.0;
280	Int_t nLayers = 0;
281	Int_t sector[20];
282	AliTrackReference *trackRef = 0x0;
283	for(Int_t itr = 0; itr < trkInf->GetNTrackRefs(); itr++){
284	trackRef = trkInf->GetTrackRef(itr);
285	AliDebug(10, Form("%d. x[%f], y[%f], z[%f]\n", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z()));
286	x = trackRef->LocalX();
287
288	// Be Sure that we are inside TRD
289	if(x < AliTRDinfoGen::GetTPCx() \|\| x > AliTRDinfoGen::GetTOFx()) continue;
290	sector[itr] = Int_t(trackRef->Alpha()/kAlpha);
291	if(x < xmin){
292	xmin = trackRef->LocalX();
293	ymin = trackRef->LocalY();
294	zmin = trackRef->Z();
295	mom = trackRef->P();
296	} else if(x > xmax){
297	xmax = trackRef->LocalX();
298	ymax = trackRef->LocalY();
299	zmax = trackRef->Z();
300	}
301	if(itr > 0){
302	Float_t dist = TMath::Abs(x - lastx);
303	AliDebug(10, Form("x = %f, lastx = %f, dist = %f\n", x, lastx, dist));
304	if(TMath::Abs(dist - 3.7) < 0.1) nLayers++; // ref(i+1) has to be larger than ref(i)
305	}
306	lastx = x;
307	}
308
309	// Apply cuts
310	Bool_t findable = kTRUE;
311	if(trkInf->GetNTrackRefs() > 2 && xmax > xmin){
312	if(mom < 0.55) findable = kFALSE; // momentum cut at 0.6
313	Double_t yangle = (ymax -ymin)/(xmax - xmin);
314	Double_t zangle = (zmax -zmin)/(xmax - xmin);
315	AliDebug(10, Form("\ntrack: y-Angle = %f, z-Angle = %f\nnLayers = %d", yangle, zangle, nLayers));
316	if(TMath::ATan(TMath::Abs(yangle)) > 45.) findable = kFALSE;
317	if(TMath::ATan(TMath::Abs(zangle)) > 45.) findable = kFALSE;
318	if(nLayers < 4) findable = kFALSE;
319	if(!trkInf->IsPrimary()) findable = kFALSE;
320	Bool_t samesec = kTRUE;
321	for(Int_t iref = 1; iref < trkInf->GetNTrackRefs(); iref++)
322	if(sector[iref] != sector[0]) samesec = kFALSE;
323	if(!samesec) findable = kFALSE; // Discard all tracks which are passing more than one sector
324	if(DebugLevel()){
325	Double_t trackAngle = TMath::ATan(yangle);
326	Bool_t primary = trkInf->IsPrimary();
327	(*DebugStream()) << "NotFoundTrack"
328	<< "Momentum=" << mom
329	<< "trackAngle="<< trackAngle
330	<< "NLayers=" << nLayers
331	<< "Primary=" << primary
332	<< "\n";
333	}
334	}
335	else
336	findable = kFALSE;
337	return findable;
338	}
339
340	//_____________________________________________________________________________
341	void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, FillingMode_t mode){
342	//
343	// Fill Histograms in three different modes:
344	// 1st tracks which are found and accepted
345	// 2nd tracks which are not found and not already counted
346	// 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
347	//
348	const Int_t pid[AliPID::kSPECIES] = {211,13,11,321,2212};
349	Double_t trkmom = 0.; // the track momentum
350	Int_t trkpid = -1; // particle species
351	AliTRDtrackInfo *trkInf = 0x0;
352	for(Int_t itk = 0; itk < nTracks; itk++){
353	trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk]));
354	AliDebug(3, Form("Accepted MC track: %d\n", trkInf->GetLabel()));
355	if(trkInf->GetNTrackRefs()){
356	// use Monte-Carlo Information for Momentum and PID
357	trkmom = trkInf->GetTrackRef(0)->P();
358	trkpid = trkInf->GetPDG();
359	}else{
360	// Use TPC Momentum
361	trkmom = trkInf->GetTrack()->P();
362	}
363	switch(mode){
364	case kAccepted:
365	(dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 1);
366	(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
367	break;
368	case kRejected:
369	(dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))->Fill(trkmom, 0);
370	(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 0);
371	break;
372	case kContamination:
373	(dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))->Fill(trkmom, 1);
374	break;
375	}
376	// Fill species histogram
377	Int_t partSpec = -1;
378	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
379	if(trkpid == pid[ispec]) partSpec = ispec;
380	}
381	if(partSpec >= 0){
382	switch(mode){
383	case kAccepted:
384	(dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 1);
385	(dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
386	break;
387	case kRejected:
388	(dynamic_cast<TProfile >(fContainer->At(kEfficiencySpeciesHistogram + partSpec)))->Fill(trkmom, 0); (dynamic_cast<TProfile >(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 0);
389	break;
390	case kContamination:
391	(dynamic_cast<TProfile *>(fContainer->At(kContaminationSpeciesHistogram + partSpec)))->Fill(trkmom, 1);
392	break;
393	}
394	} else {
395	// The particle Type is not registered
396	(dynamic_cast<TProfile *>(fContainer->At(kEfficiencyNoPID)))->Fill(trkmom, 1);
397	(dynamic_cast<TProfile *>(fContainer->At(kContaminationNoPID)))->Fill(trkmom, 1);
398	}
399	}
400	}
401
402	//_____________________________________________________________________________
403	void AliTRDefficiencyMC::FillStreamTrackWOMC(AliTRDtrackInfo * const trkInf){
404	// fill debug stream
405	// we want to know:
406	// 1. The event number
407	// 2. The track label
408	// 3. The TRD track label
409	// 4. The frequency of the TRD Label
410	// 5. Momentum from TPC (NO STAND ALONE TRACK)
411	// 6. TPC Phi angle
412	// 7. the TRD track
413	// 8. Monte Carlo PID
414	// 9. We check the Labels of the TRD track according to them we search the maching Monte-Carlo track.
415	// From the matching Monte-Carlo track we store trackRefs, phi and momentum
416	// 10. We may also want to keep the kink index
417	Double_t mom = trkInf->GetESDinfo()->GetOuterParam()->P();
418	Int_t event = (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry();
419	Int_t label = trkInf->GetLabel();
420	Int_t kinkIndex = trkInf->GetKinkIndex();
421	Int_t pdg = trkInf->GetPDG();
422	Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
423	Int_t labelsTRD[180]; // Container for the cluster labels
424	Int_t sortlabels[360]; // Cluster Labels sorted according their occurancy
425	AliTRDseedV1 *tracklet = 0x0;
426	AliTRDcluster *c = 0x0;
427	Int_t nclusters = 0x0;
428	AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
429	for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
430	tracklet = trackTRD->GetTracklet(il);
431	if(!tracklet) continue;
432	tracklet->ResetClusterIter();
433	c = 0x0;
434	while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
435	}
436	// Determine Label and Frequency
437	AliMathBase::Freq(nclusters, const_cast<const Int_t *>(&labelsTRD[0]), &sortlabels[0], kTRUE);
438	Int_t labelTRD = sortlabels[0];
439	Int_t freqTRD = sortlabels[1];
440	// find the track info object matching to the TRD track
441	AliTRDtrackInfo *realtrack = 0;
442	TObjArrayIter rtiter(fTracks);
443	while((realtrack = (AliTRDtrackInfo *)rtiter())){
444	if(realtrack->GetLabel() != labelTRD) continue;
445	break;
446	}
447	TClonesArray trackRefs("AliTrackReference");
448	Int_t realPdg = -1;
449	Double_t realP = 0.;
450	Double_t realPhi = 0.;
451	if(realtrack){
452	// pack the track references into the trackRefsContainer
453	for(Int_t iref = 0; iref < realtrack->GetNTrackRefs(); iref++){
454	new(trackRefs[iref])AliTrackReference(*(realtrack->GetTrackRef(iref)));
455	}
456	realPdg = realtrack->GetPDG();
457	if(realtrack->GetNTrackRefs()){
458	realP = realtrack->GetTrackRef(0)->P();
459	realPhi = realtrack->GetTrackRef(0)->Phi();
460	}
461	}
462	(*DebugStream()) << "TrackingEffMCfake"
463	<< "Event=" << event
464	<< "Label=" << label
465	<< "labelTRD=" << labelTRD
466	<< "FreqTRDlabel=" << freqTRD
467	<< "TPCp=" << mom
468	<< "phiTPC=" << phiTPC
469	<< "trackTRD=" << trackTRD
470	<< "PDG=" << pdg
471	<< "TrackRefs=" << &trackRefs
472	<< "RealPDG=" << realPdg
473	<< "RealP=" << realP
474	<< "RealPhi" << realPhi
475	<< "KinkIndex=" << kinkIndex
476	<< "\n";
477	}
478
479	//_____________________________________________________________________________
480	Bool_t AliTRDefficiencyMC::IsRegistered(AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
481	//
482	// Checks if track is registered in a given mode
483	//
484	Bool_t isRegistered = kFALSE;
485	for(Int_t il = 0; il < nTracks; il++){
486	if((dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))->GetLabel() == trkInf->GetLabel()){
487	isRegistered = kTRUE;
488	break;
489	}
490	}
491	return isRegistered;
492	}
493