[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>
#include "AliAnalysisMultPbCentralitySelector.h"

using namespace std;

ClassImp(AliAnalysisTaskMultPbTracks)

AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
: AliAnalysisTaskSE("TaskMultPbTracks"),
  fESD(0),fHistoManager(0),fCentrSelector(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),fCentrSelector(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), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
{
  //copy ctor
  fESD = obj.fESD ;
  fHistoManager= obj.fHistoManager;
  fTrackCuts  = obj.fTrackCuts;
  fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
  fCentrSelector = obj.fCentrSelector;

}

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,fCentrSelector);

  // 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 );
  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]        = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );

  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->IsCentralityBinSelected(fESD,fTrackCuts);  
  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();
  Int_t acceptedTracks = 0;

  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);

    if(accepted) acceptedTracks++;

    // 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);	  

	  }


	}
      }
    }


  }
  //  cout << acceptedTracks << endl;
  
  hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]  ->Fill(acceptedTracks);


}

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
a23f7c97	22	#include <iostream>
2bbfd8c6	23	#include "AliAnalysisMultPbCentralitySelector.h"
a23f7c97	24
	25	using namespace std;
	26
	27	ClassImp(AliAnalysisTaskMultPbTracks)
	28
	29	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
54b31d6a	30	: AliAnalysisTaskSE("TaskMultPbTracks"),
2bbfd8c6	31	fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	32	{
	33	// constructor
	34
	35	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	36	DefineOutput(2, AliESDtrackCuts::Class());
2bbfd8c6	37	DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
a23f7c97	38
	39	fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
	40	if(fIsMC) fHistoManager->SetSuffix("MC");
	41
	42	}
	43	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const char * name)
54b31d6a	44	: AliAnalysisTaskSE(name),
2bbfd8c6	45	fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	46	{
	47	//
	48	// Standard constructur which should be used
	49	//
	50
	51	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	52	DefineOutput(2, AliESDtrackCuts::Class());
2bbfd8c6	53	DefineOutput(3, AliAnalysisMultPbCentralitySelector::Class());
a23f7c97	54
	55	fHistoManager = new AliAnalysisMultPbTrackHistoManager("histoManager","Hitograms, Multiplicity, Track analysis");
	56	if(fIsMC) fHistoManager->SetSuffix("MC");
	57
	58	}
	59
	60	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks(const AliAnalysisTaskMultPbTracks& obj) :
2bbfd8c6	61	AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0)
a23f7c97	62	{
	63	//copy ctor
	64	fESD = obj.fESD ;
	65	fHistoManager= obj.fHistoManager;
a23f7c97	66	fTrackCuts = obj.fTrackCuts;
a23f7c97	67	fTrackCutsNoDCA = obj.fTrackCutsNoDCA;
2bbfd8c6	68	fCentrSelector = obj.fCentrSelector;
a23f7c97	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);
2bbfd8c6	107	PostData(3,fCentrSelector);
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
2bbfd8c6	127	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
2bbfd8c6	128	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
a23f7c97	129
	130	fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
	131	if (strcmp(fESD->ClassName(),"AliESDEvent")) {
	132	AliFatal("Not processing ESDs");
	133	}
	134
	135	// FIXME: use physics selection here to keep track of events lost?
	136	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
	137
54b31d6a	138	// Centrality selection
2bbfd8c6	139	Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,fTrackCuts);
	140	if(!isCentralitySelected) return;
	141
	142	// AliESDCentrality *centrality = fESD->GetCentrality();
	143	// if(!centrality) {
	144	// AliFatal("Centrality object not available");
	145	// }
	146	// else {
	147	// Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
54b31d6a	148
2bbfd8c6	149	// if (centrBin != fCentrBin && fCentrBin != -1) return;
2bbfd8c6	150	// }
a23f7c97	151
	152	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
	153
	154	if (fIsMC) {
	155
	156
	157	if (!fMCEvent) {
	158	AliError("No MC info found");
	159	} else {
	160
e0376287	161	//loop on the MC event, only over primaries, which are always
	162	// the first in stack.
	163	Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
	164	Int_t nPhysicalPrimaries = 0;
	165	for (Int_t ipart=0; ipart<nMCTracks; ipart++) {
a23f7c97	166
	167	AliMCParticle mcPart = (AliMCParticle)fMCEvent->GetTrack(ipart);
	168
	169	// We don't care about neutrals and non-physical primaries
	170	if(mcPart->Charge() == 0) continue;
	171
54b31d6a	172	//check if current particle is a physical primary
e0376287	173	if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
54b31d6a	174
e0376287	175	nPhysicalPrimaries++;
	176	// Fill species histo
	177	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
	178
	179
a23f7c97	180	// Get MC vertex
	181	//FIXME: which vertex do I take for MC?
	182	TArrayF vertex;
	183	fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
	184	Float_t zv = vertex[2];
	185	// Float_t zv = vtxESD->GetZ();
	186	// Fill generated histo
	187	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zv);
	188
	189	}
e0376287	190	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
e0376287	191
a23f7c97	192	}
	193	}
	194
e0376287	195
a23f7c97	196	// FIXME: shall I take the primary vertex?
	197	const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
	198	if(!vtxESD) return;
	199	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
	200
	201
	202	// loop on tracks
	203	Int_t ntracks = fESD->GetNumberOfTracks();
2bbfd8c6	204	Int_t acceptedTracks = 0;
a23f7c97	205
	206	for (Int_t itrack = 0; itrack<ntracks; itrack++) {
	207	AliESDtrack * esdTrack = fESD->GetTrack(itrack);
	208
	209	// Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
	210	Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
	211	Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
	212
2bbfd8c6	213	if(accepted) acceptedTracks++;
2bbfd8c6	214
a23f7c97	215	// Compute weighted offset
	216	// TODO: other choiches for the DCA?
	217	Double_t b = fESD->GetMagneticField();
	218	Double_t dca[2];
	219	Double_t cov[3];
	220	Double_t weightedDCA = 10;
	221
	222
	223
	224	if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
	225	Double_t det = cov[0]cov[2]-cov[1]cov[1];
	226	if (det<=0) {
	227	AliError("DCA Covariance matrix is not positive definite!");
	228	}
	229	else {
	230	weightedDCA = (dca[0]dca[0]cov[2]+dca[1]dca[1]cov[0]-2dca[0]dca[1]*cov[1])/det;
	231	weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
	232	}
	233	// printf("dR:%e dZ%e sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
	234	}
	235
	236
	237
	238
	239	// Alternative: p*DCA
	240	// Float_t xz[2];
	241	// esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
	242	// Float_t dca = xz[0]*esdTrack->P();
	243
	244	// for each track
	245	if(accepted)
	246	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	247	if(acceptedNoDCA)
	248	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA);
	249
	250	// Fill separately primaries and secondaries
	251	// FIXME: fakes? Is this the correct way to account for them?
	252	// Get label and corresponding mcPart;
	253	if (fIsMC) {
	254	// Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
	255	Int_t label = esdTrack->GetLabel(); //
e0376287	256	AliMCParticle mcPart = label < 0 ? 0 : (AliMCParticle)fMCEvent->GetTrack(label);
a23f7c97	257	if (!mcPart) {
	258	if(accepted)
	259	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	260	if(acceptedNoDCA)
	261	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA);
	262	}
	263	else {
e0376287	264	if(IsPhysicalPrimaryAndTransportBit(label)) {
	265	// Fill species histo
	266	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	267	if(accepted)
	268	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	269	if(acceptedNoDCA)
	270	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA);
	271	}
	272	else {
	273	Int_t mfl=0;
	274	Int_t indexMoth=mcPart->Particle()->GetFirstMother();
	275	if(indexMoth>=0){
	276	TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
	277	Float_t codemoth = TMath::Abs(moth->GetPdgCode());
	278	mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
	279	}
	280	if(mfl==3){ // strangeness
e0376287	281	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	282	if(accepted)
	283	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	284	if(acceptedNoDCA)
	285	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);
	286	}else{ // material
e0376287	287	fHistoManager->GetHistoSpecies(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	288	if(accepted)
	289	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	290	if(acceptedNoDCA)
	291	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA);
	292
	293	}
	294
	295
	296	}
	297	}
	298	}
	299
	300
	301	}
2bbfd8c6	302	// cout << acceptedTracks << endl;
	303
	304	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] ->Fill(acceptedTracks);
a23f7c97	305
	306
	307	}
	308
	309	void AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
	310	// terminate
	311
	312	}
	313
	314
e0376287	315	Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
	316
	317	Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
	318	if (!physprim) return kFALSE;
	319	Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
	320	if(!transported) return kFALSE;
	321
	322	return kTRUE;
	323
	324	}