[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskMuonTrackingEff.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$ */

//Class to calculate the intrinsic efficiency of the detection elements of the
//MUON tracking chambers in function of the position in the detection element.
//WOrk on ESD only
//Author:  Nicolas LE BRIS - SUBATECH Nantes
// Modified by Matthieu LENHARDT - SUBATECH Nantes
// Modified by Antoine LARDEUX - SUBATECH Nantes

// ROOT includes
#include <TList.h>
#include <TH3F.h>
#include <TH2F.h>
#include <TObjArray.h>
#include <TGeoGlobalMagField.h>

// STEER includes
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliGeomManager.h"
#include "AliCDBManager.h"
#include "AliESDVZERO.h"

// ANALYSIS includes
#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskMuonTrackingEff.h"
#include "AliCentrality.h"

//MUON includes
#include "AliMUONCDB.h"
#include "AliMUONESDInterface.h"
#include "AliMUONGeometryTransformer.h"
#include "AliMUONTrack.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "AliMUONVCluster.h"
#include "AliMUONConstants.h"

//include MUON/mapping:
#include "AliMpDEManager.h"
#include "AliMpSegmentation.h"
#include "AliMpVSegmentation.h"
#include "AliMpPad.h"

ClassImp(AliAnalysisTaskMuonTrackingEff)

const Int_t AliAnalysisTaskMuonTrackingEff::fgkNbrOfDetectionElt[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26};
const Int_t AliAnalysisTaskMuonTrackingEff::fgkOffset = 100;

//________________________________________________________________________
AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff() :
  AliAnalysisTaskSE(),
  fOCDBLoaded(kFALSE),
  fOCDBpath(""),
  fMatchTrig(kFALSE),
  fApplyAccCut(kFALSE),
  fCurrentCentrality(0.),
  fTransformer(0x0),
  fDetEltTDHistList(0x0),
  fDetEltTTHistList(0x0),
  fChamberTDHistList(0x0),
  fChamberTTHistList(0x0)
{
  /// Default constructor
}

//________________________________________________________________________
AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff(TString name) :
  AliAnalysisTaskSE(name),
  fOCDBLoaded(kFALSE),
  fOCDBpath("raw://"),
  fMatchTrig(kFALSE),
  fApplyAccCut(kFALSE),
  fCurrentCentrality(100.),
  fTransformer(0x0),
  fDetEltTDHistList(0x0),
  fDetEltTTHistList(0x0),
  fChamberTDHistList(0x0),
  fChamberTTHistList(0x0)
{
  /// Constructor
  
  // Output slots 0 to 5 writes into a TClonesArray:
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  DefineOutput(3, TList::Class());
  DefineOutput(4, TList::Class());
}

//________________________________________________________________________
AliAnalysisTaskMuonTrackingEff::~AliAnalysisTaskMuonTrackingEff()
{
  /// Destructor
  if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
    delete fDetEltTDHistList;
    delete fDetEltTTHistList;
    delete fChamberTDHistList;
    delete fChamberTTHistList;
  }
  delete fTransformer;
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::NotifyRun()
{
  /// Load the OCDB and the Geometry
  
  // Load it only once
  if (fOCDBLoaded) return;
  
  // OCDB
  AliCDBManager* man = AliCDBManager::Instance();
  if (man->IsDefaultStorageSet()) printf("EfficiencyTask: CDB default storage already set!");
  else man->SetDefaultStorage(fOCDBpath.Data());
  if (man->GetRun() > -1) printf("EfficiencyTask: run number already set!");
  else man->SetRun(fCurrentRunNumber);
  
  // Geometry
  if (!AliGeomManager::GetGeometry()) {
    AliGeomManager::LoadGeometry();
    if (!AliGeomManager::GetGeometry()) return;  
    if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
  }
  fTransformer = new AliMUONGeometryTransformer();
  fTransformer->LoadGeometryData();
  
  // Magnetic field for track extrapolation
  if (!TGeoGlobalMagField::Instance()->GetField()) {
    if (!AliMUONCDB::LoadField()) return;
  }
  
  // Mapping
  if (!AliMpSegmentation::Instance(kFALSE)) {
    if (!AliMUONCDB::LoadMapping(kTRUE)) return;
  }
  
  // RecoParam for refitting
  AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
  if (!recoParam) return;
  AliMUONESDInterface::ResetTracker(recoParam);
  
  fOCDBLoaded = kTRUE;
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::UserCreateOutputObjects()
{
  /// Define efficiency histograms
  
  fDetEltTDHistList  = new TList();
  fDetEltTDHistList->SetOwner();
  fDetEltTTHistList  = new TList();
  fDetEltTTHistList->SetOwner();
  fChamberTDHistList = new TList();
  fChamberTDHistList->SetOwner();
  fChamberTTHistList = new TList();
  fChamberTTHistList->SetOwner();
  
  TH2F *h2;
  TH3F *h3;
  TString histNameTD;
  TString histNameTT;
  TString histTitle;
  Int_t nCentBins = 22;
  Double_t centRange[2] = {-5., 105.};
  Int_t iDEGlobal = 0;
  
  for (Int_t iCh = 0; iCh < 10; iCh++)
  {
    // histograms per chamber
    histTitle.Form("ChamberNbr %d", iCh+1);
    histNameTD.Form("TD_ChamberNbr%d", iCh+1);
    histNameTT.Form("TT_ChamberNbr%d",iCh+1);
    h2 = new TH2F(histNameTD, histTitle, fgkNbrOfDetectionElt[iCh], 0., fgkNbrOfDetectionElt[iCh], nCentBins, centRange[0], centRange[1]);
    fChamberTDHistList->AddAt(h2, iCh);
    h2 = new TH2F(histNameTT, histTitle, fgkNbrOfDetectionElt[iCh], 0., fgkNbrOfDetectionElt[iCh], nCentBins, centRange[0], centRange[1]);
    fChamberTTHistList->AddAt(h2, iCh);
    
    // histograms per DE
    for (Int_t iDE = 0; iDE < fgkNbrOfDetectionElt[iCh]; iDE++)
    {
      Int_t deId = FromLocalId2DetElt(iCh, iDE);
      histTitle.Form("detEltNbr %d",deId);
      histNameTD.Form("TD_detEltNbr%d",deId);
      histNameTT.Form("TT_detEltNbr%d",deId);
      if(iCh < 4)
      {// chambers 1 -> 4
	h3 = new TH3F(histNameTD, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
	fDetEltTDHistList->AddAt(h3, iDEGlobal);
	h3 = new TH3F(histNameTT, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
	fDetEltTTHistList->AddAt(h3, iDEGlobal);
      }
      else 
      {// chambers 5 -> 10
	h3 = new TH3F(histNameTD, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
	fDetEltTDHistList->AddAt(h3, iDEGlobal);	  
	h3 = new TH3F(histNameTT, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
	fDetEltTTHistList->AddAt(h3, iDEGlobal);	  
      }
      iDEGlobal++;
    }
  }
  
  // global histograms per chamber
  h2 = new TH2F("TD_Chambers 11", "Chambers 11", 10, 0.5, 10.5, nCentBins, centRange[0], centRange[1]);
  fChamberTDHistList->AddAt(h2, 10);
  h2 = new TH2F("TT_Chambers 11", "Chambers 11", 10, 0.5, 10.5, nCentBins, centRange[0], centRange[1]);
  fChamberTTHistList->AddAt(h2, 10);
  
  // post the output data at least once
  PostData(1, fDetEltTDHistList);  
  PostData(2, fDetEltTTHistList);  
  PostData(3, fChamberTDHistList);
  PostData(4, fChamberTTHistList);
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::UserExec(Option_t *)
{
  /// Main event loop
  
  // check the OCDB has been loaded properly
  if (!fOCDBLoaded) return;
  
  // get the current event
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!esd) return;
  
  // get the centrality
  fCurrentCentrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
  
  // loop over tracks
  AliMUONTrack track;
  Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
  {
    AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
    
    if(!esdTrack->ContainTrackerData()) continue;
    
    if(fMatchTrig && !esdTrack->ContainTriggerData()) continue;
    
    Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
    Double_t eta = esdTrack->Eta();
    if(fApplyAccCut && !(thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 9. && eta >= -4. && eta <= -2.5)) continue;
    
    AliMUONESDInterface::ESDToMUON(*esdTrack, track);
    
    TrackParamLoop(track.GetTrackParamAtCluster());
  }
  
  // post the output data:
  PostData(1, fDetEltTDHistList);  
  PostData(2, fDetEltTTHistList);  
  PostData(3, fChamberTDHistList);
  PostData(4, fChamberTTHistList);
}
    
//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::Terminate(Option_t *)
{
  /// final plots
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::TrackParamLoop(const TObjArray* trackParams)
{
  /// Loop on all the track params and fill the histos
  
  Bool_t trackFilter[10];
  memset(trackFilter, kFALSE, 10*sizeof(Bool_t));
  Bool_t chamberResponse[10];
  memset(chamberResponse, kFALSE, 10*sizeof(Bool_t));
  
  // check if the chamber responds
  Int_t nTrackParams = (Int_t) trackParams->GetEntriesFast();
  for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
  { 
    Int_t chamberId = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam))->GetClusterPtr()->GetChamberId();
    trackFilter[chamberId] = kTRUE;
    chamberResponse[chamberId] = kTRUE;
  }
  
  // To make sure the calculation of the efficiency of a given chamber (DE) is not biased by the tracking algorithm
  // we must make sure the track would have been reconstructed whatever this chamber (DE) has responded or not.
  // If the track is valid for a given chamber, the following code set trackFilter[chamberId] to kTRUE.
  for (Int_t station = 0; station < 4; ++station)
  {
    Int_t filter;
    Int_t ch1 = 2*station;
    Int_t ch2 = 2*station + 1;
    Int_t ch3 = 2*station + 2;
    Int_t ch4 = 2*station + 3;
    if (station < 3 )
    {
      filter           = trackFilter[ch1];
      trackFilter[ch1] = trackFilter[ch2];
      trackFilter[ch2] = filter;
    }
    else
    {
      if (chamberResponse[ch3] && chamberResponse[ch4])
      {
	filter           = trackFilter[ch1];
	trackFilter[ch1] = trackFilter[ch2];
	trackFilter[ch2] = filter;
      }
      else
      {
	trackFilter[ch1] = kFALSE;
	trackFilter[ch2] = kFALSE;
      }
      
      if (chamberResponse[ch1] && chamberResponse[ch2])
      {
	filter           = trackFilter[ch3];
	trackFilter[ch3] = trackFilter[ch4];
	trackFilter[ch4] = filter;
      }
      else
      {
	trackFilter[ch3] = kFALSE;
	trackFilter[ch4] = kFALSE;
      }
    }
  }
  
  // loop over track parameters
  Int_t oldChamber = -1;
  for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
  {
    AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam));
    AliMUONVCluster* cluster = trackParam->GetClusterPtr();
    
    Int_t newChamber = cluster->GetChamberId();
    
    // fill histograms if the track is valid for this chamber
    if(trackFilter[newChamber])
    {
      Int_t detElt = cluster->GetDetElemId();
      
      ///track position in the global coordinate system
      Double_t posXG = trackParam->GetNonBendingCoor(); 
      Double_t posYG = trackParam->GetBendingCoor(); 
      Double_t posZG = trackParam->GetZ(); 
      
      ///track position in the coordinate system of the DE
      Double_t posXL, posYL, posZL;
      fTransformer->Global2Local(detElt, posXG, posYG, posZG, posXL, posYL, posZL);
      
      // fill histograms of the cluster positions on the detection element of the TRACKS DETECTED (TD)
      FillTDHistos(newChamber, detElt, posXL, posYL);
      
      // fill histograms of the cluster positions on the detection element of ALL THE TRACKS (TT)
      FillTTHistos(newChamber, detElt, posXL, posYL);
    }
    
    // look for missing cluster(s) if any
    if (newChamber != oldChamber) 
    {
      if (newChamber > oldChamber + 1)
      {
	Int_t nbrMissChamber = newChamber - (oldChamber + 1);
	
	// find the DE(s) that should have been fired and fill the corresponding histograms
	FindAndFillMissedDetElt(trackParam, trackFilter, oldChamber+1, nbrMissChamber);
      }
      
      // in case the last chamber has not responded
      if ( iTrackParam == nTrackParams-1 && newChamber != 9) FindAndFillMissedDetElt(trackParam, trackFilter, 9, 1);
    }
    
    oldChamber = newChamber; 
  } 
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::FindAndFillMissedDetElt(const AliMUONTrackParam* trackParam,
							     const Bool_t* trackFilter,
							     Int_t firstMissCh, Int_t nbrMissCh)
{
  /// Find which detection elements should have been hit but were missed, and fill the TT histos appropriately
  
  // copy track parameters for extrapolation
  AliMUONTrackParam extrapTrackParam(*trackParam);
  
  // loop over missing chambers
  for (Int_t iCh = 0; iCh < nbrMissCh; ++iCh)
  {
    Int_t chamber = firstMissCh + iCh;
    
    // skip this chamber if the track is not valid for it
    if(!trackFilter[chamber]) continue;
    
    Int_t nbrOfDetElt =  AliMpDEManager::GetNofDEInChamber(chamber, kTRUE);
    
    Double_t pos1[6] = {0, 0, 0, 0, 0, 0};
    Double_t pos2[6] = {0, 0, 0, 0, 0, 0};
    Double_t posMiss[2] = {0, 0};
    
    // track position at the chamber z
    pos1[2] = AliMUONConstants::DefaultChamberZ(chamber);
    AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos1[2]);
    pos1[0] = extrapTrackParam.GetNonBendingCoor();
    pos1[1] = extrapTrackParam.GetBendingCoor();
    
    // track position at the chamber z + dz (where dz = distance between the 2 chamber in the station)
    pos2[2] = AliMUONConstants::DefaultChamberZ(chamber) + AliMUONConstants::DzCh();
    AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos2[2]);
    pos2[0] = extrapTrackParam.GetNonBendingCoor();
    pos2[1] = extrapTrackParam.GetBendingCoor();
    
    // loop over all the detection element of the chamber
    for (Int_t iDE = 0; iDE < nbrOfDetElt; iDE++)
    {
      Int_t deId = (chamber + 1)*fgkOffset + iDE;
      
      // track positions (at chamber z and chamber z + dz) in the local coordinate system of the DE
      fTransformer->Global2Local(deId, pos1[0], pos1[1], pos1[2], pos1[3], pos1[4], pos1[5]);
      fTransformer->Global2Local(deId, pos2[0], pos2[1], pos2[2], pos2[3], pos2[4], pos2[5]);
      
      // track position at z=0 in the local coordinate system of the DE
      CoordinatesOfMissingCluster(pos1[3], pos1[4], pos1[5], pos2[3], pos2[4], pos2[5], posMiss[0], posMiss[1]);
      
      // check if the track cross this DE and fill the corresponding histogram
      if (CoordinatesInDetElt(deId, posMiss[0], posMiss[1])) FillTTHistos(chamber, deId, posMiss[0], posMiss[1]);
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::CoordinatesOfMissingCluster(Double_t x1, Double_t y1, Double_t z1,
								 Double_t x2, Double_t y2, Double_t z2,
								 Double_t& x, Double_t& y) const
{
  /// Compute the coordinates of the missing cluster. They are defined by the intersection between
  /// the straigth line joining two extrapolated points (1 and 2) and the detection element plane.
  /// In the local coordinates, this means Z=0 in the parametric equation of the line.
  Double_t t = - z1 / (z2 - z1);
  x = t * (x2 - x1) + x1;
  y = t * (y2 - y1) + y1;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskMuonTrackingEff::CoordinatesInDetElt(Int_t DeId, Double_t x, Double_t y) const
{
  /// Return kTRUE if the coordinates are in the Detection Element.
  /// This is done by checking if a pad correspond to the (x, y) position.
  const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath0);
  const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath1);
  if (!seg1 || !seg2) return kFALSE;
  AliMpPad pad1 = seg1->PadByPosition(x, y, kFALSE);
  AliMpPad pad2 = seg2->PadByPosition(x, y, kFALSE);
  return (pad1.IsValid() && pad2.IsValid());
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::FillTDHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
{
  /// Fill the histo for detected tracks
  ((TH3F*) fDetEltTDHistList->At(FromDetElt2iDet(chamber, detElt)))->Fill(posXL, posYL, fCurrentCentrality);
  ((TH2F*) fChamberTDHistList->At(chamber))->Fill(FromDetElt2LocalId(chamber, detElt), fCurrentCentrality);
  ((TH2F*) fChamberTDHistList->At(10))->Fill(chamber+1, fCurrentCentrality);
}

//________________________________________________________________________
void AliAnalysisTaskMuonTrackingEff::FillTTHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
{
  /// Fill the histo for all tracks
  ((TH3F*) fDetEltTTHistList->At(FromDetElt2iDet(chamber, detElt))) -> Fill(posXL, posYL, fCurrentCentrality);
  ((TH2F*) fChamberTTHistList->At(chamber))->Fill(FromDetElt2LocalId(chamber, detElt), fCurrentCentrality);
  ((TH2F*) fChamberTTHistList->At(10))->Fill(chamber+1, fCurrentCentrality);
}

//________________________________________________________________________
Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2iDet(Int_t chamber, Int_t detElt) const
{
  /// Connexion between the detection element Id and its position in the list of histograms
  Int_t iDet = FromDetElt2LocalId(chamber, detElt);
  for (Int_t iCh = chamber-1; iCh >=0; iCh--) iDet += fgkNbrOfDetectionElt[iCh];
  return iDet;
}

//________________________________________________________________________
Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2LocalId(Int_t chamber, Int_t detElt) const
{
  /// Connexion between the detection element Id and its number in the chamber
  return detElt - fgkOffset*(chamber+1);    
}

//________________________________________________________________________
Int_t AliAnalysisTaskMuonTrackingEff::FromLocalId2DetElt(Int_t chamber, Int_t iDet) const
{
  /// Connexion between the number of the detection element in the chamber and its Id
  return iDet + fgkOffset*(chamber+1);    
}
Commit	Line	Data
128a8042	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
27de2dfb	16	/* $Id$ */
27de2dfb	17
128a8042	18	//Class to calculate the intrinsic efficiency of the detection elements of the
128a8042	19	//MUON tracking chambers in function of the position in the detection element.
79720f5b	20	//WOrk on ESD only
128a8042	21	//Author: Nicolas LE BRIS - SUBATECH Nantes
79720f5b	22	// Modified by Matthieu LENHARDT - SUBATECH Nantes
e1ec553d	23	// Modified by Antoine LARDEUX - SUBATECH Nantes
128a8042	24
e1ec553d	25	// ROOT includes
79720f5b	26	#include <TList.h>
e1ec553d	27	#include <TH3F.h>
128a8042	28	#include <TH2F.h>
e1ec553d	29	#include <TObjArray.h>
e1ec553d	30	#include <TGeoGlobalMagField.h>
128a8042	31
e1ec553d	32	// STEER includes
128a8042	33	#include "AliESDEvent.h"
128a8042	34	#include "AliESDMuonTrack.h"
fc7a3fd3	35	#include "AliGeomManager.h"
79720f5b	36	#include "AliCDBManager.h"
e1ec553d	37	#include "AliESDVZERO.h"
fc7a3fd3	38
e1ec553d	39	// ANALYSIS includes
	40	#include "AliInputEventHandler.h"
	41	#include "AliAnalysisManager.h"
128a8042	42	#include "AliAnalysisTaskMuonTrackingEff.h"
e1ec553d	43	#include "AliCentrality.h"
128a8042	44
128a8042	45	//MUON includes
e1ec553d	46	#include "AliMUONCDB.h"
e1ec553d	47	#include "AliMUONESDInterface.h"
27f15548	48	#include "AliMUONGeometryTransformer.h"
e1ec553d	49	#include "AliMUONTrack.h"
	50	#include "AliMUONTrackParam.h"
	51	#include "AliMUONTrackExtrap.h"
	52	#include "AliMUONVCluster.h"
79720f5b	53	#include "AliMUONConstants.h"
128a8042	54
e1ec553d	55	//include MUON/mapping:
	56	#include "AliMpDEManager.h"
	57	#include "AliMpSegmentation.h"
	58	#include "AliMpVSegmentation.h"
	59	#include "AliMpPad.h"
	60
128a8042	61	ClassImp(AliAnalysisTaskMuonTrackingEff)
128a8042	62
e1ec553d	63	const Int_t AliAnalysisTaskMuonTrackingEff::fgkNbrOfDetectionElt[10] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26};
e1ec553d	64	const Int_t AliAnalysisTaskMuonTrackingEff::fgkOffset = 100;
128a8042	65
128a8042	66	//________________________________________________________________________
e1ec553d	67	AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff() :
	68	AliAnalysisTaskSE(),
	69	fOCDBLoaded(kFALSE),
	70	fOCDBpath(""),
	71	fMatchTrig(kFALSE),
	72	fApplyAccCut(kFALSE),
	73	fCurrentCentrality(0.),
fc7a3fd3	74	fTransformer(0x0),
e1ec553d	75	fDetEltTDHistList(0x0),
	76	fDetEltTTHistList(0x0),
	77	fChamberTDHistList(0x0),
	78	fChamberTTHistList(0x0)
128a8042	79	{
e1ec553d	80	/// Default constructor
128a8042	81	}
128a8042	82
128a8042	83	//________________________________________________________________________
e1ec553d	84	AliAnalysisTaskMuonTrackingEff::AliAnalysisTaskMuonTrackingEff(TString name) :
	85	AliAnalysisTaskSE(name),
	86	fOCDBLoaded(kFALSE),
	87	fOCDBpath("raw://"),
	88	fMatchTrig(kFALSE),
	89	fApplyAccCut(kFALSE),
	90	fCurrentCentrality(100.),
fc7a3fd3	91	fTransformer(0x0),
e1ec553d	92	fDetEltTDHistList(0x0),
	93	fDetEltTTHistList(0x0),
	94	fChamberTDHistList(0x0),
	95	fChamberTTHistList(0x0)
128a8042	96	{
e1ec553d	97	/// Constructor
	98
	99	// Output slots 0 to 5 writes into a TClonesArray:
	100	DefineOutput(1, TList::Class());
	101	DefineOutput(2, TList::Class());
	102	DefineOutput(3, TList::Class());
	103	DefineOutput(4, TList::Class());
128a8042	104	}
128a8042	105
e1ec553d	106	//________________________________________________________________________
128a8042	107	AliAnalysisTaskMuonTrackingEff::~AliAnalysisTaskMuonTrackingEff()
128a8042	108	{
e1ec553d	109	/// Destructor
	110	if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
	111	delete fDetEltTDHistList;
	112	delete fDetEltTTHistList;
	113	delete fChamberTDHistList;
	114	delete fChamberTTHistList;
	115	}
fc7a3fd3	116	delete fTransformer;
128a8042	117	}
128a8042	118
128a8042	119	//________________________________________________________________________
e1ec553d	120	void AliAnalysisTaskMuonTrackingEff::NotifyRun()
128a8042	121	{
e1ec553d	122	/// Load the OCDB and the Geometry
	123
	124	// Load it only once
	125	if (fOCDBLoaded) return;
	126
	127	// OCDB
	128	AliCDBManager* man = AliCDBManager::Instance();
	129	if (man->IsDefaultStorageSet()) printf("EfficiencyTask: CDB default storage already set!");
	130	else man->SetDefaultStorage(fOCDBpath.Data());
	131	if (man->GetRun() > -1) printf("EfficiencyTask: run number already set!");
	132	else man->SetRun(fCurrentRunNumber);
	133
	134	// Geometry
	135	if (!AliGeomManager::GetGeometry()) {
	136	AliGeomManager::LoadGeometry();
	137	if (!AliGeomManager::GetGeometry()) return;
	138	if (!AliGeomManager::ApplyAlignObjsFromCDB("MUON")) return;
	139	}
	140	fTransformer = new AliMUONGeometryTransformer();
	141	fTransformer->LoadGeometryData();
	142
	143	// Magnetic field for track extrapolation
	144	if (!TGeoGlobalMagField::Instance()->GetField()) {
	145	if (!AliMUONCDB::LoadField()) return;
	146	}
	147
	148	// Mapping
	149	if (!AliMpSegmentation::Instance(kFALSE)) {
	150	if (!AliMUONCDB::LoadMapping(kTRUE)) return;
	151	}
	152
	153	// RecoParam for refitting
	154	AliMUONRecoParam* recoParam = AliMUONCDB::LoadRecoParam();
	155	if (!recoParam) return;
	156	AliMUONESDInterface::ResetTracker(recoParam);
	157
	158	fOCDBLoaded = kTRUE;
128a8042	159	}
128a8042	160
128a8042	161	//________________________________________________________________________
e1ec553d	162	void AliAnalysisTaskMuonTrackingEff::UserCreateOutputObjects()
128a8042	163	{
e1ec553d	164	/// Define efficiency histograms
	165
	166	fDetEltTDHistList = new TList();
	167	fDetEltTDHistList->SetOwner();
	168	fDetEltTTHistList = new TList();
	169	fDetEltTTHistList->SetOwner();
	170	fChamberTDHistList = new TList();
	171	fChamberTDHistList->SetOwner();
	172	fChamberTTHistList = new TList();
	173	fChamberTTHistList->SetOwner();
	174
	175	TH2F *h2;
	176	TH3F *h3;
	177	TString histNameTD;
	178	TString histNameTT;
	179	TString histTitle;
	180	Int_t nCentBins = 22;
	181	Double_t centRange[2] = {-5., 105.};
	182	Int_t iDEGlobal = 0;
	183
	184	for (Int_t iCh = 0; iCh < 10; iCh++)
	185	{
	186	// histograms per chamber
	187	histTitle.Form("ChamberNbr %d", iCh+1);
	188	histNameTD.Form("TD_ChamberNbr%d", iCh+1);
	189	histNameTT.Form("TT_ChamberNbr%d",iCh+1);
	190	h2 = new TH2F(histNameTD, histTitle, fgkNbrOfDetectionElt[iCh], 0., fgkNbrOfDetectionElt[iCh], nCentBins, centRange[0], centRange[1]);
	191	fChamberTDHistList->AddAt(h2, iCh);
	192	h2 = new TH2F(histNameTT, histTitle, fgkNbrOfDetectionElt[iCh], 0., fgkNbrOfDetectionElt[iCh], nCentBins, centRange[0], centRange[1]);
	193	fChamberTTHistList->AddAt(h2, iCh);
	194
	195	// histograms per DE
	196	for (Int_t iDE = 0; iDE < fgkNbrOfDetectionElt[iCh]; iDE++)
	197	{
	198	Int_t deId = FromLocalId2DetElt(iCh, iDE);
	199	histTitle.Form("detEltNbr %d",deId);
	200	histNameTD.Form("TD_detEltNbr%d",deId);
	201	histNameTT.Form("TT_detEltNbr%d",deId);
	202	if(iCh < 4)
	203	{// chambers 1 -> 4
	204	h3 = new TH3F(histNameTD, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
	205	fDetEltTDHistList->AddAt(h3, iDEGlobal);
	206	h3 = new TH3F(histNameTT, histTitle, 12, -10.0 , 110.0, 12, -10.0, 110.0, nCentBins, centRange[0], centRange[1]);
	207	fDetEltTTHistList->AddAt(h3, iDEGlobal);
	208	}
	209	else
	210	{// chambers 5 -> 10
	211	h3 = new TH3F(histNameTD, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
	212	fDetEltTDHistList->AddAt(h3, iDEGlobal);
	213	h3 = new TH3F(histNameTT, histTitle, 28, -140.0, 140.0, 8, -40.0, 40.0, nCentBins, centRange[0], centRange[1]);
	214	fDetEltTTHistList->AddAt(h3, iDEGlobal);
	215	}
	216	iDEGlobal++;
	217	}
	218	}
	219
	220	// global histograms per chamber
	221	h2 = new TH2F("TD_Chambers 11", "Chambers 11", 10, 0.5, 10.5, nCentBins, centRange[0], centRange[1]);
	222	fChamberTDHistList->AddAt(h2, 10);
	223	h2 = new TH2F("TT_Chambers 11", "Chambers 11", 10, 0.5, 10.5, nCentBins, centRange[0], centRange[1]);
	224	fChamberTTHistList->AddAt(h2, 10);
	225
	226	// post the output data at least once
	227	PostData(1, fDetEltTDHistList);
228	PostData(2, fDetEltTTHistList);
229	PostData(3, fChamberTDHistList);
230	PostData(4, fChamberTTHistList);
128a8042	231	}
128a8042	232
128a8042	233	//________________________________________________________________________
e1ec553d	234	void AliAnalysisTaskMuonTrackingEff::UserExec(Option_t *)
128a8042	235	{
e1ec553d	236	/// Main event loop
	237
	238	// check the OCDB has been loaded properly
	239	if (!fOCDBLoaded) return;
	240
	241	// get the current event
	242	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	243	if (!esd) return;
	244
	245	// get the centrality
	246	fCurrentCentrality = esd->GetCentrality()->GetCentralityPercentileUnchecked("V0M");
	247
	248	// loop over tracks
	249	AliMUONTrack track;
	250	Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks();
	251	for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
	252	{
	253	AliESDMuonTrack* esdTrack = esd->GetMuonTrack(iTrack);
	254
	255	if(!esdTrack->ContainTrackerData()) continue;
	256
	257	if(fMatchTrig && !esdTrack->ContainTriggerData()) continue;
	258
	259	Double_t thetaTrackAbsEnd = TMath::ATan(esdTrack->GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
	260	Double_t eta = esdTrack->Eta();
	261	if(fApplyAccCut && !(thetaTrackAbsEnd >= 2. && thetaTrackAbsEnd <= 9. && eta >= -4. && eta <= -2.5)) continue;
	262
	263	AliMUONESDInterface::ESDToMUON(*esdTrack, track);
	264
	265	TrackParamLoop(track.GetTrackParamAtCluster());
	266	}
	267
	268	// post the output data:
	269	PostData(1, fDetEltTDHistList);
	270	PostData(2, fDetEltTTHistList);
	271	PostData(3, fChamberTDHistList);
	272	PostData(4, fChamberTTHistList);
79720f5b	273	}
e1ec553d	274
79720f5b	275	//________________________________________________________________________
e1ec553d	276	void AliAnalysisTaskMuonTrackingEff::Terminate(Option_t *)
79720f5b	277	{
e1ec553d	278	/// final plots
79720f5b	279	}
0246246b	280
79720f5b	281	//________________________________________________________________________
e1ec553d	282	void AliAnalysisTaskMuonTrackingEff::TrackParamLoop(const TObjArray* trackParams)
79720f5b	283	{
e1ec553d	284	/// Loop on all the track params and fill the histos
79720f5b	285
e1ec553d	286	Bool_t trackFilter[10];
	287	memset(trackFilter, kFALSE, 10*sizeof(Bool_t));
	288	Bool_t chamberResponse[10];
	289	memset(chamberResponse, kFALSE, 10*sizeof(Bool_t));
	290
	291	// check if the chamber responds
	292	Int_t nTrackParams = (Int_t) trackParams->GetEntriesFast();
	293	for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
	294	{
	295	Int_t chamberId = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam))->GetClusterPtr()->GetChamberId();
	296	trackFilter[chamberId] = kTRUE;
	297	chamberResponse[chamberId] = kTRUE;
	298	}
79720f5b	299
e1ec553d	300	// To make sure the calculation of the efficiency of a given chamber (DE) is not biased by the tracking algorithm
	301	// we must make sure the track would have been reconstructed whatever this chamber (DE) has responded or not.
	302	// If the track is valid for a given chamber, the following code set trackFilter[chamberId] to kTRUE.
	303	for (Int_t station = 0; station < 4; ++station)
	304	{
	305	Int_t filter;
	306	Int_t ch1 = 2*station;
	307	Int_t ch2 = 2*station + 1;
	308	Int_t ch3 = 2*station + 2;
	309	Int_t ch4 = 2*station + 3;
	310	if (station < 3 )
79720f5b	311	{
e1ec553d	312	filter = trackFilter[ch1];
	313	trackFilter[ch1] = trackFilter[ch2];
	314	trackFilter[ch2] = filter;
	315	}
	316	else
	317	{
	318	if (chamberResponse[ch3] && chamberResponse[ch4])
	319	{
	320	filter = trackFilter[ch1];
	321	trackFilter[ch1] = trackFilter[ch2];
	322	trackFilter[ch2] = filter;
	323	}
	324	else
	325	{
	326	trackFilter[ch1] = kFALSE;
	327	trackFilter[ch2] = kFALSE;
	328	}
79720f5b	329
e1ec553d	330	if (chamberResponse[ch1] && chamberResponse[ch2])
	331	{
	332	filter = trackFilter[ch3];
	333	trackFilter[ch3] = trackFilter[ch4];
	334	trackFilter[ch4] = filter;
	335	}
	336	else
	337	{
	338	trackFilter[ch3] = kFALSE;
	339	trackFilter[ch4] = kFALSE;
	340	}
79720f5b	341	}
e1ec553d	342	}
	343
	344	// loop over track parameters
	345	Int_t oldChamber = -1;
	346	for (Int_t iTrackParam = 0; iTrackParam < nTrackParams; ++iTrackParam)
	347	{
	348	AliMUONTrackParam* trackParam = static_cast<AliMUONTrackParam*>(trackParams->UncheckedAt(iTrackParam));
	349	AliMUONVCluster* cluster = trackParam->GetClusterPtr();
	350
	351	Int_t newChamber = cluster->GetChamberId();
	352
	353	// fill histograms if the track is valid for this chamber
	354	if(trackFilter[newChamber])
79720f5b	355	{
e1ec553d	356	Int_t detElt = cluster->GetDetElemId();
79720f5b	357
e1ec553d	358	///track position in the global coordinate system
	359	Double_t posXG = trackParam->GetNonBendingCoor();
	360	Double_t posYG = trackParam->GetBendingCoor();
	361	Double_t posZG = trackParam->GetZ();
79720f5b	362
e1ec553d	363	///track position in the coordinate system of the DE
	364	Double_t posXL, posYL, posZL;
	365	fTransformer->Global2Local(detElt, posXG, posYG, posZG, posXL, posYL, posZL);
	366
	367	// fill histograms of the cluster positions on the detection element of the TRACKS DETECTED (TD)
	368	FillTDHistos(newChamber, detElt, posXL, posYL);
	369
	370	// fill histograms of the cluster positions on the detection element of ALL THE TRACKS (TT)
	371	FillTTHistos(newChamber, detElt, posXL, posYL);
	372	}
	373
	374	// look for missing cluster(s) if any
	375	if (newChamber != oldChamber)
	376	{
	377	if (newChamber > oldChamber + 1)
	378	{
	379	Int_t nbrMissChamber = newChamber - (oldChamber + 1);
	380
	381	// find the DE(s) that should have been fired and fill the corresponding histograms
	382	FindAndFillMissedDetElt(trackParam, trackFilter, oldChamber+1, nbrMissChamber);
	383	}
	384
	385	// in case the last chamber has not responded
	386	if ( iTrackParam == nTrackParams-1 && newChamber != 9) FindAndFillMissedDetElt(trackParam, trackFilter, 9, 1);
79720f5b	387	}
e1ec553d	388
	389	oldChamber = newChamber;
	390	}
	391	}
	392
	393	//________________________________________________________________________
	394	void AliAnalysisTaskMuonTrackingEff::FindAndFillMissedDetElt(const AliMUONTrackParam* trackParam,
	395	const Bool_t* trackFilter,
	396	Int_t firstMissCh, Int_t nbrMissCh)
	397	{
	398	/// Find which detection elements should have been hit but were missed, and fill the TT histos appropriately
79720f5b	399
e1ec553d	400	// copy track parameters for extrapolation
	401	AliMUONTrackParam extrapTrackParam(*trackParam);
	402
	403	// loop over missing chambers
	404	for (Int_t iCh = 0; iCh < nbrMissCh; ++iCh)
	405	{
	406	Int_t chamber = firstMissCh + iCh;
	407
	408	// skip this chamber if the track is not valid for it
	409	if(!trackFilter[chamber]) continue;
	410
	411	Int_t nbrOfDetElt = AliMpDEManager::GetNofDEInChamber(chamber, kTRUE);
	412
	413	Double_t pos1[6] = {0, 0, 0, 0, 0, 0};
	414	Double_t pos2[6] = {0, 0, 0, 0, 0, 0};
	415	Double_t posMiss[2] = {0, 0};
	416
	417	// track position at the chamber z
	418	pos1[2] = AliMUONConstants::DefaultChamberZ(chamber);
	419	AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos1[2]);
	420	pos1[0] = extrapTrackParam.GetNonBendingCoor();
	421	pos1[1] = extrapTrackParam.GetBendingCoor();
	422
	423	// track position at the chamber z + dz (where dz = distance between the 2 chamber in the station)
	424	pos2[2] = AliMUONConstants::DefaultChamberZ(chamber) + AliMUONConstants::DzCh();
	425	AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, pos2[2]);
	426	pos2[0] = extrapTrackParam.GetNonBendingCoor();
	427	pos2[1] = extrapTrackParam.GetBendingCoor();
	428
	429	// loop over all the detection element of the chamber
	430	for (Int_t iDE = 0; iDE < nbrOfDetElt; iDE++)
	431	{
	432	Int_t deId = (chamber + 1)*fgkOffset + iDE;
	433
	434	// track positions (at chamber z and chamber z + dz) in the local coordinate system of the DE
	435	fTransformer->Global2Local(deId, pos1[0], pos1[1], pos1[2], pos1[3], pos1[4], pos1[5]);
	436	fTransformer->Global2Local(deId, pos2[0], pos2[1], pos2[2], pos2[3], pos2[4], pos2[5]);
	437
	438	// track position at z=0 in the local coordinate system of the DE
	439	CoordinatesOfMissingCluster(pos1[3], pos1[4], pos1[5], pos2[3], pos2[4], pos2[5], posMiss[0], posMiss[1]);
	440
	441	// check if the track cross this DE and fill the corresponding histogram
	442	if (CoordinatesInDetElt(deId, posMiss[0], posMiss[1])) FillTTHistos(chamber, deId, posMiss[0], posMiss[1]);
	443	}
	444	}
128a8042	445	}
0246246b	446
e1ec553d	447	//________________________________________________________________________
	448	void AliAnalysisTaskMuonTrackingEff::CoordinatesOfMissingCluster(Double_t x1, Double_t y1, Double_t z1,
	449	Double_t x2, Double_t y2, Double_t z2,
	450	Double_t& x, Double_t& y) const
	451	{
	452	/// Compute the coordinates of the missing cluster. They are defined by the intersection between
	453	/// the straigth line joining two extrapolated points (1 and 2) and the detection element plane.
	454	/// In the local coordinates, this means Z=0 in the parametric equation of the line.
	455	Double_t t = - z1 / (z2 - z1);
	456	x = t * (x2 - x1) + x1;
	457	y = t * (y2 - y1) + y1;
	458	}
0246246b	459
0246246b	460	//________________________________________________________________________
e1ec553d	461	Bool_t AliAnalysisTaskMuonTrackingEff::CoordinatesInDetElt(Int_t DeId, Double_t x, Double_t y) const
0246246b	462	{
e1ec553d	463	/// Return kTRUE if the coordinates are in the Detection Element.
	464	/// This is done by checking if a pad correspond to the (x, y) position.
	465	const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath0);
	466	const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(DeId, AliMp::kCath1);
	467	if (!seg1 \|\| !seg2) return kFALSE;
	468	AliMpPad pad1 = seg1->PadByPosition(x, y, kFALSE);
	469	AliMpPad pad2 = seg2->PadByPosition(x, y, kFALSE);
	470	return (pad1.IsValid() && pad2.IsValid());
	471	}
0246246b	472
e1ec553d	473	//________________________________________________________________________
	474	void AliAnalysisTaskMuonTrackingEff::FillTDHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
	475	{
	476	/// Fill the histo for detected tracks
	477	((TH3F*) fDetEltTDHistList->At(FromDetElt2iDet(chamber, detElt)))->Fill(posXL, posYL, fCurrentCentrality);
	478	((TH2F*) fChamberTDHistList->At(chamber))->Fill(FromDetElt2LocalId(chamber, detElt), fCurrentCentrality);
	479	((TH2F*) fChamberTDHistList->At(10))->Fill(chamber+1, fCurrentCentrality);
	480	}
0246246b	481
e1ec553d	482	//________________________________________________________________________
	483	void AliAnalysisTaskMuonTrackingEff::FillTTHistos(Int_t chamber, Int_t detElt, Double_t posXL, Double_t posYL)
	484	{
	485	/// Fill the histo for all tracks
	486	((TH3F*) fDetEltTTHistList->At(FromDetElt2iDet(chamber, detElt))) -> Fill(posXL, posYL, fCurrentCentrality);
	487	((TH2F*) fChamberTTHistList->At(chamber))->Fill(FromDetElt2LocalId(chamber, detElt), fCurrentCentrality);
	488	((TH2F*) fChamberTTHistList->At(10))->Fill(chamber+1, fCurrentCentrality);
	489	}
0246246b	490
e1ec553d	491	//________________________________________________________________________
	492	Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2iDet(Int_t chamber, Int_t detElt) const
	493	{
	494	/// Connexion between the detection element Id and its position in the list of histograms
	495	Int_t iDet = FromDetElt2LocalId(chamber, detElt);
	496	for (Int_t iCh = chamber-1; iCh >=0; iCh--) iDet += fgkNbrOfDetectionElt[iCh];
	497	return iDet;
	498	}
79720f5b	499
e1ec553d	500	//________________________________________________________________________
	501	Int_t AliAnalysisTaskMuonTrackingEff::FromDetElt2LocalId(Int_t chamber, Int_t detElt) const
	502	{
	503	/// Connexion between the detection element Id and its number in the chamber
	504	return detElt - fgkOffset*(chamber+1);
	505	}
	506
	507	//________________________________________________________________________
	508	Int_t AliAnalysisTaskMuonTrackingEff::FromLocalId2DetElt(Int_t chamber, Int_t iDet) const
	509	{
	510	/// Connexion between the number of the detection element in the chamber and its Id
	511	return iDet + fgkOffset*(chamber+1);
0246246b	512	}
79720f5b	513