[u/mrichter/AliRoot.git] / PWG0 / multPbPb / AliAnalysisTaskMultPbTracks.cxx

// AliAnalysisTaskMultPbTracks

// Author: Michele Floris, CERN
// TODO:
// - Add chi2/cluster plot for primary, secondaries and fakes
// FIXME:
// - re-implement centrality estimator to cut on tracks and multiplicity

#include "AliAnalysisTaskMultPbTracks.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisMultPbTrackHistoManager.h"
#include "AliAnalysisManager.h"
#include "AliESDtrackCuts.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "TH1I.h"
#include "TH3D.h"
#include "AliMCParticle.h"
#include "AliGenEventHeader.h"
#include "AliESDCentrality.h"

#include <iostream>

using namespace std;

ClassImp(AliAnalysisTaskMultPbTracks)

AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
: AliAnalysisTaskSE("TaskMultPbTracks"),
  fESD(0),fHistoManager(0),fCentrBin(0),fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
{
  // constructor

  DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
  DefineOutput(2, AliESDtrackCuts::Class());
  //  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());

  fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
  if(fIsMC) fHistoManager->SetSuffix("MC");

}
AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
  : AliAnalysisTaskSE(name),
    fESD(0),fHistoManager(0),fCentrBin(0),fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
{
  //
  // Standard constructur which should be used
  //

  DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
  DefineOutput(2, AliESDtrackCuts::Class());
  //  DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());

  fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
  if(fIsMC) fHistoManager->SetSuffix("MC");

}

AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) : 
  AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrBin(0), fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
{
  //copy ctor
  fESD = obj.fESD ;
  fHistoManager= obj.fHistoManager;
  fCentrBin = obj.fCentrBin;
  fCentralityEstimator = obj.fCentralityEstimator;
  fTrackCuts  = obj.fTrackCuts;
  fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;

}

AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
  // destructor

  if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
    if(fHistoManager) {
      delete fHistoManager;
      fHistoManager = 0;
    }
  }
  // Histo list should not be destroyed: fListWrapper is owner!

}
void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
{
  // Called once

  // For the DCA distributions, we want to use exactly the same cuts
  // as for the other distributions, with the exception of the DCA cut
  fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
  // Reset all DCA cuts; FIXME: is this all?
  fTrackCutsNoDCA->SetMaxDCAToVertexXY();
  fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
  fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
  fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();

}


void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
{
  // User code

  /* PostData(0) is taken care of by AliAnalysisTaskSE */
  PostData(1,fHistoManager);
  PostData(2,fTrackCuts);
  //  PostData(3,fCentralityEstimator);

  // Cache histogram pointers
  static TH3D * hTracks  [AliAnalysisMultPbTrackHistoManager::kNHistos];
  static TH1D * hDCA     [AliAnalysisMultPbTrackHistoManager::kNHistos];
  static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);

  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen       );
  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec       );
  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim   );
  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]    = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake   );
  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]  = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);

  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen       );

  fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
  if (strcmp(fESD->ClassName(),"AliESDEvent")) {
    AliFatal("Not processing ESDs");
  }
  
  // FIXME: use physics selection here to keep track of events lost?
  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);

  // Centrality selection
  // Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD);  
  // if(!isCentralitySelected) return;

  AliESDCentrality *centrality = fESD->GetCentrality();
  if(!centrality) {
    AliFatal("Centrality object not available"); 
  }
  else {
    Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
    
    if (centrBin != fCentrBin && fCentrBin != -1) return;
  }

  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);

  if (fIsMC) {
    

    if (!fMCEvent) {
      AliError("No MC info found");
    } else {
      
      //loop on the MC event, only over primaries, which are always
      //      the first in stack.
      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
      Int_t nPhysicalPrimaries = 0;
      for (Int_t ipart=0; ipart<nMCTracks; ipart++) { 
	
	AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
	
	// We don't care about neutrals and non-physical primaries
	if(mcPart->Charge() == 0) continue;

	//check if current particle is a physical primary
	if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
 
	nPhysicalPrimaries++;
	// Fill species histo
	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());

	
	// Get MC vertex
	//FIXME: which vertex do I take for MC?
	TArrayF   vertex;
	fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
	Float_t zv = vertex[2];
	//	Float_t zv = vtxESD->GetZ();
	// Fill generated histo
	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
	
      }
      hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);

    }
  }
  

  // FIXME: shall I take the primary vertex?
  const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
  if(!vtxESD) return;
  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);


  // loop on tracks
  Int_t ntracks = fESD->GetNumberOfTracks();


  for (Int_t itrack = 0; itrack<ntracks; itrack++) {    
    AliESDtrack * esdTrack = fESD->GetTrack(itrack);

    // Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
    Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
    Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);

    // Compute weighted offset
    // TODO: other choiches for the DCA?
    Double_t b = fESD->GetMagneticField();
    Double_t dca[2];
    Double_t cov[3];
    Double_t weightedDCA = 10;
    

    if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
      Double_t det = cov[0]*cov[2]-cov[1]*cov[1]; 
      if (det<=0) {
	AliError("DCA Covariance matrix is not positive definite!");
      }
      else {
	weightedDCA = (dca[0]*dca[0]*cov[2]+dca[1]*dca[1]*cov[0]-2*dca[0]*dca[1]*cov[1])/det; 
	weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
      }
      //      printf("dR:%e dZ%e  sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
    }
    

    // Alternative: p*DCA
    // Float_t xz[2]; 
    // esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
    // Float_t dca = xz[0]*esdTrack->P();

    // for each track
    if(accepted) 
      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
    if(acceptedNoDCA)
      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA);

    // Fill separately primaries and secondaries
    // FIXME: fakes? Is this the correct way to account for them?
    // Get label and corresponding mcPart;
    if (fIsMC) {
      //      Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
      Int_t label = esdTrack->GetLabel(); // 
      AliMCParticle *mcPart  = label < 0 ? 0 : (AliMCParticle*)fMCEvent->GetTrack(label);
      if (!mcPart)  {
	if(accepted)
	  hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	if(acceptedNoDCA)
	  hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
      }
      else {
	if(IsPhysicalPrimaryAndTransportBit(label)) {
	  // Fill species histo
	  fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
	  if(accepted)
	    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	  if(acceptedNoDCA)
	    hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA);
	} 
	else {
	  Int_t mfl=0;
	  Int_t indexMoth=mcPart->Particle()->GetFirstMother();
	  if(indexMoth>=0){
	    TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
	    Float_t codemoth = TMath::Abs(moth->GetPdgCode());
	    mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
	  }
	  if(mfl==3){ // strangeness
	    fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
	    if(accepted)
	      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	    if(acceptedNoDCA)
	      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);	  
	  }else{ // material
	    fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
	    if(accepted)
	      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	    if(acceptedNoDCA)
	      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA);	  

	  }


	}
      }
    }


  }


}

void   AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
  // terminate

}


Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {

  Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
  if (!physprim) return kFALSE;
  Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
  if(!transported) return kFALSE;

  return kTRUE;

}
Commit	Line	Data
a23f7c97	1	// AliAnalysisTaskMultPbTracks
	2
	3	// Author: Michele Floris, CERN
	4	// TODO:
	5	// - Add chi2/cluster plot for primary, secondaries and fakes
e0376287	6	// FIXME:
e0376287	7	// - re-implement centrality estimator to cut on tracks and multiplicity
a23f7c97	8
	9	#include "AliAnalysisTaskMultPbTracks.h"
	10	#include "AliESDInputHandler.h"
	11	#include "AliAnalysisMultPbTrackHistoManager.h"
	12	#include "AliAnalysisManager.h"
	13	#include "AliESDtrackCuts.h"
a23f7c97	14	#include "AliMCEvent.h"
	15	#include "AliStack.h"
	16	#include "TH1I.h"
	17	#include "TH3D.h"
	18	#include "AliMCParticle.h"
	19	#include "AliGenEventHeader.h"
54b31d6a	20	#include "AliESDCentrality.h"
a23f7c97	21
	22	#include <iostream>
	23
	24	using namespace std;
	25
	26	ClassImp(AliAnalysisTaskMultPbTracks)
	27
	28	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
54b31d6a	29	: AliAnalysisTaskSE("TaskMultPbTracks"),
54b31d6a	30	fESD(0),fHistoManager(0),fCentrBin(0),fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	31	{
	32	// constructor
	33
	34	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	35	DefineOutput(2, AliESDtrackCuts::Class());
e0376287	36	// DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
a23f7c97	37
	38	fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
	39	if(fIsMC) fHistoManager->SetSuffix("MC");
	40
	41	}
	42	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
54b31d6a	43	: AliAnalysisTaskSE(name),
54b31d6a	44	fESD(0),fHistoManager(0),fCentrBin(0),fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	45	{
	46	//
	47	// Standard constructur which should be used
	48	//
	49
	50	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	51	DefineOutput(2, AliESDtrackCuts::Class());
e0376287	52	// DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
a23f7c97	53
	54	fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
	55	if(fIsMC) fHistoManager->SetSuffix("MC");
	56
	57	}
	58
	59	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) :
54b31d6a	60	AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrBin(0), fCentralityEstimator(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	61	{
	62	//copy ctor
	63	fESD = obj.fESD ;
	64	fHistoManager= obj.fHistoManager;
54b31d6a	65	fCentrBin = obj.fCentrBin;
54b31d6a	66	fCentralityEstimator = obj.fCentralityEstimator;
a23f7c97	67	fTrackCuts = obj.fTrackCuts;
	68	fTrackCutsNoDCA = obj.fTrackCutsNoDCA;
	69
	70	}
	71
	72	AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
	73	// destructor
	74
	75	if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
	76	if(fHistoManager) {
	77	delete fHistoManager;
	78	fHistoManager = 0;
	79	}
	80	}
	81	// Histo list should not be destroyed: fListWrapper is owner!
	82
	83	}
	84	void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
	85	{
	86	// Called once
	87
	88	// For the DCA distributions, we want to use exactly the same cuts
	89	// as for the other distributions, with the exception of the DCA cut
	90	fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
	91	// Reset all DCA cuts; FIXME: is this all?
	92	fTrackCutsNoDCA->SetMaxDCAToVertexXY();
	93	fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
	94	fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
	95	fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();
	96
	97	}
	98
	99
	100	void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
	101	{
	102	// User code
	103
	104	/* PostData(0) is taken care of by AliAnalysisTaskSE */
	105	PostData(1,fHistoManager);
54b31d6a	106	PostData(2,fTrackCuts);
e0376287	107	// PostData(3,fCentralityEstimator);
a23f7c97	108
54b31d6a	109	// Cache histogram pointers
e0376287	110	static TH3D * hTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
	111	static TH1D * hDCA [AliAnalysisMultPbTrackHistoManager::kNHistos];
	112	static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
a23f7c97	113	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen );
	114	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec );
	115	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
	116	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
	117	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
	118	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
	119
	120	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen );
	121	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
	122	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
	123	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
	124	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
	125	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
	126
e0376287	127	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
a23f7c97	128
	129	fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
	130	if (strcmp(fESD->ClassName(),"AliESDEvent")) {
	131	AliFatal("Not processing ESDs");
	132	}
	133
	134	// FIXME: use physics selection here to keep track of events lost?
	135	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
	136
54b31d6a	137	// Centrality selection
e0376287	138	// Bool_t isCentralitySelected = fCentrSelector->IsSelected(fESD);
	139	// if(!isCentralitySelected) return;
	140
54b31d6a	141	AliESDCentrality *centrality = fESD->GetCentrality();
54b31d6a	142	if(!centrality) {
e0376287	143	AliFatal("Centrality object not available");
54b31d6a	144	}
	145	else {
	146	Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
54b31d6a	147
a28a49f6	148	if (centrBin != fCentrBin && fCentrBin != -1) return;
54b31d6a	149	}
a23f7c97	150
	151	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
	152
	153	if (fIsMC) {
	154
	155
	156	if (!fMCEvent) {
	157	AliError("No MC info found");
	158	} else {
	159
e0376287	160	//loop on the MC event, only over primaries, which are always
	161	// the first in stack.
	162	Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
	163	Int_t nPhysicalPrimaries = 0;
	164	for (Int_t ipart=0; ipart<nMCTracks; ipart++) {
a23f7c97	165
	166	AliMCParticle mcPart = (AliMCParticle)fMCEvent->GetTrack(ipart);
	167
	168	// We don't care about neutrals and non-physical primaries
	169	if(mcPart->Charge() == 0) continue;
	170
54b31d6a	171	//check if current particle is a physical primary
e0376287	172	if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
54b31d6a	173
e0376287	174	nPhysicalPrimaries++;
	175	// Fill species histo
	176	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
	177
	178
a23f7c97	179	// Get MC vertex
	180	//FIXME: which vertex do I take for MC?
	181	TArrayF vertex;
	182	fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
	183	Float_t zv = vertex[2];
	184	// Float_t zv = vtxESD->GetZ();
	185	// Fill generated histo
	186	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
	187
	188	}
e0376287	189	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
e0376287	190
a23f7c97	191	}
	192	}
	193
e0376287	194
a23f7c97	195	// FIXME: shall I take the primary vertex?
	196	const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
	197	if(!vtxESD) return;
	198	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
	199
	200
	201	// loop on tracks
	202	Int_t ntracks = fESD->GetNumberOfTracks();
	203
	204
	205	for (Int_t itrack = 0; itrack<ntracks; itrack++) {
	206	AliESDtrack * esdTrack = fESD->GetTrack(itrack);
	207
	208	// Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
	209	Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
	210	Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
	211
	212	// Compute weighted offset
	213	// TODO: other choiches for the DCA?
	214	Double_t b = fESD->GetMagneticField();
	215	Double_t dca[2];
	216	Double_t cov[3];
	217	Double_t weightedDCA = 10;
	218
	219
	220
	221	if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
	222	Double_t det = cov[0]cov[2]-cov[1]cov[1];
	223	if (det<=0) {
	224	AliError("DCA Covariance matrix is not positive definite!");
	225	}
	226	else {
	227	weightedDCA = (dca[0]dca[0]cov[2]+dca[1]dca[1]cov[0]-2dca[0]dca[1]*cov[1])/det;
	228	weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
	229	}
	230	// printf("dR:%e dZ%e sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
	231	}
	232
	233
	234
	235
	236	// Alternative: p*DCA
	237	// Float_t xz[2];
	238	// esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
	239	// Float_t dca = xz[0]*esdTrack->P();
	240
	241	// for each track
	242	if(accepted)
	243	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	244	if(acceptedNoDCA)
	245	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA);
	246
	247	// Fill separately primaries and secondaries
	248	// FIXME: fakes? Is this the correct way to account for them?
	249	// Get label and corresponding mcPart;
	250	if (fIsMC) {
	251	// Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
	252	Int_t label = esdTrack->GetLabel(); //
e0376287	253	AliMCParticle mcPart = label < 0 ? 0 : (AliMCParticle)fMCEvent->GetTrack(label);
a23f7c97	254	if (!mcPart) {
	255	if(accepted)
	256	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	257	if(acceptedNoDCA)
	258	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
	259	}
	260	else {
e0376287	261	if(IsPhysicalPrimaryAndTransportBit(label)) {
	262	// Fill species histo
	263	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	264	if(accepted)
	265	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	266	if(acceptedNoDCA)
	267	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA);
	268	}
	269	else {
	270	Int_t mfl=0;
	271	Int_t indexMoth=mcPart->Particle()->GetFirstMother();
	272	if(indexMoth>=0){
	273	TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
	274	Float_t codemoth = TMath::Abs(moth->GetPdgCode());
	275	mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
	276	}
	277	if(mfl==3){ // strangeness
e0376287	278	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	279	if(accepted)
	280	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	281	if(acceptedNoDCA)
	282	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);
	283	}else{ // material
e0376287	284	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	285	if(accepted)
	286	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	287	if(acceptedNoDCA)
	288	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA);
	289
	290	}
	291
	292
	293	}
	294	}
	295	}
	296
	297
	298	}
	299
	300
	301	}
	302
	303	void AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
	304	// terminate
	305
	306	}
	307
	308
e0376287	309	Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
	310
	311	Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
	312	if (!physprim) return kFALSE;
	313	Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
	314	if(!transported) return kFALSE;
	315
	316	return kTRUE;
	317
	318	}