+++ /dev/null
-#include "AliAnalysisTaskHighPtDeDx.h"
-
-// ROOT includes
-#include <TList.h>
-#include <TTree.h>
-#include <TMath.h>
-#include <TH1.h>
-#include <TParticle.h>
-#include <TFile.h>
-
-// AliRoot includes
-#include <AliAnalysisManager.h>
-#include <AliAnalysisFilter.h>
-#include <AliESDInputHandler.h>
-#include <AliESDEvent.h>
-#include <AliESDVertex.h>
-#include <AliLog.h>
-#include <AliExternalTrackParam.h>
-#include <AliESDtrackCuts.h>
-#include <AliESDVZERO.h>
-#include <AliAODVZERO.h>
-
-#include <AliMCEventHandler.h>
-#include <AliMCEvent.h>
-#include <AliStack.h>
-
-#include <TTreeStream.h>
-
-#include <AliHeader.h>
-#include <AliGenPythiaEventHeader.h>
-#include <AliGenDPMjetEventHeader.h>
-
-#include "AliCentrality.h"
-
-#include <AliAODTrack.h>
-#include <AliAODPid.h>
-#include <AliAODMCHeader.h>
-
-
-// STL includes
-#include <iostream>
-using namespace std;
-
-
-//
-// Responsible:
-// Alexandru Dobrin (Wayne State)
-// Peter Christiansen (Lund)
-//
-
-/*
-To do:
-
-Separate the code into two
-
-*/
-
-
-
-
-ClassImp(AliAnalysisTaskHighPtDeDx)
-
-const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2;
-
-//_____________________________________________________________________________
-AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx():
- AliAnalysisTaskSE(),
- fESD(0x0),
- fAOD(0x0),
- fMC(0x0),
- fMCStack(0x0),
- fMCArray(0x0),
- fTrackFilter(0x0),
- fTrackFilterGolden(0x0),
- fTrackFilterTPC(0x0),
- fAnalysisType("ESD"),
- fAnalysisMC(kFALSE),
- fAnalysisPbPb(kFALSE),
- fTPCBranch(kFALSE),
- ftrigBit1(0x0),
- ftrigBit2(0x0),
- fRandom(0x0),
- fEvent(0x0),
- fTrackArrayGlobalPar(0x0),
- fTrackArrayTPCPar(0x0),
- fTrackArrayMC(0x0),
- fVZEROArray(0x0),
- fVtxCut(10.0),
- fEtaCut(0.9),
- fMinPt(0.1),
- fMinCent(0.0),
- fMaxCent(100.0),
- fLowPtFraction(0.01),
- fTreeOption(0),
- fMcProcessType(-999),
- fTriggeredEventMB(-999),
- fVtxStatus(-999),
- fZvtx(-999),
- fZvtxMC(-999),
- fRun(-999),
- fEventId(-999),
- fListOfObjects(0),
- fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
- fTree(0x0),
- fn1(0),
- fn2(0),
- fcent(0)
-
-
-{
- // Default constructor (should not be used)
-
- // fRandom = new TRandom(0); // 0 means random seed
-}
-
-//______________________________________________________________________________
-AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name):
- AliAnalysisTaskSE(name),
- fESD(0x0),
- fAOD(0x0),
- fMC(0x0),
- fMCStack(0x0),
- fMCArray(0x0),
- fTrackFilter(0x0),
- fTrackFilterGolden(0x0),
- fTrackFilterTPC(0x0),
- fAnalysisType("ESD"),
- fAnalysisMC(kFALSE),
- fAnalysisPbPb(kFALSE),
- fTPCBranch(kFALSE),
- ftrigBit1(0x0),
- ftrigBit2(0x0),
- fRandom(0x0),
- fEvent(0x0),
- fTrackArrayGlobalPar(0x0),
- fTrackArrayMC(0x0),
- fVtxCut(10.0),
- fEtaCut(0.9),
- fMinPt(0.1),
- fLowPtFraction(0.01),
- fTreeOption(0),
- fMcProcessType(-999),
- fTriggeredEventMB(-999),
- fVtxStatus(-999),
- fZvtx(-999),
- fZvtxMC(-999),
- fRun(-999),
- fEventId(-999),
- fListOfObjects(0),
- fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
- fTree(0x0),
- fn1(0),
- fn2(0),
- fcent(0)
-
-{
-
- if(fTree)fTree=0;
- // Output slot #1 writes into a TList
- DefineOutput(1, TList::Class());
-}
-
-//_____________________________________________________________________________
-AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx()
-{
- // Destructor
- // histograms are in the output list and deleted when the output
- // list is deleted by the TSelector dtor
- if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
- delete fListOfObjects;
- fListOfObjects = 0;
- }
- if (fRandom) delete fRandom;
- fRandom=0;
-
- // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse
- // if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
-
-
-
-}
-
-//______________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects()
-{
- // This method is called once per worker node
- // Here we define the output: histograms and debug tree if requested
- // We also create the random generator here so it might get different seeds...
- fRandom = new TRandom(0); // 0 means random seed
-
- OpenFile(1);
- fListOfObjects = new TList();
- fListOfObjects->SetOwner();
-
- //
- // Histograms
- //
- fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
- fListOfObjects->Add(fEvents);
-
- fn1=new TH1F("fn1","fn1",5001,-1,5000);
- fListOfObjects->Add(fn1);
-
- fn2=new TH1F("fn2","fn2",5001,-1,5000);
- fListOfObjects->Add(fn2);
-
- fcent=new TH1F("fcent","fcent",104,-2,102);
- fListOfObjects->Add(fcent);
-
- fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
- fListOfObjects->Add(fVtx);
-
- fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
- fListOfObjects->Add(fVtxBeforeCuts);
-
- fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
- fListOfObjects->Add(fVtxAfterCuts);
-
- if (fAnalysisMC) {
- fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
- fListOfObjects->Add(fVtxMC);
- }
-
- if (fTreeOption) {
-
- fTree = new TTree("tree","Event data");
- fEvent = new DeDxEvent();
- fTree->Branch("event", &fEvent);
-
- fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000);
- fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar);
- if(fTPCBranch){
- fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000);
- fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar);
- }
-
- fVZEROArray = new TClonesArray("VZEROCell", 1000);
- fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray);
-
- if (fAnalysisMC) {
- fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000);
- fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC);
- }
-
- fTree->SetDirectory(0);
-
- fListOfObjects->Add(fTree);
-
- }
-
- // Post output data.
- PostData(1, fListOfObjects);
-}
-
-//______________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *)
-{
- // Main loop
-
- //
- // First we make sure that we have valid input(s)!
- //
- AliVEvent *event = InputEvent();
- if (!event) {
- Error("UserExec", "Could not retrieve event");
- return;
- }
-
-
- if (fAnalysisType == "ESD"){
- fESD = dynamic_cast<AliESDEvent*>(event);
- if(!fESD){
- Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
- this->Dump();
- return;
- }
- } else {
- fAOD = dynamic_cast<AliAODEvent*>(event);
- if(!fAOD){
- Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
- this->Dump();
- return;
- }
- }
-
-
-
- if (fAnalysisMC) {
-
- if (fAnalysisType == "ESD"){
- fMC = dynamic_cast<AliMCEvent*>(MCEvent());
- if(!fMC){
- Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
- this->Dump();
- return;
- }
-
- fMCStack = fMC->Stack();
-
- if(!fMCStack){
- Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
- this->Dump();
- return;
- }
- } else { // AOD
-
- fMC = dynamic_cast<AliMCEvent*>(MCEvent());
- if(fMC)
- fMC->Dump();
-
- fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
- if(!fMCArray){
- Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
- this->Dump();
- return;
- }
- }
- }
-
-
- // Get trigger decision
- fTriggeredEventMB = 0; //init
-
-
- if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
- ->IsEventSelected() & ftrigBit1 ){
- fn1->Fill(1);
- fTriggeredEventMB = 1; //event triggered as minimum bias
- }
- if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
- ->IsEventSelected() & ftrigBit2 ){
- // From AliVEvent:
- // kINT7 = BIT(1), // V0AND trigger, offline V0 selection
- fTriggeredEventMB += 2;
- fn2->Fill(1);
- }
-
-
- // Get process type for MC
- fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
-
- // real data that are not triggered we skip
- if(!fAnalysisMC && !fTriggeredEventMB)
- return;
-
- if (fAnalysisMC) {
-
- if (fAnalysisType == "ESD"){
-
- AliHeader* headerMC = fMC->Header();
- if (headerMC) {
-
- AliGenEventHeader* genHeader = headerMC->GenEventHeader();
- TArrayF vtxMC(3); // primary vertex MC
- vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
- if (genHeader) {
- genHeader->PrimaryVertex(vtxMC);
- }
- fZvtxMC = vtxMC[2];
-
- // PYTHIA:
- AliGenPythiaEventHeader* pythiaGenHeader =
- dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
- if (pythiaGenHeader) { //works only for pythia
- fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
- }
- // PHOJET:
- AliGenDPMjetEventHeader* dpmJetGenHeader =
- dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
- if (dpmJetGenHeader) {
- fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
- }
- }
- } else { // AOD
-
- AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
- if(mcHeader) {
- fZvtxMC = mcHeader->GetVtxZ();
-
- if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
- fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
- } else {
- fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
- }
- }
- }
- }
-
- // There are 3 cases
- // Vertex: NO - status -1
- // Vertex: YES : outside cut - status 0
- // : inside cut - status 1
- // We have to be careful how we normalize because we probably want to
- // normalize to:
- // Nevents=(No vertex + outside + inside)/(out + in)*in
-
-
- if (fAnalysisType == "ESD")
- fZvtx = GetVertex(fESD);
- else // AOD
- fZvtx = GetVertex(fAOD);
-
- fVtxStatus = -999;
-
- if(fZvtx<-990) {
-
- fVtxStatus = -1;
- if(fTriggeredEventMB)
- fVtx->Fill(0);
- if(fAnalysisMC)
- fVtxMC->Fill(0);
- } else {
-
- if(fTriggeredEventMB)
- fVtx->Fill(1);
- if(fAnalysisMC)
- fVtxMC->Fill(1);
- fVtxBeforeCuts->Fill(fZvtx);
- fVtxStatus = 0;
- if (TMath::Abs(fZvtx) < fVtxCut) {
- fVtxAfterCuts->Fill(fZvtx);
- fVtxStatus = 1;
- }
- }
-
- // store MC event data no matter what
- if(fAnalysisMC) {
-
- if (fAnalysisType == "ESD") {
- ProcessMCTruthESD();
- } else { // AOD
- ProcessMCTruthAOD();
- }
- }
-
-
-
- Float_t centrality = -10;
- // only analyze triggered events
- if(fTriggeredEventMB) {
-
- if (fAnalysisType == "ESD"){
- if(fAnalysisPbPb){
- AliCentrality *centObject = fESD->GetCentrality();
- centrality = centObject->GetCentralityPercentile("V0M");
- if((centrality>fMaxCent)||(centrality<fMinCent))return;
- }
- fcent->Fill(centrality);
- AnalyzeESD(fESD);
- } else { // AOD
- if(fAnalysisPbPb){
- AliCentrality *centObject = fAOD->GetCentrality();
- if(centObject)
- centrality = centObject->GetCentralityPercentile("V0M");
- if((centrality>fMaxCent)||(centrality<fMinCent))return;
- }
- fcent->Fill(centrality);
- AnalyzeAOD(fAOD);
- }
- }
-
- if( fTreeOption) {
-
- fEvent->process = fMcProcessType;
- fEvent->trig = fTriggeredEventMB;
- fEvent->zvtxMC = fZvtxMC;
- fEvent->cent = centrality;
-
- fTree->Fill();
- fTrackArrayGlobalPar->Clear();
- if(fTPCBranch)
- fTrackArrayTPCPar->Clear();
- fVZEROArray->Clear();
-
- if (fAnalysisMC)
- fTrackArrayMC->Clear();
- }
-
- // Post output data.
- PostData(1, fListOfObjects);
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
-{
- fRun = esdEvent->GetRunNumber();
- fEventId = 0;
- if(esdEvent->GetHeader())
- fEventId = GetEventIdAsLong(esdEvent->GetHeader());
-
- Short_t isPileup = esdEvent->IsPileupFromSPD();
-
- // Int_t event = esdEvent->GetEventNumberInFile();
- UInt_t time = esdEvent->GetTimeStamp();
- // ULong64_t trigger = esdEvent->GetTriggerMask();
- Float_t magf = esdEvent->GetMagneticField();
-
-
-
-
-
- if(fTriggeredEventMB) {// Only MC case can we have not triggered events
-
- // accepted event
- fEvents->Fill(0);
-
-
- if(fVtxStatus!=1) return; // accepted vertex
- Int_t nESDTracks = esdEvent->GetNumberOfTracks();
-
- if(nESDTracks<1)return;
- cout<<"Multiplicity="<<nESDTracks<<endl;
-
- ProduceArrayTrksESD( esdEvent, kGlobalTrk );//produce array with global track parameters
- if(fTPCBranch)
- ProduceArrayTrksESD( esdEvent, kTPCTrk );//array with tpc track parametes
-
- fEvents->Fill(1);
- AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
-
- for (Int_t iCh=0; iCh<64; ++iCh) {
- Float_t multv=esdV0->GetMultiplicity(iCh);
- Int_t intexv=iCh;
- VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
- cellv0->cellindex=intexv;
- cellv0->cellmult= multv;
- }
-
-
-
- } // end if triggered
-
- if(fTreeOption) {
-
- fEvent->run = fRun;
- fEvent->eventid = fEventId;
- fEvent->time = time;
- //fEvent->cent = centrality;
- fEvent->mag = magf;
- fEvent->zvtx = fZvtx;
- fEvent->vtxstatus = fVtxStatus;
- fEvent->pileup = isPileup;
-
- }
-
-
-
-
-}
-
-//________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
-{
- fRun = aodEvent->GetRunNumber();
- fEventId = 0;
- if(aodEvent->GetHeader())
- fEventId = GetEventIdAsLong(aodEvent->GetHeader());
-
- UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
- Float_t magf = aodEvent->GetMagneticField();
-
- //Int_t trackmult = 0; // no pt cuts
- //Int_t nadded = 0;
-
- Short_t isPileup = aodEvent->IsPileupFromSPD();
-
-
-
-
- if(fTriggeredEventMB) {// Only MC case can we have not triggered events
-
- // accepted event
- fEvents->Fill(0);
-
- if(fVtxStatus!=1) return; // accepted vertex
- Int_t nAODTracks = aodEvent->GetNumberOfTracks();
- if(nAODTracks<1)return;
-
- ProduceArrayTrksAOD( aodEvent, kGlobalTrk );
- if(fTPCBranch)
- ProduceArrayTrksAOD( aodEvent, kTPCTrk );
- fEvents->Fill(1);
-
- AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad
-
- for (Int_t iCh=0; iCh<64; ++iCh) {
- Float_t multv=esdV0->GetMultiplicity(iCh);
- Int_t intexv=iCh;
- VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell();
- cellv0->cellindex=intexv;
- cellv0->cellmult= multv;
- }
-
-
-
-
- } // end if triggered
-
- if(fTreeOption) {
-
- //Sort(fTrackArrayGlobalPar, kFALSE);
-
- fEvent->run = fRun;
- fEvent->eventid = fEventId;
- fEvent->time = time;
- //fEvent->cent = centrality;
- fEvent->mag = magf;
- fEvent->zvtx = fZvtx;
- fEvent->vtxstatus = fVtxStatus;
- //fEvent->trackmult = trackmult;
- //fEvent->n = nadded;
- fEvent->pileup = isPileup;
- }
-}
-
-//_____________________________________________________________________________
-Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const
-{
- Float_t zvtx = -999;
-
- const AliVVertex* primaryVertex = event->GetPrimaryVertex();
-
- if(primaryVertex->GetNContributors()>0)
- zvtx = primaryVertex->GetZ();
-
- return zvtx;
-}
-
-//_____________________________________________________________________________
-Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const
-{
- // return our internal code for pions, kaons, and protons
-
- Short_t pidCode = 6;
-
- switch (TMath::Abs(pdgCode)) {
- case 211:
- pidCode = 1; // pion
- break;
- case 321:
- pidCode = 2; // kaon
- break;
- case 2212:
- pidCode = 3; // proton
- break;
- case 11:
- pidCode = 4; // electron
- break;
- case 13:
- pidCode = 5; // muon
- break;
- default:
- pidCode = 6; // something else?
- };
-
- return pidCode;
-}
-
-
-//_____________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD()
-{
- // Fill the special MC histogram with the MC truth info
-
- Short_t trackmult = 0;
- Short_t nadded = 0;
- const Int_t nTracksMC = fMCStack->GetNtrack();
-
- for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
-
- //Cuts
- if(!(fMCStack->IsPhysicalPrimary(iTracks)))
- continue;
-
- TParticle* trackMC = fMCStack->Particle(iTracks);
-
- Double_t chargeMC = trackMC->GetPDG()->Charge();
- if (chargeMC == 0)
- continue;
-
- if (TMath::Abs(trackMC->Eta()) > fEtaCut )
- continue;
-
- trackmult++;
-
- // "charge:pt:p:eta:phi:pidCode"
- Float_t ptMC = trackMC->Pt();
- Float_t pMC = trackMC->P();
- Float_t etaMC = trackMC->Eta();
- Float_t phiMC = trackMC->Phi();
-
- Int_t pdgCode = trackMC->GetPdgCode();
- Short_t pidCodeMC = 0;
- pidCodeMC = GetPidCode(pdgCode);
-
- // Here we want to add some of the MC histograms!
-
- // And therefore we first cut here!
- if (trackMC->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the MC histograms!
- // }
-
- if(fTreeOption) {
-
- DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
- nadded++;
-
- track->pMC = pMC;
- track->ptMC = ptMC;
- track->etaMC = etaMC;
- track->phiMC = phiMC;
- track->qMC = Short_t(chargeMC);
- track->pidMC = pidCodeMC;
- track->pdgMC = pdgCode;
- }
-
- }//MC track loop
-
- if(fTreeOption) {
-
- Sort(fTrackArrayMC, kTRUE);
-
- fEvent->trackmultMC = trackmult;
- fEvent->nMC = nadded;
- }
-
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD()
-{
- // Fill the special MC histogram with the MC truth info
-
- Short_t trackmult = 0;
- Short_t nadded = 0;
- const Int_t nTracksMC = fMCArray->GetEntriesFast();
-
- for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
-
- AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
-
- //Cuts
- if(!(trackMC->IsPhysicalPrimary()))
- continue;
-
-
- Double_t chargeMC = trackMC->Charge();
- if (chargeMC == 0)
- continue;
-
- if (TMath::Abs(trackMC->Eta()) > fEtaCut )
- continue;
-
- trackmult++;
-
- // "charge:pt:p:eta:phi:pidCode"
- Float_t ptMC = trackMC->Pt();
- Float_t pMC = trackMC->P();
- Float_t etaMC = trackMC->Eta();
- Float_t phiMC = trackMC->Phi();
-
- Int_t pdgCode = trackMC->PdgCode();
- Short_t pidCodeMC = 0;
- pidCodeMC = GetPidCode(pdgCode);
-
- // Here we want to add some of the MC histograms!
-
- // And therefore we first cut here!
- if (trackMC->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the MC histograms!
- // }
-
- if(fTreeOption) {
-
- DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC();
- nadded++;
-
- track->pMC = pMC;
- track->ptMC = ptMC;
- track->etaMC = etaMC;
- track->phiMC = phiMC;
- track->qMC = Short_t(chargeMC);
- track->pidMC = pidCodeMC;
- track->pdgMC = pdgCode;
- }
-
- }//MC track loop
-
- if(fTreeOption) {
-
- Sort(fTrackArrayMC, kTRUE);
-
- fEvent->trackmultMC = trackmult;
- fEvent->nMC = nadded;
- }
-
-}
-
-//_____________________________________________________________________________
-Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
- //
- // Get the process type of the event. PYTHIA
- //
- // source PWG0 dNdpt
-
- Short_t globalType = -1; //init
-
- if(pythiaType==92||pythiaType==93){
- globalType = 2; //single diffractive
- }
- else if(pythiaType==94){
- globalType = 3; //double diffractive
- }
- //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
- else {
- globalType = 1; //non diffractive
- }
- return globalType;
-}
-
-//_____________________________________________________________________________
-Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
- //
- // get the process type of the event. PHOJET
- //
- //source PWG0 dNdpt
- // can only read pythia headers, either directly or from cocktalil header
- Short_t globalType = -1;
-
- if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
- globalType = 1;
- }
- else if (dpmJetType==5 || dpmJetType==6) {
- globalType = 2;
- }
- else if (dpmJetType==7) {
- globalType = 3;
- }
- return globalType;
-}
-
-//_____________________________________________________________________________
-ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
-{
- // To have a unique id for each event in a run!
- // Modified from AliRawReader.h
- return ((ULong64_t)header->GetBunchCrossNumber()+
- (ULong64_t)header->GetOrbitNumber()*3564+
- (ULong64_t)header->GetPeriodNumber()*16777215*3564);
-}
-
-
-//____________________________________________________________________
-TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
-{
- //
- // Finds the first mother among the primary particles of the particle identified by <label>,
- // i.e. the primary that "caused" this particle
- //
- // Taken from AliPWG0Helper class
- //
-
- Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
- if (motherLabel < 0)
- return 0;
-
- return stack->Particle(motherLabel);
-}
-
-//____________________________________________________________________
-Int_t AliAnalysisTaskHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
-{
- //
- // Finds the first mother among the primary particles of the particle identified by <label>,
- // i.e. the primary that "caused" this particle
- //
- // returns its label
- //
- // Taken from AliPWG0Helper class
- //
- const Int_t nPrim = stack->GetNprimary();
-
- while (label >= nPrim) {
-
- //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
-
- TParticle* particle = stack->Particle(label);
- if (!particle) {
-
- AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
- return -1;
- }
-
- // find mother
- if (particle->GetMother(0) < 0) {
-
- AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
- return -1;
- }
-
- label = particle->GetMother(0);
- }
-
- return label;
-}
-
-//____________________________________________________________________
-AliAODMCParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
-{
- //
- // Finds the first mother among the primary particles of the particle identified by <label>,
- // i.e. the primary that "caused" this particle
- //
- // Taken from AliPWG0Helper class
- //
-
- AliAODMCParticle* mcPart = startParticle;
-
- while (mcPart)
- {
-
- if(mcPart->IsPrimary())
- return mcPart;
-
- Int_t mother = mcPart->GetMother();
-
- mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
- }
-
- return 0;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::Sort(TClonesArray* array, Bool_t isMC)
-{
- const Int_t n = array->GetEntriesFast();
- if(n==0) {
-
- if(isMC)
- fEvent->ptmaxMC = 0;
- else
- fEvent->ptmax = 0;
-
- return;
- }
-
- Float_t ptArray[n];
- Int_t indexArray[n];
-
- for(Int_t i = 0; i < n; i++) {
-
- Float_t pt = 0;
- if(isMC) {
- DeDxTrackMC* track = (DeDxTrackMC*)array->At(i);
- pt = track->ptMC;
- } else {
- DeDxTrack* track = (DeDxTrack*)array->At(i);
- pt = track->pt;
- }
- ptArray[i] = pt;
- indexArray[i] = i;
- }
-
- TMath::Sort(n, ptArray, indexArray);
-
- // set max pt
- if(isMC) {
- DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[0]);
- fEvent->ptmaxMC = track->ptMC;
- } else {
- DeDxTrack* track = (DeDxTrack*)array->At(indexArray[0]);
- fEvent->ptmax = track->pt;
- }
-
- // set order of each track
- for(Int_t i = 0; i < n; i++) {
-
- if(isMC) {
- DeDxTrackMC* track = (DeDxTrackMC*)array->At(indexArray[i]);
- track->orderMC = i;
- } else {
- DeDxTrack* track = (DeDxTrack*)array->At(indexArray[i]);
- track->order = i;
- }
- }
-}
-//__________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent, AnalysisMode analysisMode ){
-
- const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
- Int_t trackmult=0;
- Int_t nadded=0;
- if( analysisMode == kGlobalTrk ){
- if(fTrackArrayGlobalPar)
- fTrackArrayGlobalPar->Clear();
- } else if( analysisMode == kTPCTrk ){
- if(fTrackArrayTPCPar)
- fTrackArrayTPCPar->Clear();
- }
-
- if( analysisMode==kGlobalTrk ){
-
- for(Int_t iT = 0; iT < nESDTracks; iT++) {
-
- AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
-
-
- if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
- continue;
-
- UShort_t filterFlag = 0;
- Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
- Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
- Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
-
- UInt_t selectDebug = 0;
- if (fTrackFilterGolden) {
- selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
- if (selectDebug) {
- filterFlag +=1;
- filterCut_Set3=kTRUE;
- }
- }
-
- if (fTrackFilterTPC) {
-
- selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
- if (selectDebug){//only tracks which pass the TPC-only track cuts
- trackmult++;
- filterFlag +=2;
- filterCut_Set1=kTRUE;
-
- }
-
- }
-
- if (fTrackFilter) {
- selectDebug = fTrackFilter->IsSelected(esdTrack);
- if (selectDebug) {
- filterCut_Set2=kTRUE;
- }
- }
-
-
- if(filterFlag==0)
- continue;
-
-
- //
- // Track was accepted
- //
-
- // Here we want to add histograms!
-
- if (esdTrack->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the histograms!
- // }
-
- Short_t charge = esdTrack->Charge();
- Float_t pt = esdTrack->Pt();
- Float_t p = esdTrack->P();
- Float_t eta = esdTrack->Eta();
- Float_t phi = esdTrack->Phi();
- Short_t ncl = esdTrack->GetTPCsignalN();
- Short_t neff = Short_t(esdTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
- // Short_t nclf = esdTrack->GetTPCNclsF();
- Float_t dedx = esdTrack->GetTPCsignal();
- Float_t tpcchi = 0;
- if(esdTrack->GetTPCNcls() > 0)
- tpcchi = esdTrack->GetTPCchi2()/Float_t(esdTrack->GetTPCNcls());
- Float_t b[2];
- Float_t bCov[3];
- esdTrack->GetImpactParameters(b,bCov);
- Float_t dcaxy = b[0];
- Float_t dcaz = b[1];
- Double_t p_con[3] = {0, 0, 0};
- esdTrack->GetConstrainedPxPyPz(p_con);
-
-
- // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
- // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
- // Float_t pttpc = tpcParam->Pt();
- // Float_t ptpc = tpcParam->P();
-
- Float_t ptMC = 0;
- Short_t pidCode = 0; // 0 = real data / no mc track!
- Short_t primaryFlag = 0; // 0 = real data / not primary mc track
- Int_t pdgMother = 0;
-
-
- //with Globaltrack parameters????
- if(fAnalysisMC) {
-
- const Int_t label = TMath::Abs(esdTrack->GetLabel());
- TParticle* mcTrack = fMCStack->Particle(label);
- if (mcTrack){
-
- if(fMCStack->IsPhysicalPrimary(label))
- primaryFlag = 1;
-
- Int_t pdgCode = mcTrack->GetPdgCode();
- pidCode = GetPidCode(pdgCode);
-
- ptMC = mcTrack->Pt();
-
- TParticle* mother = FindPrimaryMother(fMCStack, label);
- pdgMother = mother->GetPdgCode();
- }
- }
-
-
- //TOF
- Float_t beta = -99;
- if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
- if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
- beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
- }
-
- if(fTreeOption) {
-
- DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
- nadded++;
-
- track->p = p;
- track->pt = pt;
- // track->ptcon = pt_con;
- // track->tpcchi = tpcchi;
- track->eta = eta;
- track->phi = phi;
- track->q = charge;
- track->filter = filterFlag;
- track->ncl = ncl;
- track->neff = neff;
- track->dedx = dedx;
- track->beta = beta;
- track->dcaxy = dcaxy;
- track->dcaz = dcaz;
- track->pid = pidCode;
- track->primary = primaryFlag;
- track->pttrue = ptMC;
- track->mother = pdgMother;
- track->filterset1 = filterCut_Set1;
- track->filterset2 = filterCut_Set2;
- track->filterset3 = filterCut_Set3;
-
-
- }
- }//end of track loop
-
- }//end global
-
- else if( analysisMode==kTPCTrk ){
- const AliESDVertex *vtxSPD = ESDevent->GetPrimaryVertexSPD();
- if( vtxSPD->GetNContributors() < 1 || TMath::Abs(vtxSPD->GetZ()) > 10.0 ) return;
-
-
- for(Int_t iT = 0; iT < nESDTracks; iT++) {
-
- AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
-
- AliESDtrack *trackTPC = AliESDtrackCuts::GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(ESDevent),esdTrack->GetID());
-
- if(!trackTPC) continue;
-
- UInt_t selectDebug = 0;
- selectDebug = fTrackFilterTPC->IsSelected(trackTPC);
- if(selectDebug==0) continue;
-
- if(trackTPC->Pt()>0.){
- // only constrain tracks above threshold
- AliExternalTrackParam exParam;
- // take the B-field from the ESD, no 3D fieldMap available at this point
- Bool_t relate = false;
- relate = trackTPC->RelateToVertexTPC(vtxSPD,ESDevent->GetMagneticField(),
- kVeryBig,&exParam);
- if(!relate){
- delete trackTPC;
- continue;
- }
- trackTPC->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),
- exParam.GetCovariance());
- }
- else continue;
-
-
-
-
- if(TMath::Abs(trackTPC->Eta()) > fEtaCut)
- continue;
-
- //
- // Track was accepted
- //
-
- // Here we want to add histograms!
-
- if (trackTPC->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the histograms!
- // }
-
- Short_t charge = trackTPC->Charge();
- Float_t pt = trackTPC->Pt();
- Float_t p = trackTPC->P();
- Float_t eta = trackTPC->Eta();
- Float_t phi = trackTPC->Phi();
- Short_t ncl = trackTPC->GetTPCsignalN();
- Short_t neff = Short_t(trackTPC->GetTPCClusterInfo(2, 1)); // effective track length for pT res
- // Short_t nclf = esdTrack->GetTPCNclsF();
- Float_t dedx = trackTPC->GetTPCsignal();
- Float_t tpcchi = 0;
- if(trackTPC->GetTPCNcls() > 0)
- tpcchi = trackTPC->GetTPCchi2()/Float_t(trackTPC->GetTPCNcls());
- //Float_t b[2];
- //Float_t bCov[3];
- //trackTPC->GetImpactParameters(b,bCov);
- //Float_t dcaxy = b[0];
- //Float_t dcaz = b[1];
-
- Float_t dcaxy = 0.;
- Float_t dcaz = 0.;
- trackTPC->GetImpactParameters(dcaxy,dcaz);
-
-
- Double_t p_con[3] = {0, 0, 0};
- trackTPC->GetConstrainedPxPyPz(p_con);
-
-
- // Float_t pt_con = (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
- // const AliExternalTrackParam* tpcParam = esdTrack->GetTPCInnerParam();
- // Float_t pttpc = tpcParam->Pt();
- // Float_t ptpc = tpcParam->P();
-
- Float_t ptMC = 0;
- Short_t pidCode = 0; // 0 = real data / no mc track!
- Short_t primaryFlag = 0; // 0 = real data / not primary mc track
- Int_t pdgMother = 0;
-
-
- //with Globaltrack parameters????
- if(fAnalysisMC) {
-
- const Int_t label = TMath::Abs(esdTrack->GetLabel());
- TParticle* mcTrack = fMCStack->Particle(label);
- if (mcTrack){
-
- if(fMCStack->IsPhysicalPrimary(label))
- primaryFlag = 1;
-
- Int_t pdgCode = mcTrack->GetPdgCode();
- pidCode = GetPidCode(pdgCode);
-
- ptMC = mcTrack->Pt();
-
- TParticle* mother = FindPrimaryMother(fMCStack, label);
- pdgMother = mother->GetPdgCode();
- }
- }
-
-
- //TOF
- Float_t beta = -99;
- if (esdTrack->GetStatus()&AliESDtrack::kTOFpid){
- if ((esdTrack->GetIntegratedLength() != 0) && (esdTrack->GetTOFsignal() != 0))
- beta = esdTrack->GetIntegratedLength()/esdTrack->GetTOFsignal()/fgkClight;
- }
-
- if(fTreeOption) {
-
- DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
- nadded++;
-
- tracktpc->p = p;
- tracktpc->pt = pt;
- // track->ptcon = pt_con;
- // track->tpcchi = tpcchi;
- tracktpc->eta = eta;
- tracktpc->phi = phi;
- tracktpc->q = charge;
- tracktpc->filter = 1;
- tracktpc->ncl = ncl;
- tracktpc->neff = neff;
- tracktpc->dedx = dedx;
- tracktpc->beta = beta;//computed with Global tracks
- tracktpc->dcaxy = dcaxy;
- tracktpc->dcaz = dcaz;
- tracktpc->pid = pidCode;
- tracktpc->primary = primaryFlag;
- tracktpc->pttrue = ptMC;
- tracktpc->mother = pdgMother;
- tracktpc->filterset1 = 0;
- tracktpc->filterset2 = 0;
- tracktpc->filterset3 = 0;
-
- }
-
-
- if(trackTPC)
- delete trackTPC;
-
-
- }//end of track loop
-
- }//end of: if isglobal==kFALSE
-
-
- if(fTreeOption) {
-
- if( analysisMode==kGlobalTrk ){
- Sort(fTrackArrayGlobalPar, kFALSE);
-
- fEvent->trackmult = trackmult;
- fEvent->n = nadded;
- }
- else if( analysisMode==kTPCTrk ){
- Sort(fTrackArrayTPCPar, kFALSE);
- }
-
- }
-
-
-}
-//__________________________________________________________________
-void AliAnalysisTaskHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent, AnalysisMode analysisMode ){
-
-
- Int_t trackmult = 0; // no pt cuts
- Int_t nadded = 0;
- Int_t nAODTracks = AODevent->GetNumberOfTracks();
- if( analysisMode == kGlobalTrk ){
-
- if(fTrackArrayGlobalPar)
- fTrackArrayGlobalPar->Clear();
-
- }
- if( analysisMode == kTPCTrk ){
- if(fTrackArrayTPCPar)
- fTrackArrayTPCPar->Clear();
- }
-
-
-
-
- //const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
- //if( vertexSPD->GetNContributors() < 1 || TMath::Abs( vertexSPD->GetZ()) > 10.0 ) return;
-
-
-
- if( analysisMode==kGlobalTrk ){
-
-
- const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
- for(Int_t iT = 0; iT < nAODTracks; iT++) {
-
- AliAODTrack* aodTrack = AODevent->GetTrack(iT);
-
-
-
- UShort_t filterFlag = 0;
- Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
- Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
- Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
-
-
-
- if (fTrackFilterGolden) {
-
- // ITSTPC2010 cuts is bit 32 according to above macro, new versions of aliroot includes the golden cuts
- if(aodTrack->TestFilterBit(32)) {
- filterFlag +=1;
- filterCut_Set3 = kTRUE;
- }
- }
-
-
- if (fTrackFilterTPC) {
- // TPC only cuts is bit 1 according to above macro
- // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
- if(aodTrack->TestFilterBit(1)){
- filterFlag +=2;
- filterCut_Set1 = kTRUE;
- trackmult++;
-
- }
- }
-
-
-
- if(filterFlag==0)
- continue;
-
-
- Double_t b[2], cov[3];
- if (!(aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov)))
- filterFlag = 32; // propagation failed!!!!!;
- Float_t dcaxy = b[0];
- Float_t dcaz = b[1];
-
-
-
-
-
- // As I understand this routine recalculates the momentum so it should be called first!
- //Double_t b[2], cov[3];
-
-
- //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
- // filterFlag = 32; // propagation failed!!!!!
-
- if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
- continue;
- if (aodTrack->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the histograms!
- // }
-
-
-
-
- Short_t charge = aodTrack->Charge();
- Float_t pt = aodTrack->Pt();
- Float_t p = aodTrack->P();
- Float_t eta = aodTrack->Eta();
- Float_t phi = aodTrack->Phi();
- AliAODPid* aodPid = aodTrack->GetDetPid();
- Short_t ncl = -10;
- Short_t neff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
- // Short_t nclf = aodTrack->GetTPCNclsF();
- Float_t dedx = -10;
- Float_t beta = -99;
- if(aodPid) {
- ncl = aodPid->GetTPCsignalN();
- dedx = aodPid->GetTPCsignal();
- //TOF
- if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
- Double_t tof[5];
- aodPid->GetIntegratedTimes(tof);
- beta = tof[0]/aodPid->GetTOFsignal();
- }
- }
- // Float_t tpcchi = aodTrack->Chi2perNDF();
-
- // Double_t p_con[3] = {0, 0, 0};
- // aodTrack->GetConstrainedPxPyPz(p_con);
- // Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
- // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
- // Float_t pttpc = tpcParam->Pt();
- // Float_t ptpc = tpcParam->P();
-
- Float_t ptMC = 0;
- Short_t pidCode = 0; // 0 = real data / no mc track!
- Short_t primaryFlag = 0; // 0 = real data / not primary mc track
- Int_t pdgMother = 0;
-
- if(fAnalysisMC) {
-
- const Int_t label = TMath::Abs(aodTrack->GetLabel());
-
- AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
- if (mcTrack){
-
- if(mcTrack->IsPhysicalPrimary())
- primaryFlag = 1;
-
- Int_t pdgCode = mcTrack->GetPdgCode();
- pidCode = GetPidCode(pdgCode);
-
- ptMC = mcTrack->Pt();
-
- AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
- pdgMother = mother->GetPdgCode();
- }
- }
-
-
- if(fTreeOption) {
-
- DeDxTrack* track = new((*fTrackArrayGlobalPar)[nadded]) DeDxTrack();
- nadded++;
-
- track->p = p;
- track->pt = pt;
- // track->ptcon = pt_con;
- // track->tpcchi = tpcchi;
- track->eta = eta;
- track->phi = phi;
- track->q = charge;
- track->filter = filterFlag;
- track->ncl = ncl;
- track->neff = neff;
- track->dedx = dedx;
- track->beta = beta;
- track->dcaxy = dcaxy;
- track->dcaz = dcaz;
- track->pid = pidCode;
- track->primary = primaryFlag;
- track->pttrue = ptMC;
- track->mother = pdgMother;
- track->filterset1 = filterCut_Set1;
- track->filterset2 = filterCut_Set2;
- track->filterset3 = filterCut_Set3;
- }
- }//end of track loop
-
-
-
- }//end of global track analysis
-
-
-
- else if( analysisMode==kTPCTrk ){
-
- const AliAODVertex* vertexSPD= (AliAODVertex*)AODevent->GetPrimaryVertexSPD();//GetPrimaryVertex()
- if( vertexSPD->GetNContributors() < 1 || TMath::Abs(vertexSPD->GetZ()) > 10.0 ) return;
-
-
- //const AliAODVertex* vertex = AODevent->GetPrimaryVertex();
- for(Int_t iT = 0; iT < nAODTracks; iT++) {
-
- AliAODTrack* aodTrack = AODevent->GetTrack(iT);
-
-
-
- UShort_t filterFlag = 0;
- Bool_t filterCut_Set1 = kFALSE;//parameters from global tracks, with TPC cuts (filter bit =1 in AOD)
- Bool_t filterCut_Set2 = kFALSE;//parameters from global tracks, cuts tpc+its 2010 W/O golden cuts
- Bool_t filterCut_Set3 = kFALSE;//parameters from global tracks, cuts its+tpc 2010 WITH golden cuts
-
-
-
- // TPC only cuts is bit 1 according to above macro
- // Alex always uses 128, NOTE: FILTER 128 ARE TPC TRACKS (TPC PARAMETERS) CONTRAINED TO THE SPD VERTEX,
- if(aodTrack->TestFilterBit(128)){
- filterFlag +=2;
- trackmult++;
- }
-
-
- if(filterFlag==0)
- continue;
-
-
- Double_t b[2], cov[3];
- //AliAODVertex* vertex = AODevent->GetPrimaryVertex();
- //Double_t b[2] = {-99., -99.};
- //Double_t bCov[3] = {-99., -99., -99.};
- //aodTrack->PropagateToDCA(aod->GetPrimaryVertex(), aod->GetMagneticField(), 100., b, bCov);
- if (!(aodTrack->PropagateToDCA(vertexSPD, AODevent->GetMagneticField(), kVeryBig, b, cov)))
- filterFlag = 32; // propagation failed!!!!!;
- Float_t dcaxy = b[0];
- Float_t dcaz = b[1];
-
-
-
-
- // As I understand this routine recalculates the momentum so it should be called first!
- //Double_t b[2], cov[3];
-
-
- //if(!aodTrack->PropagateToDCA(vertex, AODevent->GetMagneticField(), kVeryBig, b, cov))
- // filterFlag = 32; // propagation failed!!!!!
-
- if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
- continue;
- if (aodTrack->Pt() < fMinPt) {
-
- // Keep small fraction of low pT tracks
- if(fRandom->Rndm() > fLowPtFraction)
- continue;
- } // else {
- // Here we want to add the high pt part of the histograms!
- // }
-
-
-
-
- Short_t charge = aodTrack->Charge();
- Float_t pt = aodTrack->Pt();
- Float_t p = aodTrack->P();
- Float_t eta = aodTrack->Eta();
- Float_t phi = aodTrack->Phi();
- AliAODPid* aodPid = aodTrack->GetDetPid();
- Short_t ncl = -10;
- Short_t neff = 0; // This is not yet there! Short_t(aodTrack->GetTPCClusterInfo(2, 1)); // effective track length for pT res
- // Short_t nclf = aodTrack->GetTPCNclsF();
- Float_t dedx = -10;
- Float_t beta = -99;
- if(aodPid) {
- ncl = aodPid->GetTPCsignalN();
- dedx = aodPid->GetTPCsignal();
- //TOF
- if (aodTrack->GetStatus()&AliESDtrack::kTOFpid){
- Double_t tof[5];
- aodPid->GetIntegratedTimes(tof);
- beta = tof[0]/aodPid->GetTOFsignal();
- }
- }
- // Float_t tpcchi = aodTrack->Chi2perNDF();
-
- // Double_t p_con[3] = {0, 0, 0};
- // aodTrack->GetConstrainedPxPyPz(p_con);
- // Float_t pt_con = 0; // This is not there! (Float_t)TMath::Sqrt(p_con[0]*p_con[0] + p_con[1]*p_con[1]);
- // const AliExternalTrackParam* tpcParam = aodTrack->GetTPCInnerParam();
- // Float_t pttpc = tpcParam->Pt();
- // Float_t ptpc = tpcParam->P();
-
- Float_t ptMC = 0;
- Short_t pidCode = 0; // 0 = real data / no mc track!
- Short_t primaryFlag = 0; // 0 = real data / not primary mc track
- Int_t pdgMother = 0;
-
- if(fAnalysisMC) {
-
- const Int_t label = TMath::Abs(aodTrack->GetLabel());
-
- AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
- if (mcTrack){
-
- if(mcTrack->IsPhysicalPrimary())
- primaryFlag = 1;
-
- Int_t pdgCode = mcTrack->GetPdgCode();
- pidCode = GetPidCode(pdgCode);
-
- ptMC = mcTrack->Pt();
-
- AliAODMCParticle* mother = FindPrimaryMotherAOD(mcTrack);
- pdgMother = mother->GetPdgCode();
- }
- }
-
-
- if(fTreeOption) {
-
- DeDxTrack* tracktpc = new((*fTrackArrayTPCPar)[nadded]) DeDxTrack();
- nadded++;
-
- tracktpc->p = p;
- tracktpc->pt = pt;
- // track->ptcon = pt_con;
- // track->tpcchi = tpcchi;
- tracktpc->eta = eta;
- tracktpc->phi = phi;
- tracktpc->q = charge;
- tracktpc->filter = filterFlag;
- tracktpc->ncl = ncl;
- tracktpc->neff = neff;
- tracktpc->dedx = dedx;
- tracktpc->beta = beta;
- tracktpc->dcaxy = dcaxy;
- tracktpc->dcaz = dcaz;
- tracktpc->pid = pidCode;
- tracktpc->primary = primaryFlag;
- tracktpc->pttrue = ptMC;
- tracktpc->mother = pdgMother;
- tracktpc->filterset1 = filterCut_Set1;
- tracktpc->filterset2 = filterCut_Set2;
- tracktpc->filterset3 = filterCut_Set3;
- }
- }//end of track loop
-
-
- }//end of global track analysis
-
-
- if(fTreeOption) {
-
- if( analysisMode==kGlobalTrk ){
- Sort(fTrackArrayGlobalPar, kFALSE);
-
- fEvent->trackmult = trackmult;
- fEvent->n = nadded;
-
-
- }
- if( analysisMode==kTPCTrk ){
- Sort(fTrackArrayTPCPar, kFALSE);
- }
-
- }
-
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff