[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 "TH2D.h"
#include "AliMCParticle.h"
#include "AliGenEventHeader.h"
#include "AliESDCentrality.h"
#include "AliMultiplicity.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),fIsTPCOnly(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),fIsTPCOnly(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),fIsTPCOnly(0)
{
  //copy ctor
  fESD = obj.fESD ;
  fHistoManager= obj.fHistoManager;
  fTrackCuts  = obj.fTrackCuts;
  fTrackCutsNoDCA  = obj.fTrackCutsNoDCA;
  fCentrSelector = obj.fCentrSelector;
  fIsMC = obj.fIsMC;
  fIsTPCOnly = obj.fIsTPCOnly;

}

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 TH2D * 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");
  }
  
  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);

  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fOfflineTrigger);

  if(!isSelected) return;
  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatPhysSel);


  if (fIsMC) {
    // Get MC vertex
    //FIXME: which vertex do I take for MC?
    TArrayF   vertex;
    fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
    Float_t zvGen = vertex[2];    

    if (!fMCEvent) {
      AliError("No MC info found");
    } else {
      
      //loop on the MC event, only over primaries, which are always
      //      the first in stack.
      // WRONG: D&B decays are produced in the transport... Need To Loop over all particles
      //      Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
      Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
      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;
	if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) {
	  // This cures a bug in Hijing
	  // A little hack here: I put those in the underflow bin of the process type to keep track of them
	  fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(-1);
	  continue;
	}

	nPhysicalPrimaries++;
	// Fill species histo and particle species
	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
	fHistoManager->FillParticleID(AliAnalysisMultPbTrackHistoManager::kHistoGen, mcPart);
	
	//	Float_t zv = vtxESD->GetZ();
	// Fill generated histo
	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
	Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
	cout << "Part " << partCode << endl;
	if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
	    partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
	    partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
	    partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
	    partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
	    partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar  
	    ){
	  cout << "Found " << partCode << endl;
	  
	  fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
	}
      }
      hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
      fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(zvGen);
  
    
    }
  }
  

  const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
  if(!vtxESD) return;
  // TODO: leave the cuts here or find a better place?)
  // Quality cut on vertexer Z, as suggested by Francesco Prino
  if(vtxESD->IsFromVertexerZ()) {
    if (vtxESD->GetNContributors() <= 0) return;
    if (vtxESD->GetDispersion() >= 0.04) return;
    if (vtxESD->GetZRes() >= 0.25) return;
  }
  // "Beam gas" vertex cut
  const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC(); 
  if(vtxESDTPC->GetNContributors()<1) return;
  if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0)))     return;

  // Fill vertex and statistics
  fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
  fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(vtxESD->GetZ());

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

  for (Int_t itrack = 0; itrack<ntracks; itrack++) {    
    AliESDtrack * esdTrack = fIsTPCOnly ? fTrackCuts->GetTPCOnlyTrack(fESD,itrack) : fESD->GetTrack(itrack);// FIXMEL TrackCuts or TrackCutsNoDCA for the TPC?
    if (!esdTrack) continue;
          
    // 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 weightedDCA = 10;
    

#if defined WEIGHTED_DCA
    Double_t b = fESD->GetMagneticField();
    Double_t dca[2];
    Double_t cov[3];
    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]);
    }
    
#elif defined TRANSVERSE_DCA
    Float_t xz[2]; 
    esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz); 
    weightedDCA = xz[0];
#endif
    

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

  // This cures a bug in Hijing (displaced primaries)
  if (fIsMC) {
    Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!    
    if (IsPhysicalPrimaryAndTransportBit(label)) {
      AliMCParticle *mcPart  =  (AliMCParticle*)fMCEvent->GetTrack(label);
      if(!mcPart) AliFatal("No particle");
      TArrayF   vertex;
      fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
      Float_t zvGen = vertex[2];    
      if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) continue;    
    }
  }
  // for each track
  if(accepted) 
    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
  if(acceptedNoDCA)
    hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());

    // 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,esdTrack->Pt());
      }
      else {
	if(IsPhysicalPrimaryAndTransportBit(label)) {
	  // Fill species histo
	  fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
	  if(accepted) {
	    hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	    Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
	    if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus  || 
		partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus ||
		partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus   || 
		partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus  ||
		partCode == AliAnalysisMultPbTrackHistoManager::kPartP       || 
		partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar  
		)
	      fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, partCode);
	  }
	  if(acceptedNoDCA)
	    hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA,esdTrack->Pt());
	} 
	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->GetHistoProcess(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,esdTrack->Pt());	  
	  }else{ // material
	    fHistoManager->GetHistoProcess(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,esdTrack->Pt());	  

	  }


	}
      }
    }


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


}

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

}


Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
  // Check if it is a primary and if it was transported
  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;

}

// void AliAnalysisTaskEvil::PrintProcInfo()
// {
//   ProcInfo_t info;
//   gSystem->GetProcInfo(&info);
//   AliInfo(Form("fMemResident=%10ld kB  fMemVirtual=%10ld kB",info.fMemResident,info.fMemVirtual));
// }
Commit	Line	Data
a23f7c97	1	// AliAnalysisTaskMultPbTracks
	2
	3	// Author: Michele Floris, CERN
	4	// TODO:
	5	// - Add chi2/cluster plot for primary, secondaries and fakes
a23f7c97	6
	7	#include "AliAnalysisTaskMultPbTracks.h"
	8	#include "AliESDInputHandler.h"
	9	#include "AliAnalysisMultPbTrackHistoManager.h"
	10	#include "AliAnalysisManager.h"
	11	#include "AliESDtrackCuts.h"
a23f7c97	12	#include "AliMCEvent.h"
	13	#include "AliStack.h"
	14	#include "TH1I.h"
	15	#include "TH3D.h"
bcc49144	16	#include "TH2D.h"
a23f7c97	17	#include "AliMCParticle.h"
a23f7c97	18	#include "AliGenEventHeader.h"
54b31d6a	19	#include "AliESDCentrality.h"
7f5f2e0c	20	#include "AliMultiplicity.h"
a23f7c97	21	#include <iostream>
2bbfd8c6	22	#include "AliAnalysisMultPbCentralitySelector.h"
a23f7c97	23
	24	using namespace std;
	25
	26	ClassImp(AliAnalysisTaskMultPbTracks)
	27
	28	AliAnalysisTaskMultPbTracks::AliAnalysisTaskMultPbTracks()
54b31d6a	29	: AliAnalysisTaskSE("TaskMultPbTracks"),
9441cdaa	30	fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0),fIsTPCOnly(0)
a23f7c97	31	{
	32	// constructor
	33
	34	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	35	DefineOutput(2, AliESDtrackCuts::Class());
2bbfd8c6	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),
9441cdaa	44	fESD(0),fHistoManager(0),fCentrSelector(0),fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0),fIsTPCOnly(0)
a23f7c97	45	{
	46	//
	47	// Standard constructur which should be used
	48	//
	49
	50	DefineOutput(1, AliAnalysisMultPbTrackHistoManager::Class());
54b31d6a	51	DefineOutput(2, AliESDtrackCuts::Class());
2bbfd8c6	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) :
9441cdaa	60	AliAnalysisTaskSE(obj) ,fESD (0), fHistoManager(0), fCentrSelector(0), fTrackCuts(0),fTrackCutsNoDCA(0),fIsMC(0),fIsTPCOnly(0)
a23f7c97	61	{
	62	//copy ctor
	63	fESD = obj.fESD ;
	64	fHistoManager= obj.fHistoManager;
a23f7c97	65	fTrackCuts = obj.fTrackCuts;
a23f7c97	66	fTrackCutsNoDCA = obj.fTrackCutsNoDCA;
2bbfd8c6	67	fCentrSelector = obj.fCentrSelector;
9441cdaa	68	fIsMC = obj.fIsMC;
9441cdaa	69	fIsTPCOnly = obj.fIsTPCOnly;
a23f7c97	70
	71	}
	72
	73	AliAnalysisTaskMultPbTracks::~AliAnalysisTaskMultPbTracks(){
	74	// destructor
	75
	76	if(!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
	77	if(fHistoManager) {
	78	delete fHistoManager;
	79	fHistoManager = 0;
	80	}
	81	}
	82	// Histo list should not be destroyed: fListWrapper is owner!
	83
	84	}
	85	void AliAnalysisTaskMultPbTracks::UserCreateOutputObjects()
	86	{
	87	// Called once
	88
	89	// For the DCA distributions, we want to use exactly the same cuts
	90	// as for the other distributions, with the exception of the DCA cut
	91	fTrackCutsNoDCA = new AliESDtrackCuts(*fTrackCuts); // clone cuts
	92	// Reset all DCA cuts; FIXME: is this all?
	93	fTrackCutsNoDCA->SetMaxDCAToVertexXY();
	94	fTrackCutsNoDCA->SetMaxDCAToVertexZ ();
	95	fTrackCutsNoDCA->SetMaxDCAToVertexXYPtDep();
	96	fTrackCutsNoDCA->SetMaxDCAToVertexZPtDep();
	97
	98	}
	99
	100
	101	void AliAnalysisTaskMultPbTracks::UserExec(Option_t *)
	102	{
	103	// User code
	104
	105	/* PostData(0) is taken care of by AliAnalysisTaskSE */
	106	PostData(1,fHistoManager);
54b31d6a	107	PostData(2,fTrackCuts);
2bbfd8c6	108	PostData(3,fCentrSelector);
a23f7c97	109
54b31d6a	110	// Cache histogram pointers
e0376287	111	static TH3D * hTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
bcc49144	112	static TH2D * hDCA [AliAnalysisMultPbTrackHistoManager::kNHistos];
e0376287	113	static TH1D * hNTracks [AliAnalysisMultPbTrackHistoManager::kNHistos];
a23f7c97	114	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen );
	115	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRec );
	116	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
	117	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
	118	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
	119	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
	120
	121	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoGen );
	122	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRec );
	123	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim );
	124	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecFake );
	125	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat );
	126	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak] = fHistoManager->GetHistoDCA(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak);
	127
2bbfd8c6	128	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoGen );
2bbfd8c6	129	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] = fHistoManager->GetHistoMult(AliAnalysisMultPbTrackHistoManager::kHistoRec );
a23f7c97	130
	131	fESD = dynamic_cast<AliESDEvent*>(fInputEvent);
	132	if (strcmp(fESD->ClassName(),"AliESDEvent")) {
	133	AliFatal("Not processing ESDs");
	134	}
	135
a23f7c97	136	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatAll);
a23f7c97	137
8b628a20	138
54b31d6a	139	// Centrality selection
2bbfd8c6	140	Bool_t isCentralitySelected = fCentrSelector->IsCentralityBinSelected(fESD,fTrackCuts);
	141	if(!isCentralitySelected) return;
	142
	143	// AliESDCentrality *centrality = fESD->GetCentrality();
	144	// if(!centrality) {
	145	// AliFatal("Centrality object not available");
	146	// }
	147	// else {
	148	// Int_t centrBin = centrality->GetCentralityClass5(fCentralityEstimator.Data()) ;
54b31d6a	149
2bbfd8c6	150	// if (centrBin != fCentrBin && fCentrBin != -1) return;
2bbfd8c6	151	// }
a23f7c97	152
	153	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatCentr);
	154
72491d7c	155	Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fOfflineTrigger);
	156
	157	if(!isSelected) return;
	158	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatPhysSel);
4d0aa70f	159
4d0aa70f	160
a23f7c97	161	if (fIsMC) {
4d0aa70f	162	// Get MC vertex
	163	//FIXME: which vertex do I take for MC?
	164	TArrayF vertex;
	165	fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
	166	Float_t zvGen = vertex[2];
a23f7c97	167
	168	if (!fMCEvent) {
	169	AliError("No MC info found");
	170	} else {
	171
e0376287	172	//loop on the MC event, only over primaries, which are always
e0376287	173	// the first in stack.
bcc49144	174	// WRONG: D&B decays are produced in the transport... Need To Loop over all particles
	175	// Int_t nMCTracks = fMCEvent->GetNumberOfPrimaries();
	176	Int_t nMCTracks = fMCEvent->GetNumberOfTracks();
e0376287	177	Int_t nPhysicalPrimaries = 0;
e0376287	178	for (Int_t ipart=0; ipart<nMCTracks; ipart++) {
a23f7c97	179
	180	AliMCParticle mcPart = (AliMCParticle)fMCEvent->GetTrack(ipart);
	181
	182	// We don't care about neutrals and non-physical primaries
	183	if(mcPart->Charge() == 0) continue;
	184
54b31d6a	185	//check if current particle is a physical primary
e0376287	186	if(!IsPhysicalPrimaryAndTransportBit(ipart)) continue;
bcc49144	187	if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) {
	188	// This cures a bug in Hijing
	189	// A little hack here: I put those in the underflow bin of the process type to keep track of them
	190	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(-1);
	191	continue;
	192	}
	193
e0376287	194	nPhysicalPrimaries++;
3b8cbf2d	195	// Fill species histo and particle species
	196	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(mcPart->Particle()->GetUniqueID());
	197	fHistoManager->FillParticleID(AliAnalysisMultPbTrackHistoManager::kHistoGen, mcPart);
e0376287	198
a23f7c97	199	// Float_t zv = vtxESD->GetZ();
a23f7c97	200	// Fill generated histo
4d0aa70f	201	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(mcPart->Pt(),mcPart->Eta(),zvGen);
4d0aa70f	202	Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
bcc49144	203	cout << "Part " << partCode << endl;
4d0aa70f	204	if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus \|\|
bcc49144	205	partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus \|\|
	206	partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus \|\|
	207	partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus \|\|
	208	partCode == AliAnalysisMultPbTrackHistoManager::kPartP \|\|
	209	partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar
	210	){
	211	cout << "Found " << partCode << endl;
	212
4d0aa70f	213	fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoGen, partCode);
bcc49144	214	}
a23f7c97	215	}
e0376287	216	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoGen]->Fill(nPhysicalPrimaries);
4d0aa70f	217	fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoGen)->Fill(zvGen);
	218
	219
a23f7c97	220	}
	221	}
	222
e0376287	223
a23f7c97	224	const AliESDVertex* vtxESD = fESD->GetPrimaryVertex();
a23f7c97	225	if(!vtxESD) return;
bcc49144	226	// TODO: leave the cuts here or find a better place?)
bcc49144	227	// Quality cut on vertexer Z, as suggested by Francesco Prino
4d0aa70f	228	if(vtxESD->IsFromVertexerZ()) {
	229	if (vtxESD->GetNContributors() <= 0) return;
	230	if (vtxESD->GetDispersion() >= 0.04) return;
	231	if (vtxESD->GetZRes() >= 0.25) return;
	232	}
55e8e56a	233	// "Beam gas" vertex cut
55e8e56a	234	const AliESDVertex * vtxESDTPC= fESD->GetPrimaryVertexTPC();
bcc49144	235	if(vtxESDTPC->GetNContributors()<1) return;
55e8e56a	236	if (vtxESDTPC->GetNContributors()<(-10.+0.25*fESD->GetMultiplicity()->GetNumberOfITSClusters(0))) return;
bcc49144	237
bcc49144	238	// Fill vertex and statistics
a23f7c97	239	fHistoManager->GetHistoStats()->Fill(AliAnalysisMultPbTrackHistoManager::kStatVtx);
4d0aa70f	240	fHistoManager->GetHistoVzEvent(AliAnalysisMultPbTrackHistoManager::kHistoRec)->Fill(vtxESD->GetZ());
a23f7c97	241
	242	// loop on tracks
	243	Int_t ntracks = fESD->GetNumberOfTracks();
2bbfd8c6	244	Int_t acceptedTracks = 0;
a23f7c97	245
a23f7c97	246	for (Int_t itrack = 0; itrack<ntracks; itrack++) {
9441cdaa	247	AliESDtrack * esdTrack = fIsTPCOnly ? fTrackCuts->GetTPCOnlyTrack(fESD,itrack) : fESD->GetTrack(itrack);// FIXMEL TrackCuts or TrackCutsNoDCA for the TPC?
	248	if (!esdTrack) continue;
	249
a23f7c97	250	// Fill DCA distibutions, without the DCA cut, Fill the other stuff, with all cuts!
	251	Bool_t accepted = fTrackCuts->AcceptTrack(esdTrack);
	252	Bool_t acceptedNoDCA = fTrackCutsNoDCA->AcceptTrack(esdTrack);
	253
2bbfd8c6	254	if(accepted) acceptedTracks++;
2bbfd8c6	255
a23f7c97	256	// Compute weighted offset
a23f7c97	257	// TODO: other choiches for the DCA?
a23f7c97	258	Double_t weightedDCA = 10;
	259
	260
bcc49144	261	#if defined WEIGHTED_DCA
	262	Double_t b = fESD->GetMagneticField();
	263	Double_t dca[2];
	264	Double_t cov[3];
a23f7c97	265	if (esdTrack->PropagateToDCA(vtxESD, b,3., dca, cov)) {
	266	Double_t det = cov[0]cov[2]-cov[1]cov[1];
	267	if (det<=0) {
	268	AliError("DCA Covariance matrix is not positive definite!");
	269	}
	270	else {
	271	weightedDCA = (dca[0]dca[0]cov[2]+dca[1]dca[1]cov[0]-2dca[0]dca[1]*cov[1])/det;
	272	weightedDCA = weightedDCA>0 ? TMath::Sqrt(weightedDCA) : 10;
	273	}
	274	// printf("dR:%e dZ%e sigR2:%e sigRZ:%e sigZ2:%e\n",dca[0],dca[1],cov[0],cov[1],cov[2]);
	275	}
	276
bcc49144	277	#elif defined TRANSVERSE_DCA
	278	Float_t xz[2];
	279	esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
	280	weightedDCA = xz[0];
	281	#endif
a23f7c97	282
	283
	284	// Alternative: p*DCA
	285	// Float_t xz[2];
	286	// esdTrack->GetDZ(vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(), fESD->GetMagneticField(), xz);
	287	// Float_t dca = xz[0]*esdTrack->P();
	288
bcc49144	289	// This cures a bug in Hijing (displaced primaries)
	290	if (fIsMC) {
	291	Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
	292	if (IsPhysicalPrimaryAndTransportBit(label)) {
	293	AliMCParticle mcPart = (AliMCParticle)fMCEvent->GetTrack(label);
	294	if(!mcPart) AliFatal("No particle");
	295	TArrayF vertex;
	296	fMCEvent->GenEventHeader()->PrimaryVertex(vertex);
	297	Float_t zvGen = vertex[2];
	298	if(TMath::Abs(mcPart->Zv()-zvGen)>1e-6) continue;
	299	}
	300	}
	301	// for each track
	302	if(accepted)
	303	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	304	if(acceptedNoDCA)
	305	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRec]->Fill(weightedDCA,esdTrack->Pt());
a23f7c97	306
	307	// Fill separately primaries and secondaries
	308	// FIXME: fakes? Is this the correct way to account for them?
	309	// Get label and corresponding mcPart;
	310	if (fIsMC) {
72491d7c	311	Int_t label = TMath::Abs(esdTrack->GetLabel()); // no fakes!!!
72491d7c	312	//Int_t label = esdTrack->GetLabel(); //
e0376287	313	AliMCParticle mcPart = label < 0 ? 0 : (AliMCParticle)fMCEvent->GetTrack(label);
a23f7c97	314	if (!mcPart) {
	315	if(accepted)
	316	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	317	if(acceptedNoDCA)
bcc49144	318	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecFake]->Fill(weightedDCA,esdTrack->Pt());
a23f7c97	319	}
a23f7c97	320	else {
e0376287	321	if(IsPhysicalPrimaryAndTransportBit(label)) {
e0376287	322	// Fill species histo
3b8cbf2d	323	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim)->Fill(mcPart->Particle()->GetUniqueID());
4d0aa70f	324	if(accepted) {
a23f7c97	325	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
4d0aa70f	326	Int_t partCode = fHistoManager->GetLocalParticleID(mcPart);
4d0aa70f	327	if (partCode == AliAnalysisMultPbTrackHistoManager::kPartPiPlus \|\|
bcc49144	328	partCode == AliAnalysisMultPbTrackHistoManager::kPartPiMinus \|\|
	329	partCode == AliAnalysisMultPbTrackHistoManager::kPartKPlus \|\|
	330	partCode == AliAnalysisMultPbTrackHistoManager::kPartKMinus \|\|
	331	partCode == AliAnalysisMultPbTrackHistoManager::kPartP \|\|
	332	partCode == AliAnalysisMultPbTrackHistoManager::kPartPBar
	333	)
4d0aa70f	334	fHistoManager->GetHistoPtEtaVz(AliAnalysisMultPbTrackHistoManager::kHistoRecPrim, partCode);
4d0aa70f	335	}
a23f7c97	336	if(acceptedNoDCA)
bcc49144	337	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecPrim]->Fill(weightedDCA,esdTrack->Pt());
a23f7c97	338	}
	339	else {
	340	Int_t mfl=0;
	341	Int_t indexMoth=mcPart->Particle()->GetFirstMother();
	342	if(indexMoth>=0){
	343	TParticle* moth = fMCEvent->Stack()->Particle(indexMoth);
	344	Float_t codemoth = TMath::Abs(moth->GetPdgCode());
	345	mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
	346	}
	347	if(mfl==3){ // strangeness
3b8cbf2d	348	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	349	if(accepted)
	350	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	351	if(acceptedNoDCA)
bcc49144	352	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecWeak]->Fill(weightedDCA,esdTrack->Pt());
a23f7c97	353	}else{ // material
3b8cbf2d	354	fHistoManager->GetHistoProcess(AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat)->Fill(mcPart->Particle()->GetUniqueID());
a23f7c97	355	if(accepted)
	356	hTracks[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(esdTrack->Pt(),esdTrack->Eta(),vtxESD->GetZ());
	357	if(acceptedNoDCA)
bcc49144	358	hDCA[AliAnalysisMultPbTrackHistoManager::kHistoRecSecMat]->Fill(weightedDCA,esdTrack->Pt());
a23f7c97	359
	360	}
	361
	362
	363	}
	364	}
	365	}
	366
	367
	368	}
2bbfd8c6	369	// cout << acceptedTracks << endl;
	370
	371	hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] ->Fill(acceptedTracks);
7f5f2e0c	372	// FIXME
7f5f2e0c	373	// hNTracks[AliAnalysisMultPbTrackHistoManager::kHistoRec] ->Fill(fESD->GetMultiplicity()->GetNumberOfTracklets());
a23f7c97	374
	375
	376	}
	377
	378	void AliAnalysisTaskMultPbTracks::Terminate(Option_t *){
	379	// terminate
	380
	381	}
	382
	383
e0376287	384	Bool_t AliAnalysisTaskMultPbTracks::IsPhysicalPrimaryAndTransportBit(Int_t ipart) {
bcc49144	385	// Check if it is a primary and if it was transported
e0376287	386	Bool_t physprim=fMCEvent->IsPhysicalPrimary(ipart);
	387	if (!physprim) return kFALSE;
	388	Bool_t transported = ((AliMCParticle*)fMCEvent->GetTrack(ipart))->Particle()->TestBit(kTransportBit);
	389	if(!transported) return kFALSE;
	390
	391	return kTRUE;
	392
	393	}
eef42d18	394
	395	// void AliAnalysisTaskEvil::PrintProcInfo()
	396	// {
	397	// ProcInfo_t info;
	398	// gSystem->GetProcInfo(&info);
	399	// AliInfo(Form("fMemResident=%10ld kB fMemVirtual=%10ld kB",info.fMemResident,info.fMemVirtual));
	400	// }