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

// AliAnalysisTaskMultPbTracks

// Author: Michele Floris, CERN
// TODO:
// - Add chi2/cluster plot for primary, secondaries and fakes


#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(2, TH1I::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());

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

  // Cache histogram pointers
  static TH3D * hTracks[AliAnalysisMultPbTrackHistoManager::kNHistos];
  static TH1D * hDCA   [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);


  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
  AliESDCentrality *centrality = fESD->GetCentrality();
  if(!centrality) {
    AliError("Centrality object not available"); // FIXME AliFatal here after debug
  }
  else {
    Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
    cout << fCentralityEstimator.Data() << " BIN: " << centrBin << endl;
    
    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
      //      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
      Int_t offset    = fMCEvent->GetPrimaryOffset();
      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries()+offset;
      for (Int_t ipart=offset; ipart<nMCTracks; ipart++) { 
	
	AliMCParticle *mcPart  = (AliMCParticle*)fMCEvent->GetTrack(ipart);
	
	// We don't care about neutrals and non-physical primaries
	if(mcPart->Charge() == 0) continue;

	// FIXME: add kTransportBit (uncomment below)
	if(!fMCEvent->Stack()->IsPhysicalPrimary(ipart)) continue;

	//check if current particle is a physical primary
	// Bool_t physprim=fMCEvent->IsPhysicalPrimary(label);
	// if (!physprim) continue;
	// if (!track) return kFALSE;
	// Bool_t transported = mcPart->Particle()->TestBit(kTransportBit);
	// if(!transported) return kFALSE;
 
	// 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);
	
      }
    }
  }
  
  // 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  = (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(fMCEvent->Stack()->IsPhysicalPrimary(label)) {
	  // FIXME add kTransportBit
	  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
	    if(accepted)
	      hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	    if(acceptedNoDCA)
	      hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA);	  
	  }else{ // material
	    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

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