-/**************************************************************************\r
-* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
-* *\r
-* Author: The ALICE Off-line Project. *\r
-* Contributors are mentioned in the code where appropriate. *\r
-* *\r
-* Permission to use, copy, modify and distribute this software and its *\r
-* documentation strictly for non-commercial purposes is hereby granted *\r
-* without fee, provided that the above copyright notice appears in all *\r
-* copies and that both the copyright notice and this permission notice *\r
-* appear in the supporting documentation. The authors make no claims *\r
-* about the suitability of this software for any purpose. It is *\r
-* provided "as is" without express or implied warranty. *\r
-**************************************************************************/\r
-\r
-#include "iostream"\r
-\r
-#include <TPDGCode.h>\r
-#include <TDatabasePDG.h>\r
-\r
-#include "TChain.h"\r
-#include "TTreeStream.h"\r
-#include "TTree.h"\r
-#include "TH1F.h"\r
-#include "TCanvas.h"\r
-#include "TList.h"\r
-#include "TObjArray.h"\r
-#include "TFile.h"\r
-#include "TMatrixD.h"\r
-#include "TRandom3.h"\r
-\r
-#include "AliHeader.h" \r
-#include "AliGenEventHeader.h" \r
-#include "AliInputEventHandler.h" \r
-#include "AliStack.h" \r
-#include "AliTrackReference.h" \r
-\r
-#include "AliPhysicsSelection.h"\r
-#include "AliAnalysisTask.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliESDEvent.h"\r
-#include "AliESDfriend.h"\r
-#include "AliMCEvent.h"\r
-#include "AliESDInputHandler.h"\r
-#include "AliESDVertex.h"\r
-#include "AliTracker.h"\r
-#include "AliGeomManager.h"\r
-\r
-#include "AliCentrality.h"\r
-#include "AliESDVZERO.h"\r
-#include "AliMultiplicity.h"\r
-\r
-#include "AliESDtrackCuts.h"\r
-#include "AliMCEventHandler.h"\r
-#include "AlidNdPt.h"\r
-#include "AlidNdPtEventCuts.h"\r
-#include "AlidNdPtAcceptanceCuts.h"\r
-\r
-#include "AlidNdPtTrackDumpTask.h"\r
-#include "AliKFParticle.h"\r
-#include "AliESDv0.h"\r
-\r
-using namespace std;\r
-\r
-ClassImp(AlidNdPtTrackDumpTask)\r
-\r
-//_____________________________________________________________________________\r
-AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
- : AliAnalysisTaskSE(name)\r
- , fESD(0)\r
- , fMC(0)\r
- , fESDfriend(0)\r
- , fOutput(0)\r
- , fPitList(0)\r
- , fUseMCInfo(kFALSE)\r
- , fUseESDfriends(kFALSE)\r
- , fdNdPtEventCuts(0)\r
- , fdNdPtAcceptanceCuts(0)\r
- , fdNdPtRecAcceptanceCuts(0)\r
- , fEsdTrackCuts(0)\r
- , fTrigger(AliTriggerAnalysis::kMB1) \r
- , fAnalysisMode(AlidNdPtHelper::kTPC) \r
- , fTreeSRedirector(0)\r
- , fCentralityEstimator(0)\r
- , fLowPtTrackDownscaligF(0)\r
- , fLowPtV0DownscaligF(0)\r
- , fProcessAll(kFALSE)\r
-{\r
- // Constructor\r
-\r
- // Define input and output slots here\r
- DefineOutput(1, TList::Class());\r
-\r
-}\r
-\r
-//_____________________________________________________________________________\r
-AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
-{\r
- if(fOutput) delete fOutput; fOutput =0; \r
- if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0; \r
-\r
- if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
- if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
- if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL; \r
- if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::Notify()\r
-{\r
- static Int_t count = 0;\r
- count++;\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(chain)\r
- {\r
- Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
- }\r
-\r
-return kTRUE;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
-{\r
- // Create histograms\r
- // Called once\r
- fOutput = new TList;\r
- fOutput->SetOwner();\r
-\r
- //\r
- // create temporary file for output tree\r
- fTreeSRedirector = new TTreeSRedirector("jotwinow_Temp_Trees.root");\r
-\r
- PostData(1, fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
-{\r
- //\r
- // Called for each event\r
- //\r
-\r
- // ESD event\r
- fESD = (AliESDEvent*) (InputEvent());\r
- if (!fESD) {\r
- Printf("ERROR: ESD event not available");\r
- return;\r
- }\r
-\r
- // MC event\r
- if(fUseMCInfo) {\r
- fMC = MCEvent();\r
- if (!fMC) {\r
- Printf("ERROR: MC event not available");\r
- return;\r
- }\r
- }\r
-\r
- if(fUseESDfriends) {\r
- fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
- if(!fESDfriend) {\r
- Printf("ERROR: ESD friends not available");\r
- }\r
- }\r
-\r
- //\r
- if(fProcessAll) { \r
- ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC\r
- ProcessV0(fESD,fMC,fESDfriend);\r
- ProcessLaser(fESD,fMC,fESDfriend);\r
- ProcessdEdx(fESD,fMC,fESDfriend);\r
- if(IsUseMCInfo())\r
- ProcessMCEff(fESD,fMC,fESDfriend);\r
- }\r
- else {\r
- Process(fESD,fMC,fESDfriend);\r
- ProcessV0(fESD,fMC,fESDfriend);\r
- ProcessLaser(fESD,fMC,fESDfriend);\r
- ProcessdEdx(fESD,fMC,fESDfriend);\r
- if(IsUseMCInfo())\r
- ProcessMCEff(fESD,fMC,fESDfriend);\r
- }\r
-\r
- // Post output data.\r
- PostData(1, fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Select real events with high-pT tracks \r
- //\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
- AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
- AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
-\r
- if(!evtCuts || !accCuts || !esdTrackCuts) {\r
- Printf("ERROR cuts not available");\r
- return;\r
- }\r
-\r
- // trigger selection\r
- Bool_t isEventTriggered = kTRUE;\r
- AliPhysicsSelection *physicsSelection = NULL;\r
- AliTriggerAnalysis* triggerAnalysis = NULL;\r
-\r
- // \r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
- if (!inputHandler)\r
- {\r
- Printf("ERROR: Could not receive input handler");\r
- return;\r
- }\r
- \r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // trigger\r
- if(evtCuts->IsTriggerRequired()) \r
- {\r
- // always MB\r
- isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
-\r
- physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
- if(!physicsSelection) return;\r
- //SetPhysicsTriggerSelection(physicsSelection);\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) return;\r
- isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
- }\r
- }\r
-\r
- // centrality determination\r
- Float_t centralityF = -1;\r
- AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
- centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // use MC information\r
- AliHeader* header = 0;\r
- AliGenEventHeader* genHeader = 0;\r
- AliStack* stack = 0;\r
- TArrayF vtxMC(3);\r
-\r
- Int_t multMCTrueTracks = 0;\r
- if(IsUseMCInfo())\r
- {\r
- //\r
- if(!mcEvent) {\r
- AliDebug(AliLog::kError, "mcEvent not available");\r
- return;\r
- }\r
- // get MC event header\r
- header = mcEvent->Header();\r
- if (!header) {\r
- AliDebug(AliLog::kError, "Header not available");\r
- return;\r
- }\r
- // MC particle stack\r
- stack = mcEvent->Stack();\r
- if (!stack) {\r
- AliDebug(AliLog::kError, "Stack not available");\r
- return;\r
- }\r
-\r
- // get MC vertex\r
- genHeader = header->GenEventHeader();\r
- if (!genHeader) {\r
- AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
- return;\r
- }\r
- genHeader->PrimaryVertex(vtxMC);\r
-\r
- // multipliticy of all MC primary tracks\r
- // in Zv, pt and eta ranges)\r
- multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
-\r
- } // end bUseMC\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- const AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
- vtxESD = esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
- vtxESD = esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- return;\r
- }\r
-\r
- if(!vtxESD) return;\r
-\r
- Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
- //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
-\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered)\r
- {\r
-\r
- //\r
- // get IR information\r
- //\r
- AliESDHeader *esdHeader = 0;\r
- esdHeader = esdEvent->GetHeader();\r
- if(!esdHeader) return;\r
- //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s\r
- //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval\r
-\r
- // Use when Peter commit the changes in the header\r
- Int_t ir1 = 0;\r
- Int_t ir2 = 0;\r
-\r
- //\r
- Double_t vert[3] = {0}; \r
- vert[0] = vtxESD->GetXv();\r
- vert[1] = vtxESD->GetYv();\r
- vert[2] = vtxESD->GetZv();\r
- Int_t mult = vtxESD->GetNContributors();\r
- Double_t bz = esdEvent->GetMagneticField();\r
- Double_t runNumber = esdEvent->GetRunNumber();\r
- Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
-\r
- // high pT tracks\r
- for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
- {\r
- AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
- if(!track) continue;\r
- if(track->Charge()==0) continue;\r
- if(!esdTrackCuts->AcceptTrack(track)) continue;\r
- if(!accCuts->AcceptTrack(track)) continue;\r
- \r
- // downscale low-pT tracks\r
- Double_t scalempt= TMath::Min(track->Pt(),10.);\r
- Double_t downscaleF = gRandom->Rndm();\r
- downscaleF *= fLowPtTrackDownscaligF;\r
- if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
- //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
-\r
- AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
- if (!tpcInner) continue;\r
- // transform to the track reference frame \r
- Bool_t isOK = kFALSE;\r
- isOK = tpcInner->Rotate(track->GetAlpha());\r
- isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
- if(!isOK) continue;\r
-\r
- // Dump to the tree \r
- // vertex\r
- // TPC-ITS tracks\r
- //\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "Bz="<<bz<<\r
- "vertX="<<vert[0]<<\r
- "vertY="<<vert[1]<<\r
- "vertZ="<<vert[2]<<\r
- "IRtot="<<ir1<<\r
- "IRint2="<<ir2<<\r
- "mult="<<mult<<\r
- "esdTrack.="<<track<<\r
- "centralityF="<<centralityF<<\r
- "\n";\r
- }\r
- }\r
- \r
- PostData(1, fOutput);\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Process laser events\r
- //\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // laser events \r
- const AliESDHeader* esdHeader = esdEvent->GetHeader();\r
- if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) \r
- {\r
- Int_t countLaserTracks = 0;\r
- for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
- {\r
- AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
- if(!track) continue;\r
-\r
- if(track->GetTPCInnerParam()) countLaserTracks++;\r
- }\r
- \r
- if(countLaserTracks > 100) \r
- { \r
- Double_t runNumber = esdEvent->GetRunNumber();\r
- Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
- Double_t bz = esdEvent->GetMagneticField();\r
-\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"Laser"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "Bz="<<bz<<\r
- "multTPCtracks="<<countLaserTracks<<\r
- "\n";\r
- }\r
- }\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
-{\r
- //\r
- // Select real events with high-pT tracks\r
- // Calculate and stor in the output tree:\r
- // TPC constrained tracks\r
- // InnerParams constrained tracks\r
- // TPC-ITS tracks\r
- // ITSout-InnerParams tracks\r
- // chi2 distance between TPC constrained and TPC-ITS tracks\r
- // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
- // chi2 distance between ITSout and InnerParams tracks\r
- // MC information: \r
- // track references at ITSin, TPCin; InnerParam at first TPC track reference, \r
- // particle ID, mother ID, production mechanism ...\r
- // \r
-\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
- AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
- AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
-\r
- if(!evtCuts || !accCuts || !esdTrackCuts) {\r
- AliDebug(AliLog::kError, "cuts not available");\r
- return;\r
- }\r
-\r
- // trigger selection\r
- Bool_t isEventTriggered = kTRUE;\r
- AliPhysicsSelection *physicsSelection = NULL;\r
- AliTriggerAnalysis* triggerAnalysis = NULL;\r
-\r
- // \r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
- if (!inputHandler)\r
- {\r
- Printf("ERROR: Could not receive input handler");\r
- return;\r
- }\r
- \r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // trigger\r
- if(evtCuts->IsTriggerRequired()) \r
- {\r
- // always MB\r
- isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
-\r
- physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
- if(!physicsSelection) return;\r
- //SetPhysicsTriggerSelection(physicsSelection);\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) return;\r
- isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
- }\r
- }\r
-\r
- // centrality determination\r
- Float_t centralityF = -1;\r
- AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
- centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // use MC information\r
- AliHeader* header = 0;\r
- AliGenEventHeader* genHeader = 0;\r
- AliStack* stack = 0;\r
- TArrayF vtxMC(3);\r
-\r
- Int_t multMCTrueTracks = 0;\r
- if(IsUseMCInfo())\r
- {\r
- //\r
- if(!mcEvent) {\r
- AliDebug(AliLog::kError, "mcEvent not available");\r
- return;\r
- }\r
- // get MC event header\r
- header = mcEvent->Header();\r
- if (!header) {\r
- AliDebug(AliLog::kError, "Header not available");\r
- return;\r
- }\r
- // MC particle stack\r
- stack = mcEvent->Stack();\r
- if (!stack) {\r
- AliDebug(AliLog::kError, "Stack not available");\r
- return;\r
- }\r
-\r
- // get MC vertex\r
- genHeader = header->GenEventHeader();\r
- if (!genHeader) {\r
- AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
- return;\r
- }\r
- genHeader->PrimaryVertex(vtxMC);\r
-\r
- // multipliticy of all MC primary tracks\r
- // in Zv, pt and eta ranges)\r
- multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
-\r
- } // end bUseMC\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- const AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
- vtxESD = esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
- vtxESD = esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- return;\r
- }\r
-\r
- if(!vtxESD) return;\r
-\r
- Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
- //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
-\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered)\r
- {\r
- //\r
- // get IR information\r
- //\r
- AliESDHeader *esdHeader = 0;\r
- esdHeader = esdEvent->GetHeader();\r
- if(!esdHeader) return;\r
- //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s\r
- //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval\r
- //\r
- Int_t ir1 = 0;\r
- Int_t ir2 = 0;\r
-\r
- //\r
- Double_t vert[3] = {0}; \r
- vert[0] = vtxESD->GetXv();\r
- vert[1] = vtxESD->GetYv();\r
- vert[2] = vtxESD->GetZv();\r
- Int_t mult = vtxESD->GetNContributors();\r
- Double_t bz = esdEvent->GetMagneticField();\r
- Double_t runNumber = esdEvent->GetRunNumber();\r
- Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
-\r
- // high pT tracks\r
- for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
- {\r
- AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
- if(!track) continue;\r
- if(track->Charge()==0) continue;\r
- if(!esdTrackCuts->AcceptTrack(track)) continue;\r
- if(!accCuts->AcceptTrack(track)) continue;\r
- \r
- // downscale low-pT tracks\r
- Double_t scalempt= TMath::Min(track->Pt(),10.);\r
- Double_t downscaleF = gRandom->Rndm();\r
- downscaleF *= fLowPtTrackDownscaligF;\r
- if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
- //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
-\r
- // Dump to the tree \r
- // vertex\r
- // TPC constrained tracks\r
- // InnerParams constrained tracks\r
- // TPC-ITS tracks\r
- // ITSout-InnerParams tracks\r
- // chi2 distance between TPC constrained and TPC-ITS tracks\r
- // chi2 distance between TPC refitted constrained and TPC-ITS tracks\r
- // chi2 distance between ITSout and InnerParams tracks\r
- // MC information\r
- \r
- Double_t x[3]; track->GetXYZ(x);\r
- Double_t b[3]; AliTracker::GetBxByBz(x,b);\r
-\r
- //\r
- // Transform TPC inner params to track reference frame\r
- //\r
- Bool_t isOKtpcInner = kFALSE;\r
- AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
- if (tpcInner) {\r
- // transform to the track reference frame \r
- isOKtpcInner = tpcInner->Rotate(track->GetAlpha());\r
- isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
- }\r
-\r
- //\r
- // Constrain TPC inner to vertex\r
- // clone TPCinner has to be deleted\r
- //\r
- Bool_t isOKtpcInnerC = kFALSE;\r
- AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
- if (tpcInnerC) {\r
- isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
- isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());\r
- isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
- }\r
-\r
- //\r
- // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex \r
- // Clone track InnerParams has to be deleted\r
- //\r
- Bool_t isOKtrackInnerC = kFALSE;\r
- AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));\r
- if (trackInnerC) {\r
- isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
- isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());\r
- isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
- } \r
- \r
- //\r
- // calculate chi2 between vi and vj vectors\r
- // with covi and covj covariance matrices\r
- // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
- //\r
- TMatrixD deltaT(5,1), deltaTtrackC(5,1);\r
- TMatrixD delta(1,5), deltatrackC(1,5);\r
- TMatrixD covarM(5,5), covarMtrackC(5,5);\r
- TMatrixD chi2(1,1);\r
- TMatrixD chi2trackC(1,1);\r
-\r
- if(isOKtpcInnerC && isOKtrackInnerC) \r
- {\r
- for (Int_t ipar=0; ipar<5; ipar++) {\r
- deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
- delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
-\r
- deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
- deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];\r
-\r
- for (Int_t jpar=0; jpar<5; jpar++) {\r
- Int_t index=track->GetIndex(ipar,jpar);\r
- covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];\r
- covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];\r
- }\r
- }\r
-\r
- // chi2 distance TPC constrained and TPC+ITS\r
- TMatrixD covarMInv = covarM.Invert();\r
- TMatrixD mat2 = covarMInv*deltaT;\r
- chi2 = delta*mat2; \r
- //chi2.Print();\r
-\r
- // chi2 distance TPC refitted constrained and TPC+ITS\r
- TMatrixD covarMInvtrackC = covarMtrackC.Invert();\r
- TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;\r
- chi2trackC = deltatrackC*mat2trackC; \r
- //chi2trackC.Print();\r
- }\r
-\r
-\r
- //\r
- // Propagate ITSout to TPC inner wall \r
- // and calculate chi2 distance to track (InnerParams)\r
- //\r
- const Double_t kTPCRadius=85; \r
- const Double_t kStep=3; \r
-\r
- // clone track InnerParams has to be deleted\r
- Bool_t isOKtrackInnerC2 = kFALSE;\r
- AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
- if (trackInnerC2) {\r
- isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
- }\r
-\r
- Bool_t isOKouterITSc = kFALSE;\r
- AliExternalTrackParam *outerITSc = NULL;\r
- TMatrixD chi2OuterITS(1,1);\r
-\r
- if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
- {\r
- // propagate ITSout to TPC inner wall\r
- AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
-\r
- if(friendTrack) \r
- {\r
- outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut());\r
- if(outerITSc) \r
- {\r
- isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);\r
- isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
- isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
-\r
- //\r
- // calculate chi2 between outerITS and innerParams\r
- // cov without z-coordinate at the moment\r
- //\r
- TMatrixD deltaTouterITS(4,1);\r
- TMatrixD deltaouterITS(1,4);\r
- TMatrixD covarMouterITS(4,4);\r
-\r
- if(isOKtrackInnerC2 && isOKouterITSc) {\r
- Int_t kipar = 0;\r
- Int_t kjpar = 0;\r
- for (Int_t ipar=0; ipar<5; ipar++) {\r
- if(ipar!=1) {\r
- deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
- deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
- }\r
-\r
- kjpar=0;\r
- for (Int_t jpar=0; jpar<5; jpar++) {\r
- Int_t index=outerITSc->GetIndex(ipar,jpar);\r
- if(ipar !=1 || jpar!=1) {\r
- covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
- }\r
- if(jpar!=1) kjpar++;\r
- }\r
- if(ipar!=1) kipar++;\r
- }\r
-\r
- // chi2 distance ITSout and InnerParams\r
- TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
- TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
- chi2OuterITS = deltaouterITS*mat2outerITS; \r
- //chi2OuterITS.Print();\r
- } \r
- }\r
- }\r
- }\r
-\r
- //\r
- // MC info\r
- //\r
- TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
- TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
- Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
- Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
- Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
- Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
- Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
-\r
- AliTrackReference *refTPCIn = NULL;\r
- AliTrackReference *refITS = NULL;\r
-\r
- Bool_t isOKtrackInnerC3 = kFALSE;\r
- AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
-\r
- if(IsUseMCInfo()) \r
- {\r
- if(!stack) return;\r
-\r
- //\r
- // global track\r
- //\r
- Int_t label = TMath::Abs(track->GetLabel()); \r
- particle = stack->Particle(label);\r
- if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)\r
- {\r
- particleMother = GetMother(particle,stack);\r
- mech = particle->GetUniqueID();\r
- isPrim = stack->IsPhysicalPrimary(label);\r
- isFromStrangess = IsFromStrangeness(label,stack);\r
- isFromConversion = IsFromConversion(label,stack);\r
- isFromMaterial = IsFromMaterial(label,stack);\r
- }\r
-\r
- //\r
- // TPC track\r
- //\r
- Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
- particleTPC = stack->Particle(labelTPC);\r
- if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
- {\r
- particleMotherTPC = GetMother(particleTPC,stack);\r
- mechTPC = particleTPC->GetUniqueID();\r
- isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
- isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);\r
- isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
- isFromMaterialTPC = IsFromMaterial(labelTPC,stack);\r
- }\r
-\r
- //\r
- // store first track reference\r
- // for TPC track\r
- //\r
- TParticle *part=0;\r
- TClonesArray *trefs=0;\r
- Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
-\r
- if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
- {\r
- Int_t nTrackRef = trefs->GetEntries();\r
- //printf("nTrackRef %d \n",nTrackRef);\r
-\r
- Int_t countITS = 0;\r
- for (Int_t iref = 0; iref < nTrackRef; iref++) \r
- {\r
- AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
-\r
- // Ref. in the middle ITS \r
- if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
- {\r
- if(!refITS && countITS==2) {\r
- refITS = ref;\r
- //printf("refITS %p \n",refITS);\r
- }\r
- countITS++;\r
- }\r
-\r
- // TPC\r
- if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
- {\r
- if(!refTPCIn) {\r
- refTPCIn = ref;\r
- //printf("refTPCIn %p \n",refTPCIn);\r
- //break;\r
- }\r
- }\r
- }\r
-\r
- // transform inner params to TrackRef\r
- // reference frame\r
- if(refTPCIn && trackInnerC3) \r
- {\r
- Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
- isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);\r
- isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
- }\r
- }\r
-\r
- //\r
- // ITS track\r
- //\r
- Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
- particleITS = stack->Particle(labelITS);\r
- if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
- {\r
- particleMotherITS = GetMother(particleITS,stack);\r
- mechITS = particleITS->GetUniqueID();\r
- isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
- isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
- isFromConversionITS = IsFromConversion(labelITS,stack);\r
- isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
- }\r
- }\r
-\r
- //\r
- Bool_t dumpToTree=kFALSE;\r
- \r
- if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;\r
- if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;\r
- if(fUseMCInfo && isOKtrackInnerC3) dumpToTree = kTRUE;\r
-\r
- //\r
- if(fTreeSRedirector && dumpToTree) \r
- {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "Bz="<<bz<<\r
- "vertX="<<vert[0]<<\r
- "vertY="<<vert[1]<<\r
- "vertZ="<<vert[2]<<\r
- "IRtot="<<ir1<<\r
- "IRint2="<<ir2<<\r
- "mult="<<mult<<\r
- "esdTrack.="<<track<<\r
- "extTPCInnerC.="<<tpcInnerC<<\r
- "extInnerParamC.="<<trackInnerC<<\r
- "extInnerParam.="<<trackInnerC2<<\r
- "extOuterITS.="<<outerITSc<<\r
- "extInnerParamRef.="<<trackInnerC3<<\r
- "refTPCIn.="<<refTPCIn<<\r
- "refITS.="<<refITS<<\r
- "chi2TPCInnerC="<<chi2(0,0)<<\r
- "chi2InnerC="<<chi2trackC(0,0)<<\r
- "chi2OuterITS="<<chi2OuterITS(0,0)<<\r
- "centralityF="<<centralityF<<\r
- "particle.="<<particle<<\r
- "particleMother.="<<particleMother<<\r
- "mech="<<mech<<\r
- "isPrim="<<isPrim<<\r
- "isFromStrangess="<<isFromStrangess<<\r
- "isFromConversion="<<isFromConversion<<\r
- "isFromMaterial="<<isFromMaterial<<\r
- "particleTPC.="<<particleTPC<<\r
- "particleMotherTPC.="<<particleMotherTPC<<\r
- "mechTPC="<<mechTPC<<\r
- "isPrimTPC="<<isPrimTPC<<\r
- "isFromStrangessTPC="<<isFromStrangessTPC<<\r
- "isFromConversionTPC="<<isFromConversionTPC<<\r
- "isFromMaterialTPC="<<isFromMaterialTPC<<\r
- "particleITS.="<<particleITS<<\r
- "particleMotherITS.="<<particleMotherITS<<\r
- "mechITS="<<mechITS<<\r
- "isPrimITS="<<isPrimITS<<\r
- "isFromStrangessITS="<<isFromStrangessITS<<\r
- "isFromConversionITS="<<isFromConversionITS<<\r
- "isFromMaterialITS="<<isFromMaterialITS<<\r
- "\n";\r
- }\r
- \r
- if(tpcInnerC) delete tpcInnerC;\r
- if(trackInnerC) delete trackInnerC;\r
- if(trackInnerC2) delete trackInnerC2;\r
- if(outerITSc) delete outerITSc;\r
- if(trackInnerC3) delete trackInnerC3;\r
- }\r
- }\r
- \r
- PostData(1, fOutput);\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Fill tree for efficiency studies MC only\r
-\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- if(!mcEvent) {\r
- AliDebug(AliLog::kError, "mcEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
- AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
- AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
-\r
- if(!evtCuts || !accCuts || !esdTrackCuts) {\r
- AliDebug(AliLog::kError, "cuts not available");\r
- return;\r
- }\r
-\r
- // trigger selection\r
- Bool_t isEventTriggered = kTRUE;\r
- AliPhysicsSelection *physicsSelection = NULL;\r
- AliTriggerAnalysis* triggerAnalysis = NULL;\r
-\r
- // \r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
- if (!inputHandler)\r
- {\r
- Printf("ERROR: Could not receive input handler");\r
- return;\r
- }\r
- \r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // trigger\r
- if(evtCuts->IsTriggerRequired()) \r
- {\r
- // always MB\r
- isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
-\r
- physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
- if(!physicsSelection) return;\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) return;\r
- isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
- }\r
- }\r
-\r
- // centrality determination\r
- Float_t centralityF = -1;\r
- AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
- centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
- // use MC information\r
- AliHeader* header = 0;\r
- AliGenEventHeader* genHeader = 0;\r
- AliStack* stack = 0;\r
- TArrayF vtxMC(3);\r
-\r
- Int_t multMCTrueTracks = 0;\r
- if(IsUseMCInfo())\r
- {\r
- //\r
- if(!mcEvent) {\r
- AliDebug(AliLog::kError, "mcEvent not available");\r
- return;\r
- }\r
- // get MC event header\r
- header = mcEvent->Header();\r
- if (!header) {\r
- AliDebug(AliLog::kError, "Header not available");\r
- return;\r
- }\r
- // MC particle stack\r
- stack = mcEvent->Stack();\r
- if (!stack) {\r
- AliDebug(AliLog::kError, "Stack not available");\r
- return;\r
- }\r
-\r
- // get MC vertex\r
- genHeader = header->GenEventHeader();\r
- if (!genHeader) {\r
- AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
- return;\r
- }\r
- genHeader->PrimaryVertex(vtxMC);\r
-\r
- // multipliticy of all MC primary tracks\r
- // in Zv, pt and eta ranges)\r
- multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
-\r
- } // end bUseMC\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- const AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
- vtxESD = esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
- vtxESD = esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- return;\r
- }\r
-\r
- if(!vtxESD) return;\r
-\r
- Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
- //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered)\r
- {\r
- if(IsUseMCInfo()) \r
- {\r
- if(!stack) return;\r
-\r
- //\r
- // MC info\r
- //\r
- TParticle *particle=NULL;\r
- TParticle *particleMother=NULL;\r
- Int_t mech=-1;\r
-\r
- // reco event info\r
- Double_t vert[3] = {0}; \r
- vert[0] = vtxESD->GetXv();\r
- vert[1] = vtxESD->GetYv();\r
- vert[2] = vtxESD->GetZv();\r
- Int_t mult = vtxESD->GetNContributors();\r
- Double_t bz = esdEvent->GetMagneticField();\r
- Double_t runNumber = esdEvent->GetRunNumber();\r
- Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
-\r
- // loop over MC stack\r
- for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) \r
- {\r
- particle = stack->Particle(iMc);\r
- if (!particle)\r
- continue;\r
-\r
- // only charged particles\r
- if(!particle->GetPDG()) continue;\r
- Double_t charge = particle->GetPDG()->Charge()/3.;\r
- if (TMath::Abs(charge) < 0.001)\r
- continue;\r
-\r
- // only primary particles\r
- Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
- if(!prim) continue;\r
-\r
- // downscale low-pT particles\r
- Double_t scalempt= TMath::Min(particle->Pt(),10.);\r
- Double_t downscaleF = gRandom->Rndm();\r
- downscaleF *= fLowPtTrackDownscaligF;\r
- if(TMath::Exp(2*scalempt)<downscaleF) continue;\r
-\r
- // is particle in acceptance\r
- if(!accCuts->AcceptTrack(particle)) continue;\r
- \r
- // check if particle reconstructed\r
- Bool_t isRec = kFALSE;\r
- Int_t trackIndex = -1;\r
- for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
- {\r
- \r
- AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
- if(!track) continue;\r
- if(track->Charge()==0) continue;\r
- if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
- {\r
- Int_t label = TMath::Abs(track->GetLabel());\r
- if(label == iMc) {\r
- isRec = kTRUE;\r
- trackIndex = iTrack;\r
- break;\r
- }\r
- } \r
- }\r
-\r
- // Store information in the output tree\r
- AliESDtrack *recTrack = NULL; \r
- if(trackIndex>-1) { \r
- recTrack = esdEvent->GetTrack(trackIndex); \r
- } else {\r
- recTrack = new AliESDtrack(); \r
- } \r
-\r
- particleMother = GetMother(particle,stack);\r
- mech = particle->GetUniqueID();\r
-\r
- //MC particle track length\r
- Double_t tpcTrackLength = 0.;\r
- AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);\r
- if(mcParticle) {\r
- Int_t counter;\r
- tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);\r
- } \r
-\r
-\r
- //\r
- if(fTreeSRedirector) {\r
- (*fTreeSRedirector)<<"MCEffTree"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "Bz="<<bz<<\r
- "vertX="<<vert[0]<<\r
- "vertY="<<vert[1]<<\r
- "vertZ="<<vert[2]<<\r
- "mult="<<mult<<\r
- "esdTrack.="<<recTrack<<\r
- "isRec="<<isRec<<\r
- "tpcTrackLength="<<tpcTrackLength<<\r
- "particle.="<<particle<<\r
- "particleMother.="<<particleMother<<\r
- "mech="<<mech<<\r
- "\n";\r
- }\r
-\r
- if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
- }\r
- }\r
- }\r
- \r
- PostData(1, fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::IsHighDeDxParticle(AliESDtrack * track) {\r
- //\r
- // check if particle is Z > 1 \r
- //\r
- if (track->GetTPCNcls() < 60) return kFALSE;\r
- Double_t mom = track->GetInnerParam()->GetP();\r
- if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization\r
- Float_t dca[2], bCov[3];\r
- track->GetImpactParameters(dca,bCov);\r
- //\r
-\r
- Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);\r
-\r
- if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;\r
-\r
- return kFALSE;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Select real events with V0 (K0s and Lambda) high-pT candidates\r
- //\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
- AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
- AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
-\r
- if(!evtCuts || !accCuts || !esdTrackCuts) {\r
- AliDebug(AliLog::kError, "cuts not available");\r
- return;\r
- }\r
-\r
- // trigger selection\r
- Bool_t isEventTriggered = kTRUE;\r
- AliPhysicsSelection *physicsSelection = NULL;\r
- AliTriggerAnalysis* triggerAnalysis = NULL;\r
-\r
- // \r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
- if (!inputHandler)\r
- {\r
- Printf("ERROR: Could not receive input handler");\r
- return;\r
- }\r
- \r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // trigger\r
- if(evtCuts->IsTriggerRequired()) \r
- {\r
- // always MB\r
- isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
-\r
- physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
- if(!physicsSelection) return;\r
- //SetPhysicsTriggerSelection(physicsSelection);\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) return;\r
- isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
- }\r
- }\r
-\r
- // centrality determination\r
- Float_t centralityF = -1;\r
- AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
- centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
-\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- const AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
- vtxESD = esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
- vtxESD = esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- return;\r
- }\r
-\r
- if(!vtxESD) return;\r
-\r
- Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
- //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered) {\r
- //\r
- // Dump the pt downscaled V0 into the tree\r
- // \r
- Int_t ntracks = esdEvent->GetNumberOfTracks();\r
- Int_t nV0s = esdEvent->GetNumberOfV0s();\r
- Int_t run = esdEvent->GetRunNumber();\r
- Int_t time= esdEvent->GetTimeStamp();\r
- Int_t evNr=esdEvent->GetEventNumberInFile();\r
- Double_t bz = esdEvent->GetMagneticField();\r
-\r
-\r
- for (Int_t iv0=0; iv0<nV0s; iv0++){\r
- AliESDv0 * v0 = esdEvent->GetV0(iv0);\r
- if (!v0) continue;\r
- AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));\r
- AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));\r
- if (!track0) continue;\r
- if (!track1) continue;\r
- if (track0->GetSign()<0) {\r
- track1 = esdEvent->GetTrack(v0->GetIndex(0));\r
- track0 = esdEvent->GetTrack(v0->GetIndex(1));\r
- }\r
- //\r
- Bool_t isDownscaled = IsV0Downscaled(v0);\r
- if (isDownscaled) continue;\r
- AliKFParticle kfparticle; //\r
- Int_t type=GetKFParticle(v0,esdEvent,kfparticle);\r
- if (type==0) continue; \r
-\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"V0s"<<\r
- "isDownscaled="<<isDownscaled<<\r
- "Bz="<<bz<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<run<<\r
- "evtTimeStamp="<<time<<\r
- "evtNumberInFile="<<evNr<<\r
- "type="<<type<<\r
- "ntracks="<<ntracks<<\r
- "v0.="<<v0<<\r
- "kf.="<<&kfparticle<<\r
- "track0.="<<track0<<\r
- "track1.="<<track1<<\r
- "centralityF="<<centralityF<<\r
- "\n";\r
- }\r
- }\r
- PostData(1, fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Select real events with large TPC dEdx signal\r
- //\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
- AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
- AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
-\r
- if(!evtCuts || !accCuts || !esdTrackCuts) {\r
- AliDebug(AliLog::kError, "cuts not available");\r
- return;\r
- }\r
-\r
- // get file name\r
- TTree *chain = (TChain*)GetInputData(0);\r
- if(!chain) { \r
- Printf("ERROR: Could not receive input chain");\r
- return;\r
- }\r
- TObjString fileName(chain->GetCurrentFile()->GetName());\r
-\r
- // trigger\r
- Bool_t isEventTriggered = kTRUE;\r
- AliPhysicsSelection *physicsSelection = NULL;\r
- AliTriggerAnalysis* triggerAnalysis = NULL;\r
-\r
- // \r
- AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
- if (!inputHandler)\r
- {\r
- Printf("ERROR: Could not receive input handler");\r
- return;\r
- }\r
-\r
- if(evtCuts->IsTriggerRequired()) \r
- {\r
- // always MB\r
- isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
-\r
- physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
- if(!physicsSelection) return;\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) return;\r
- isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
- }\r
- }\r
-\r
- // get reconstructed vertex \r
- const AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
- vtxESD = esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
- vtxESD = esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- return;\r
- }\r
- if(!vtxESD) return;\r
-\r
- Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
- //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
-\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered)\r
- {\r
- Double_t vert[3] = {0}; \r
- vert[0] = vtxESD->GetXv();\r
- vert[1] = vtxESD->GetYv();\r
- vert[2] = vtxESD->GetZv();\r
- Int_t mult = vtxESD->GetNContributors();\r
- Double_t bz = esdEvent->GetMagneticField();\r
- Double_t runNumber = esdEvent->GetRunNumber();\r
- Double_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
-\r
- // large dEdx\r
- for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
- {\r
- AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
- if(!track) continue;\r
- if(track->Charge()==0) continue;\r
- if(!esdTrackCuts->AcceptTrack(track)) continue;\r
- if(!accCuts->AcceptTrack(track)) continue;\r
-\r
- if(!IsHighDeDxParticle(track)) continue;\r
- \r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"dEdx"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "Bz="<<bz<<\r
- "vertX="<<vert[0]<<\r
- "vertY="<<vert[1]<<\r
- "vertZ="<<vert[2]<<\r
- "mult="<<mult<<\r
- "esdTrack.="<<track<<\r
- "\n";\r
- }\r
- }\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Int_t AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)\r
-{\r
- //\r
- // Create KF particle in case the V0 fullfill selection criteria\r
- //\r
- // Selection criteria\r
- // 0. algorithm cut\r
- // 1. track cut\r
- // 3. chi2 cut\r
- // 4. rough mass cut\r
- // 5. Normalized pointing angle cut\r
- //\r
- const Double_t cutMass=0.2;\r
- const Double_t kSigmaDCACut=3;\r
- //\r
- // 0.) algo cut - accept only on the fly\r
- //\r
- if (v0->GetOnFlyStatus() ==kFALSE) return 0; \r
- //\r
- // 1.) track cut\r
- // \r
- AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));\r
- AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));\r
- /*\r
- TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";\r
- TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";\r
- TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";\r
- */ \r
- if (TMath::Abs(track0->GetTgl())>1) return 0;\r
- if (TMath::Abs(track1->GetTgl())>1) return 0;\r
- if ((track0->GetTPCClusterInfo(2,1))<100) return 0;\r
- if ((track1->GetTPCClusterInfo(2,1))<100) return 0;\r
- //if ((track0->GetITSclusters(0))<2) return 0;\r
- //if ((track1->GetITSclusters(0))<2) return 0; \r
- Float_t pos0[2]={0}, cov0[3]={0};\r
- Float_t pos1[2]={0}, cov1[3]={0};\r
- track0->GetImpactParameters(pos0,cov0);\r
- track0->GetImpactParameters(pos1,cov1);\r
- //\r
- if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;\r
- if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;\r
- // \r
- //\r
- // 3.) Chi2 cut\r
- //\r
- Double_t chi2KF = v0->GetKFInfo(2,2,2);\r
- if (chi2KF>25) return 0;\r
- //\r
- // 4.) Rough mass cut - 0.200 GeV\r
- //\r
- static Double_t masses[2]={-1};\r
- if (masses[0]<0){\r
- masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();\r
- masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();\r
- }\r
- Double_t mass00= v0->GetEffMass(0,0);\r
- Double_t mass22= v0->GetEffMass(2,2);\r
- Double_t mass42= v0->GetEffMass(4,2);\r
- Double_t mass24= v0->GetEffMass(2,4);\r
- Bool_t massOK=kFALSE;\r
- Int_t type=0;\r
- Int_t ptype=0;\r
- Double_t dmass=1;\r
- Int_t p1=0, p2=0;\r
- if (TMath::Abs(mass00-0)<cutMass) {\r
- massOK=kTRUE; type+=1; \r
- if (TMath::Abs(mass00-0)<dmass) {\r
- ptype=1;\r
- dmass=TMath::Abs(mass00-0); \r
- p1=0; p2=0;\r
- } \r
- }\r
- if (TMath::Abs(mass24-masses[1])<cutMass) {\r
- massOK=kTRUE; type+=2; \r
- if (TMath::Abs(mass24-masses[1])<dmass){\r
- dmass = TMath::Abs(mass24-masses[1]);\r
- ptype=2;\r
- p1=2; p2=4;\r
- }\r
- }\r
- if (TMath::Abs(mass42-masses[1])<cutMass) {\r
- massOK=kTRUE; type+=4;\r
- if (TMath::Abs(mass42-masses[1])<dmass){\r
- dmass = TMath::Abs(mass42-masses[1]);\r
- ptype=4;\r
- p1=4; p2=2;\r
- }\r
- }\r
- if (TMath::Abs(mass22-masses[0])<cutMass) {\r
- massOK=kTRUE; type+=8;\r
- if (TMath::Abs(mass22-masses[0])<dmass){\r
- dmass = TMath::Abs(mass22-masses[0]);\r
- ptype=8;\r
- p1=2; p2=2;\r
- }\r
- }\r
- if (type==0) return 0;\r
- //\r
- const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};\r
- const AliExternalTrackParam *paramP = v0->GetParamP();\r
- const AliExternalTrackParam *paramN = v0->GetParamN();\r
- if (paramP->GetSign()<0){\r
- paramP=v0->GetParamP();\r
- paramN=v0->GetParamN();\r
- }\r
- //Double_t *pparam1 = (Double_t*)paramP->GetParameter();\r
- //Double_t *pparam2 = (Double_t*)paramN->GetParameter();\r
- //\r
- AliKFParticle kfp1( *paramP, spdg[p1] );\r
- AliKFParticle kfp2( *paramN, -1 *spdg[p2] );\r
- AliKFParticle V0KF;\r
- (V0KF)+=kfp1;\r
- (V0KF)+=kfp2;\r
- kfparticle=V0KF;\r
- //\r
- // Pointing angle\r
- //\r
- Double_t errPhi = V0KF.GetErrPhi();\r
- Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());\r
- if (pointAngle/errPhi>10) return 0; \r
- //\r
- return ptype; \r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)\r
-{\r
- //\r
- // Downscale randomly low pt V0\r
- //\r
- //return kFALSE;\r
- Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());\r
- Double_t scalempt= TMath::Min(maxPt,10.);\r
- Double_t downscaleF = gRandom->Rndm();\r
- downscaleF *= fLowPtV0DownscaligF;\r
- \r
- //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);\r
- if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;\r
- return kFALSE;\r
-\r
- /*\r
- TH1F his1("his1","his1",100,0,10);\r
- TH1F his2("his2","his2",100,0,10);\r
- {for (Int_t i=0; i<10000; i++){\r
- Double_t rnd=gRandom->Exp(1);\r
- Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();\r
- his1->Fill(rnd); \r
- if (!isDownscaled) his2->Fill(rnd); \r
- }}\r
-\r
- */\r
-\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-Bool_t AlidNdPtTrackDumpTask::MakeTPCInnerC(AliESDtrack *const track, const AliESDVertex* vtx, Double_t b[3])\r
-{\r
-//\r
-// return TPC inner constrain object (must be deleted)\r
-//\r
-\r
-if(!track) return NULL;\r
-if(!vtx) return NULL;\r
- \r
- AliExternalTrackParam * tpcInnerC = NULL;\r
- Bool_t isOK = kFALSE;\r
-\r
- tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));\r
- if (!tpcInnerC) return NULL;\r
- isOK = ConstrainTPCInner(tpcInnerC,vtx,b);\r
- isOK = tpcInnerC->Rotate(track->GetAlpha());\r
- isOK = tpcInnerC->PropagateTo(track->GetX(),b[2]);\r
- if(!isOK) {\r
- delete tpcInnerC;\r
- return NULL;\r
- }\r
-\r
- if(fTreeSRedirector) {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "esdTrack.="<<track<<\r
- "extTPCInnerC.="<<tpcInnerC<<\r
- "chi2TPCInnerC="<<chi2TPCInnerC;\r
- }\r
-\r
-return tpcInnerC;\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
-{\r
- // Constrain TPC inner params constrained\r
- //\r
- if(!tpcInnerC) return kFALSE; \r
- if(!vtx) return kFALSE;\r
-\r
- Double_t dz[2],cov[3];\r
- //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
- //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
- //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
- if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
-\r
-\r
- Double_t covar[6]; vtx->GetCovMatrix(covar);\r
- Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
- Double_t c[3]={covar[2],0.,covar[5]};\r
- Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
- if (chi2C>kVeryBig) return kFALSE; \r
-\r
- if(!tpcInnerC->Update(p,c)) return kFALSE;\r
-\r
- return kTRUE;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-AliExternalTrackParam * AlidNdPtTrackDumpTask::CalculateChi2(AliESDtrack *const track, AliExternalTrackParam *const trackParam)\r
-{\r
-//\r
-// calculate chi2 between vi and vj vectors\r
-// with covi and covj covariance matrices\r
-// chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
-//\r
-\r
-if(!track) return 0;\r
-if(!trackParam) return 0;\r
-\r
- TMatrixD deltaT(5,1);\r
- TMatrixD delta(1,5);\r
- TMatrixD covarM(5,5);\r
-\r
- for (Int_t ipar=0; ipar<5; ipar++) {\r
- deltaT(ipar,0)=trackParam->GetParameter()[ipar]-track->GetParameter()[ipar];\r
- delta(0,ipar)=trackParam->GetParameter()[ipar]-track->GetParameter()[ipar];\r
- for (Int_t jpar=0; jpar<5; jpar++) {\r
- Int_t index=track->GetIndex(ipar,jpar);\r
- covarM(ipar,jpar)=track->GetCovariance()[index]+trackParam->GetCovariance()[index];\r
- }\r
- }\r
-\r
- // chi2 distance \r
- TMatrixD covarMInv = covarM.Invert();\r
- TMatrixD mat2 = covarMInv*deltaT;\r
- TMatrixD chi2 = delta*mat2; \r
- //chi2.Print();\r
-\r
-return ((Double_t)chi(0,0));\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-AliExternalTrackParam * AlidNdPtTrackDumpTask::CalculateChi2(AliExternalTrackParam *const trackParam1, AliExternalTrackParam *const trackParam2)\r
-{\r
-//\r
-// calculate chi2 between vi and vj vectors\r
-// with covi and covj covariance matrices\r
-// chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)\r
-//\r
-\r
-if(!track) return 0;\r
-if(!trackParam) return 0;\r
-\r
-TMatrixD deltaT(5,1);\r
-TMatrixD delta(1,5);\r
-TMatrixD covarM(5,5);\r
-\r
- for (Int_t ipar=0; ipar<5; ipar++) {\r
- deltaT(ipar,0)=trackParam1->GetParameter()[ipar]-trackParam2->GetParameter()[ipar];\r
- delta(0,ipar)=trackParam1->GetParameter()[ipar]-trackParam2->GetParameter()[ipar];\r
- for (Int_t jpar=0; jpar<5; jpar++) {\r
- Int_t index=track->GetIndex(ipar,jpar);\r
- covarM(ipar,jpar)=trackParam1->GetCovariance()[index]+trackParam2->GetCovariance()[index];\r
- }\r
- }\r
-\r
- // chi2 distance \r
- TMatrixD covarMInv = covarM.Invert();\r
- TMatrixD mat2 = covarMInv*deltaT;\r
- TMatrixD chi2 = delta*mat2; \r
- //chi2.Print();\r
-\r
-return ((Double_t)chi(0,0));\r
-}\r
-*/\r
-\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-AliExternalTrackParam * AlidNdPtTrackDumpTask::MakeTrackInnerC(AliESDtrack *const track, const AliESDVertex* vtx, Double_t b[3])\r
-{\r
-//\r
-// Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
-// to vertex\r
-//\r
-if(!track) return NULL;\r
-if(!vtx) return NULL;\r
-\r
- // clone track InnerParams has to be deleted\r
- AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));\r
- if (!trackInnerC) return NULL;\r
- Bool_t isOK = ConstrainTrackInner(trackInnerC,vtx,track->GetMass(),b);\r
- isOK = trackInnerC->Rotate(track->GetAlpha());\r
- isOK = trackInnerC->PropagateTo(track->GetX(),b[2]);\r
- if(!isOK) {\r
- if(trackInnerC) delete trackInnerC;\r
- return NULL;\r
- }\r
- \r
- \r
- // Dump to tree\r
- if(fTreeSRedirector) \r
- {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "extInnerParamC.="<<trackInnerC<<\r
- "chi2InnerC="<<chi2InnerC<<\r
- "\n";\r
- }\r
-\r
-return trackInnerC;\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
-{\r
- // Constrain TPC inner params constrained\r
- //\r
- if(!trackInnerC) return kFALSE; \r
- if(!vtx) return kFALSE;\r
-\r
- const Double_t kRadius = 2.8; \r
- const Double_t kMaxStep = 1.0; \r
-\r
- Double_t dz[2],cov[3];\r
-\r
- //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
- //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
- //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
-\r
- if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
- if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
-\r
- //\r
- Double_t covar[6]; vtx->GetCovMatrix(covar);\r
- Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
- Double_t c[3]={covar[2],0.,covar[5]};\r
- Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
- if (chi2C>kVeryBig) return kFALSE; \r
-\r
- if(!trackInnerC->Update(p,c)) return kFALSE;\r
-\r
- return kTRUE;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-AliExternalTrackParam * AlidNdPtTrackDumpTask::PropagateInnerParam(AliESDtrack *const track, Double_t radius, Double_t step)\r
-{\r
-//\r
-// Propagate InnerParams to inner wall \r
-//\r
-if(!track) return NULL;\r
-\r
- Bool_t isOK = kFALSE;\r
-\r
- // clone track InnerParams has to be deleted\r
- AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
- if (!trackInnerC2) return NULL; \r
- isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC2,radius,track->GetMass(),step,kFALSE);\r
- if(!isOK) {\r
- delete trackInnerC2;\r
- return NULL;\r
- }\r
-\r
-return trackInnerC2;\r
-}\r
-*/\r
-\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-AliExternalTrackParam * AlidNdPtTrackDumpTask::PropagateITSOut(Int_t iTrack, AliExternalTrackParam * trackParam, AliESDfriend *const esdFriend, Double_t b[3], Double_t radius, Double_t step)\r
-{\r
-//\r
-// Propagate ITSout to TPC inner wall \r
-//\r
-\r
-if(!track) return NULL;\r
-if(!esdFriend) return NULL;\r
-\r
- Bool_t isOK = kFALSE;\r
- AliExternalTrackParam *outerITSc = 0;\r
-\r
- if(esdFriend && esdFriend->TestSkipBit()==kFALSE) \r
- {\r
- // propagate ITSout to TPC inner wall\r
- AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);\r
- if(friendTrack) \r
- {\r
- if( (outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut())) ) \r
- {\r
- if(AliTracker::PropagateTrackToBxByBz(outerITSc,radius,track->GetMass(),step,kFALSE))\r
- {\r
- // transform ITSout to the track InnerParams reference frame \r
- isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
- isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),b[2]);\r
- if(!isOK) {\r
- if(trackInnerC2) delete trackInnerC2;\r
- if(outerITSc) delete outerITSc;\r
- return kFALSE;\r
- }\r
- \r
- //\r
- // calculate chi2 between outerITS and innerParams\r
- // cov without z-coordinate at the moment\r
- // \r
- TMatrixD deltaTouterITS(4,1);\r
- TMatrixD deltaouterITS(1,4);\r
- TMatrixD covarMouterITS(4,4);\r
-\r
- Int_t kipar = 0;\r
- Int_t kjpar = 0;\r
- for (Int_t ipar=0; ipar<5; ipar++) {\r
- if(ipar!=1) {\r
- deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
- deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
- }\r
-\r
- kjpar=0;\r
- for (Int_t jpar=0; jpar<5; jpar++) {\r
- Int_t index=outerITSc->GetIndex(ipar,jpar);\r
- if(ipar !=1 || jpar!=1) {\r
- covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
- }\r
- if(jpar!=1) kjpar++;\r
- }\r
- if(ipar!=1) kipar++;\r
- }\r
-\r
- // chi2 distance ITSout and InnerParams\r
- TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
- TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
- TMatrixD chi2OuterITS = deltaouterITS*mat2outerITS; \r
- //chi2OuterITS.Print();\r
- chi2ITSout = chi2OuterITS(0,0);\r
-\r
- // dump to the tree\r
- if(fTreeSRedirector) \r
- {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "extInnerParam.="<<trackInnerC2<<\r
- "extOuterITS.="<<outerITSc<<\r
- "chi2OuterITS="<<chi2ITSout;\r
- }\r
- } \r
- }\r
- }\r
- }\r
-\r
- if(trackInnerC2) delete trackInnerC2;\r
- if(outerITSc) delete outerITSc;\r
-\r
-return kTRUE;\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-Bool_t AlidNdPtTrackDumpTask::UseMCInfoAndDump(AliMCEvent *const mcEvent, AliESDtrack *const track, AliStack *const stack)\r
-{\r
-//\r
-// MC info\r
-//\r
-if(!mcEvent) return kFALSE;\r
-if(!track) return kFALSE;\r
-if(!stack) return kFALSE;\r
-\r
-\r
- const Double_t kStep=3; \r
-\r
-\r
- TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;\r
- TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;\r
- static Int_t mech=-1, mechTPC=-1, mechITS=-1;\r
- static Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;\r
- static Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;\r
- static Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;\r
- static Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;\r
-\r
- AliTrackReference *refTPCIn = NULL;\r
- AliTrackReference *refITS = NULL;\r
-\r
- AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
- if (!trackInnerC3) return kFALSE; \r
-\r
- //\r
- // global track\r
- //\r
- Int_t label = TMath::Abs(track->GetLabel()); \r
- particle = stack->Particle(label);\r
- if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.) \r
- {\r
- particleMother = GetMother(particle,stack);\r
- mech = particle->GetUniqueID();\r
- isPrim = stack->IsPhysicalPrimary(label);\r
- isFromStrangess = IsFromStrangeness(label,stack);\r
- isFromConversion = IsFromConversion(label,stack);\r
- isFromMaterial = IsFromMaterial(label,stack);\r
- }\r
-\r
- //\r
- // TPC track\r
- //\r
- Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); \r
- particleTPC = stack->Particle(labelTPC);\r
- if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)\r
- {\r
- particleMotherTPC = GetMother(particleTPC,stack);\r
- mechTPC = particleTPC->GetUniqueID();\r
- isPrimTPC = stack->IsPhysicalPrimary(labelTPC);\r
- isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);\r
- isFromConversionTPC = IsFromConversion(labelTPC,stack);\r
- isFromMaterialTPC = IsFromMaterial(labelTPC,stack);\r
- }\r
-\r
- //\r
- // store first track reference\r
- // for TPC track\r
- //\r
- TParticle *part=0;\r
- TClonesArray *trefs=0;\r
- Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);\r
-\r
- if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) \r
- {\r
- Int_t nTrackRef = trefs->GetEntries();\r
- //printf("nTrackRef %d \n",nTrackRef);\r
-\r
- Int_t countITS = 0;\r
- for (Int_t iref = 0; iref < nTrackRef; iref++) \r
- {\r
- AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);\r
- //printf("ref %p \n",ref);\r
- //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
- //printf("AliTrackReference::kTPC %d \n",AliTrackReference::kTPC);\r
-\r
- // Ref. in the middle ITS \r
- if(ref && ref->DetectorId()==AliTrackReference::kITS)\r
- {\r
- if(!refITS && countITS==2) {\r
- refITS = ref;\r
- //printf("refITS %p \n",refITS);\r
- }\r
- countITS++;\r
- }\r
-\r
- // TPC\r
- if(ref && ref->DetectorId()==AliTrackReference::kTPC)\r
- {\r
- if(!refTPCIn) {\r
- refTPCIn = ref;\r
- //printf("refTPCIn %p \n",refTPCIn);\r
- //break;\r
- }\r
- }\r
- }\r
-\r
- // transform inner params to TrackRef\r
- // reference frame\r
- if(refTPCIn && trackInnerC3) {\r
- Double_t fRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());\r
- Bool_t isOK = kFALSE;\r
- isOK = trackInnerC3->Rotate(fRefPhi);\r
- isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
- if(!isOK){\r
- delete trackInnerC3;\r
- return kFALSE;\r
- }\r
- }\r
- }\r
-\r
- //\r
- // ITS track\r
- //\r
- Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
- particleITS = stack->Particle(labelITS);\r
- if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
- {\r
- particleMotherITS = GetMother(particleITS,stack);\r
- mechITS = particleITS->GetUniqueID();\r
- isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
- isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
- isFromConversionITS = IsFromConversion(labelITS,stack);\r
- isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
- }\r
-\r
- // dump to tree\r
- if(fTreeSRedirector) \r
- {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "extInnerParamRef.="<<trackInnerC3<<\r
- "refTPCIn.="<<refTPCIn<<\r
- "refITS.="<<refITS<<\r
- "particle.="<<particle<<\r
- "particleMother.="<<particleMother<<\r
- "mech="<<mech<<\r
- "isPrim="<<isPrim<<\r
- "isFromStrangess="<<isFromStrangess<<\r
- "isFromConversion="<<isFromConversion<<\r
- "isFromMaterial="<<isFromMaterial<<\r
- "particleTPC.="<<particleTPC<<\r
- "particleMotherTPC.="<<particleMotherTPC<<\r
- "mechTPC="<<mechTPC<<\r
- "isPrimTPC="<<isPrimTPC<<\r
- "isFromStrangessTPC="<<isFromStrangessTPC<<\r
- "isFromConversionTPC="<<isFromConversionTPC<<\r
- "isFromMaterialTPC="<<isFromMaterialTPC<<\r
- "particleITS.="<<particleITS<<\r
- "particleMotherITS.="<<particleMotherITS<<\r
- "mechITS="<<mechITS<<\r
- "isPrimITS="<<isPrimITS<<\r
- "isFromStrangessITS="<<isFromStrangessITS<<\r
- "isFromConversionITS="<<isFromConversionITS<<\r
- "isFromMaterialITS="<<isFromMaterialITS<<\r
- "\n";\r
- }\r
-\r
- if(trackInnerC3) delete trackInnerC3;\r
-\r
-return kTRUE;\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-/*\r
-Bool_t AlidNdPtTrackDumpTask::DumpEventInfo() \r
-{\r
-//\r
-// Dump run and event information to tree\r
-// \r
- if(fTreeSRedirector) \r
- {\r
- (*fTreeSRedirector)<<"dNdPtTree"<<\r
- "\n";\r
- }\r
-\r
-return kTRUE;\r
-}\r
-*/\r
-\r
-//_____________________________________________________________________________\r
-TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
-{\r
- if(!particle) return NULL;\r
- if(!stack) return NULL;\r
-\r
- Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
- TParticle* mother = NULL; \r
- mother = stack->Particle(motherLabel); \r
-\r
-return mother;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromConversion = kFALSE;\r
-\r
- if(stack) {\r
- TParticle* particle = stack->Particle(label);\r
-\r
- if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
- {\r
- Int_t mech = particle->GetUniqueID(); // production mechanism \r
- Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
-\r
- Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
- TParticle* mother = stack->Particle(motherLabel); \r
- if(mother) {\r
- Int_t motherPdg = mother->GetPdgCode();\r
-\r
- if(!isPrim && mech==5 && motherPdg==kGamma) { \r
- isFromConversion=kTRUE; \r
- }\r
- }\r
- } \r
- } \r
-\r
- return isFromConversion;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromMaterial = kFALSE;\r
-\r
- if(stack) {\r
- TParticle* particle = stack->Particle(label);\r
-\r
- if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
- {\r
- Int_t mech = particle->GetUniqueID(); // production mechanism \r
- Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
-\r
- Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
- TParticle* mother = stack->Particle(motherLabel); \r
- if(mother) {\r
- if(!isPrim && mech==13) { \r
- isFromMaterial=kTRUE; \r
- }\r
- }\r
- } \r
- } \r
-\r
- return isFromMaterial;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromStrangeness = kFALSE;\r
-\r
- if(stack) {\r
- TParticle* particle = stack->Particle(label);\r
-\r
- if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
- {\r
- Int_t mech = particle->GetUniqueID(); // production mechanism \r
- Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
-\r
- Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
- TParticle* mother = stack->Particle(motherLabel); \r
- if(mother) {\r
- Int_t motherPdg = mother->GetPdgCode();\r
-\r
- // K+-, lambda, antilambda, K0s decays\r
- if(!isPrim && mech==4 && \r
- (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
- {\r
- isFromStrangeness = kTRUE;\r
- } \r
- }\r
- } \r
- } \r
-\r
- return isFromStrangeness;\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
-{\r
- //\r
- // Called one at the end \r
- // locally on working node\r
- //\r
-\r
- // must be deleted to store trees\r
- if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;\r
-\r
- // open temporary file and copy trees to the ouptut container\r
-\r
- TChain* chain = 0;\r
- TTree* tree1 = 0;\r
- TTree* tree2 = 0;\r
- TTree* tree3 = 0;\r
- TTree* tree4 = 0;\r
- TTree* tree5 = 0;\r
- //\r
- chain = new TChain("dNdPtTree");\r
- if(chain) { \r
- chain->Add("jotwinow_Temp_Trees.root");\r
- tree1 = chain->CopyTree("1");\r
- delete chain; chain=0; \r
- }\r
- if(tree1) tree1->Print();\r
-\r
- //\r
- chain = new TChain("V0s");\r
- if(chain) { \r
- chain->Add("jotwinow_Temp_Trees.root");\r
- tree2 = chain->CopyTree("1");\r
- delete chain; chain=0; \r
- }\r
- if(tree2) tree2->Print();\r
-\r
- //\r
- chain = new TChain("dEdx");\r
- if(chain) { \r
- chain->Add("jotwinow_Temp_Trees.root");\r
- tree3 = chain->CopyTree("1");\r
- delete chain; chain=0; \r
- }\r
- if(tree3) tree3->Print();\r
-\r
- //\r
- chain = new TChain("Laser");\r
- if(chain) { \r
- chain->Add("jotwinow_Temp_Trees.root");\r
- tree4 = chain->CopyTree("1");\r
- delete chain; chain=0; \r
- }\r
- if(tree4) tree4->Print();\r
-\r
- //\r
- chain = new TChain("MCEffTree");\r
- if(chain) { \r
- chain->Add("jotwinow_Temp_Trees.root");\r
- tree5 = chain->CopyTree("1");\r
- delete chain; chain=0; \r
- }\r
- if(tree5) tree5->Print();\r
-\r
-\r
- OpenFile(1);\r
-\r
- if(tree1) fOutput->Add(tree1);\r
- if(tree2) fOutput->Add(tree2);\r
- if(tree3) fOutput->Add(tree3);\r
- if(tree4) fOutput->Add(tree4);\r
- if(tree5) fOutput->Add(tree5);\r
- \r
- // Post output data.\r
- PostData(1, fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
-{\r
- // Called one at the end \r
- /*\r
- fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));\r
- if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
- TChain* chain = new TChain("dNdPtTree");\r
- if(!chain) return;\r
- chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");\r
- TTree *tree = chain->CopyTree("1");\r
- if (chain) { delete chain; chain=0; }\r
- if(!tree) return;\r
- tree->Print();\r
- fOutputSummary = tree;\r
-\r
- if (!fOutputSummary) {\r
- Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
- return;\r
- }\r
- */\r
-\r
- PostData(1, fOutput);\r
-\r
-}\r
+/**************************************************************************
+* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* *
+* Author: The ALICE Off-line Project. *
+* Contributors are mentioned in the code where appropriate. *
+* *
+* Permission to use, copy, modify and distribute this software and its *
+* documentation strictly for non-commercial purposes is hereby granted *
+* without fee, provided that the above copyright notice appears in all *
+* copies and that both the copyright notice and this permission notice *
+* appear in the supporting documentation. The authors make no claims *
+* about the suitability of this software for any purpose. It is *
+* provided "as is" without express or implied warranty. *
+**************************************************************************/
+
+#include "iostream"
+
+#include <TPDGCode.h>
+#include <TDatabasePDG.h>
+
+#include "TChain.h"
+#include "TTreeStream.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include "TList.h"
+#include "TObjArray.h"
+#include "TFile.h"
+#include "TMatrixD.h"
+#include "TRandom3.h"
+
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliInputEventHandler.h"
+#include "AliStack.h"
+#include "AliTrackReference.h"
+
+#include "AliPhysicsSelection.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliESDfriend.h"
+#include "AliMCEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDVertex.h"
+#include "AliTracker.h"
+#include "AliGeomManager.h"
+
+#include "AliCentrality.h"
+#include "AliESDVZERO.h"
+#include "AliMultiplicity.h"
+
+#include "AliESDtrackCuts.h"
+#include "AliMCEventHandler.h"
+#include "AlidNdPt.h"
+#include "AlidNdPtEventCuts.h"
+#include "AlidNdPtAcceptanceCuts.h"
+
+#include "AlidNdPtTrackDumpTask.h"
+#include "AliKFParticle.h"
+#include "AliESDv0.h"
+
+using namespace std;
+
+ClassImp(AlidNdPtTrackDumpTask)
+
+//_____________________________________________________________________________
+AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name)
+ : AliAnalysisTaskSE(name)
+ , fESD(0)
+ , fMC(0)
+ , fESDfriend(0)
+ , fOutput(0)
+ , fPitList(0)
+ , fUseMCInfo(kFALSE)
+ , fUseESDfriends(kFALSE)
+ , fdNdPtEventCuts(0)
+ , fdNdPtAcceptanceCuts(0)
+ , fdNdPtRecAcceptanceCuts(0)
+ , fEsdTrackCuts(0)
+ , fTrigger(AliTriggerAnalysis::kMB1)
+ , fAnalysisMode(AlidNdPtHelper::kTPC)
+ , fTreeSRedirector(0)
+ , fCentralityEstimator(0)
+ , fLowPtTrackDownscaligF(0)
+ , fLowPtV0DownscaligF(0)
+ , fProcessAll(kFALSE)
+ , fProcessCosmics(kFALSE)
+ , fTree1(0)
+ , fTree2(0)
+ , fTree3(0)
+ , fTree4(0)
+ , fTree5(0)
+ , fTree6(0)
+{
+ // Constructor
+
+ // Define input and output slots here
+ DefineOutput(1, TTree::Class());
+ DefineOutput(2, TTree::Class());
+ DefineOutput(3, TTree::Class());
+ DefineOutput(4, TTree::Class());
+ DefineOutput(5, TTree::Class());
+ DefineOutput(6, TTree::Class());
+}
+
+//_____________________________________________________________________________
+AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()
+{
+ if(fOutput) delete fOutput; fOutput =0;
+ if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0;
+
+ if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL;
+ if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;
+ if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL;
+ if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;
+
+
+}
+
+//____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::Notify()
+{
+ static Int_t count = 0;
+ count++;
+ TTree *chain = (TChain*)GetInputData(0);
+ if(chain)
+ {
+ Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
+ }
+
+return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::UserCreateOutputObjects()
+{
+ // Create histograms
+ // Called once
+ fOutput = new TList;
+ fOutput->SetOwner();
+
+ //
+ // create temporary file for output tree
+ fTreeSRedirector = new TTreeSRedirector("jotwinow_Temp_Trees.root");
+
+ fTree1 = new TTree;
+ fTree2 = new TTree;
+ fTree3 = new TTree;
+ fTree4 = new TTree;
+ fTree5 = new TTree;
+ fTree6 = new TTree;
+
+ fOutput->Add(fTree1);
+ fOutput->Add(fTree2);
+ fOutput->Add(fTree3);
+ fOutput->Add(fTree4);
+ fOutput->Add(fTree5);
+ fOutput->Add(fTree6);
+
+ PostData(1, fTree1);
+ PostData(2, fTree2);
+ PostData(3, fTree3);
+ PostData(4, fTree4);
+ PostData(5, fTree5);
+ PostData(6, fTree6);
+
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::UserExec(Option_t *)
+{
+ //
+ // Called for each event
+ //
+
+ // ESD event
+ fESD = (AliESDEvent*) (InputEvent());
+ if (!fESD) {
+ Printf("ERROR: ESD event not available");
+ return;
+ }
+
+ // MC event
+ if(fUseMCInfo) {
+ fMC = MCEvent();
+ if (!fMC) {
+ Printf("ERROR: MC event not available");
+ return;
+ }
+ }
+
+ if(fUseESDfriends) {
+ fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
+ if(!fESDfriend) {
+ Printf("ERROR: ESD friends not available");
+ }
+ }
+
+ //
+ if(fProcessAll) {
+ ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC
+ }
+ else {
+ Process(fESD,fMC,fESDfriend); // only global and TPC tracks
+ }
+
+ //
+ ProcessV0(fESD,fMC,fESDfriend);
+ ProcessLaser(fESD,fMC,fESDfriend);
+ ProcessdEdx(fESD,fMC,fESDfriend);
+
+ if (fProcessCosmics) { ProcessCosmics(fESD); }
+ if(IsUseMCInfo()) { ProcessMCEff(fESD,fMC,fESDfriend); }
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessCosmics(AliESDEvent *const event)
+{
+ //
+ // Select real events with high-pT tracks
+ //
+ if(!event) {
+ AliDebug(AliLog::kError, "event not available");
+ return;
+ }
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+
+ // check for cosmic pairs
+ //
+ // find cosmic pairs trigger by random trigger
+ //
+ //
+ AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
+ AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC();
+ const Double_t kMinPt=0.8;
+ const Double_t kMinPtMax=0.8;
+ const Double_t kMinNcl=50;
+ const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
+ Int_t ntracks=event->GetNumberOfTracks();
+ // Float_t dcaTPC[2]={0,0};
+ // Float_t covTPC[3]={0,0,0};
+
+ UInt_t specie = event->GetEventSpecie(); // skip laser events
+ if (specie==AliRecoParam::kCalib) return;
+
+
+
+ for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
+ AliESDtrack *track0 = event->GetTrack(itrack0);
+ if (!track0) continue;
+ if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;
+
+ if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;
+ if (track0->Pt() < kMinPt) continue;
+ if (track0->GetTPCncls() < kMinNcl) continue;
+ if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue;
+ if (track0->GetKinkIndex(0)>0) continue;
+ const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
+ //rm primaries
+ //
+ //track0->GetImpactParametersTPC(dcaTPC,covTPC);
+ //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
+ //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
+ // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
+ for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
+ AliESDtrack *track1 = event->GetTrack(itrack1);
+ if (!track1) continue;
+ if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;
+ if (track1->GetKinkIndex(0)>0) continue;
+ if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;
+ if (track1->Pt() < kMinPt) continue;
+ if (track1->GetTPCncls()<kMinNcl) continue;
+ if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
+ if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
+ //track1->GetImpactParametersTPC(dcaTPC,covTPC);
+ // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
+ //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
+ //
+ const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
+ //
+ Bool_t isPair=kTRUE;
+ for (Int_t ipar=0; ipar<5; ipar++){
+ if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
+ if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
+ }
+ if (!isPair) continue;
+ if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
+ //delta with correct sign
+ /*
+ TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
+ TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
+ TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
+ */
+ if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign
+ if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
+ if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
+ if (!isPair) continue;
+ TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
+ Int_t eventNumber = event->GetEventNumberInFile();
+ //Bool_t hasFriend = kFALSE;
+ //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
+ //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
+ // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();
+ //
+ //
+ Int_t ntracksSPD = vertexSPD->GetNContributors();
+ Int_t ntracksTPC = vertexTPC->GetNContributors();
+ Int_t runNumber = event->GetRunNumber();
+ Int_t timeStamp = event->GetTimeStamp();
+ ULong64_t triggerMask = event->GetTriggerMask();
+ Float_t magField = event->GetMagneticField();
+ TObjString triggerClass = event->GetFiredTriggerClasses().Data();
+
+ //
+ // Dump to the tree
+ // vertex
+ // TPC-ITS tracks
+ //
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"CosmicPairs"<<
+ "fileName.="<<&fileName<< // file name
+ "runNumber="<<runNumber<< // run number
+ "evtTimeStamp="<<timeStamp<< // time stamp of event
+ "evtNumberInFile="<<eventNumber<< // event number
+ "trigger="<<triggerMask<< // trigger
+ "triggerClass="<<&triggerClass<< // trigger
+ "Bz="<<magField<< // magnetic field
+ //
+ "multSPD="<<ntracksSPD<<
+ "multTPC="<<ntracksTPC<<
+ "vertSPD.="<<vertexSPD<< //primary vertex -SPD
+ "vertTPC.="<<vertexTPC<< //primary vertex -TPC
+ "t0.="<<track0<< //track0
+ "t1.="<<track1<< //track1
+ "\n";
+ }
+ }
+}
+
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Select real events with high-pT tracks
+ //
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AlidNdPtEventCuts *evtCuts = GetEventCuts();
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ Printf("ERROR cuts not available");
+ return;
+ }
+
+ // trigger selection
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // trigger
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) return;
+ //SetPhysicsTriggerSelection(physicsSelection);
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) return;
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
+ }
+ }
+
+ // centrality determination
+ Float_t centralityF = -1;
+ AliCentrality *esdCentrality = esdEvent->GetCentrality();
+ centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // use MC information
+ AliHeader* header = 0;
+ AliGenEventHeader* genHeader = 0;
+ AliStack* stack = 0;
+ TArrayF vtxMC(3);
+
+ Int_t multMCTrueTracks = 0;
+ if(IsUseMCInfo())
+ {
+ //
+ if(!mcEvent) {
+ AliDebug(AliLog::kError, "mcEvent not available");
+ return;
+ }
+ // get MC event header
+ header = mcEvent->Header();
+ if (!header) {
+ AliDebug(AliLog::kError, "Header not available");
+ return;
+ }
+ // MC particle stack
+ stack = mcEvent->Stack();
+ if (!stack) {
+ AliDebug(AliLog::kError, "Stack not available");
+ return;
+ }
+
+ // get MC vertex
+ genHeader = header->GenEventHeader();
+ if (!genHeader) {
+ AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+ return;
+ }
+ genHeader->PrimaryVertex(vtxMC);
+
+ // multipliticy of all MC primary tracks
+ // in Zv, pt and eta ranges)
+ multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+
+ } // end bUseMC
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+
+ // check event cuts
+ if(isEventOK && isEventTriggered)
+ {
+
+ //
+ // get IR information
+ //
+ AliESDHeader *esdHeader = 0;
+ esdHeader = esdEvent->GetHeader();
+ if(!esdHeader) return;
+ //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
+ //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
+
+ // Use when Peter commit the changes in the header
+ Int_t ir1 = 0;
+ Int_t ir2 = 0;
+
+ //
+ Double_t vert[3] = {0};
+ vert[0] = vtxESD->GetXv();
+ vert[1] = vtxESD->GetYv();
+ vert[2] = vtxESD->GetZv();
+ Int_t mult = vtxESD->GetNContributors();
+ Double_t bz = esdEvent->GetMagneticField();
+ Double_t runNumber = esdEvent->GetRunNumber();
+ Double_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ // high pT tracks
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->Charge()==0) continue;
+ if(!esdTrackCuts->AcceptTrack(track)) continue;
+ if(!accCuts->AcceptTrack(track)) continue;
+
+ // downscale low-pT tracks
+ Double_t scalempt= TMath::Min(track->Pt(),10.);
+ Double_t downscaleF = gRandom->Rndm();
+ downscaleF *= fLowPtTrackDownscaligF;
+ if(TMath::Exp(2*scalempt)<downscaleF) continue;
+ //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
+
+ AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
+ if (!tpcInner) continue;
+ // transform to the track reference frame
+ Bool_t isOK = kFALSE;
+ isOK = tpcInner->Rotate(track->GetAlpha());
+ isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
+ if(!isOK) continue;
+
+ // Dump to the tree
+ // vertex
+ // TPC-ITS tracks
+ //
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"highPt"<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "Bz="<<bz<<
+ "vtxESD.="<<vtxESD<<
+ "IRtot="<<ir1<<
+ "IRint2="<<ir2<<
+ "mult="<<mult<<
+ "esdTrack.="<<track<<
+ "centralityF="<<centralityF<<
+ "\n";
+ }
+ }
+
+}
+
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Process laser events
+ //
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // laser events
+ const AliESDHeader* esdHeader = esdEvent->GetHeader();
+ if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib)
+ {
+ Int_t countLaserTracks = 0;
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+
+ if(track->GetTPCInnerParam()) countLaserTracks++;
+ }
+
+ if(countLaserTracks > 100)
+ {
+ Double_t runNumber = esdEvent->GetRunNumber();
+ Double_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+ Double_t bz = esdEvent->GetMagneticField();
+
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"Laser"<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "Bz="<<bz<<
+ "multTPCtracks="<<countLaserTracks<<
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessAll(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
+{
+ //
+ // Select real events with high-pT tracks
+ // Calculate and stor in the output tree:
+ // TPC constrained tracks
+ // InnerParams constrained tracks
+ // TPC-ITS tracks
+ // ITSout-InnerParams tracks
+ // chi2 distance between TPC constrained and TPC-ITS tracks
+ // chi2 distance between TPC refitted constrained and TPC-ITS tracks
+ // chi2 distance between ITSout and InnerParams tracks
+ // MC information:
+ // track references at ITSin, TPCin; InnerParam at first TPC track reference,
+ // particle ID, mother ID, production mechanism ...
+ //
+
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AlidNdPtEventCuts *evtCuts = GetEventCuts();
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ AliDebug(AliLog::kError, "cuts not available");
+ return;
+ }
+
+ // trigger selection
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // trigger
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) return;
+ //SetPhysicsTriggerSelection(physicsSelection);
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) return;
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
+ }
+ }
+
+ // centrality determination
+ Float_t centralityF = -1;
+ AliCentrality *esdCentrality = esdEvent->GetCentrality();
+ centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // use MC information
+ AliHeader* header = 0;
+ AliGenEventHeader* genHeader = 0;
+ AliStack* stack = 0;
+ TArrayF vtxMC(3);
+
+ Int_t multMCTrueTracks = 0;
+ if(IsUseMCInfo())
+ {
+ //
+ if(!mcEvent) {
+ AliDebug(AliLog::kError, "mcEvent not available");
+ return;
+ }
+ // get MC event header
+ header = mcEvent->Header();
+ if (!header) {
+ AliDebug(AliLog::kError, "Header not available");
+ return;
+ }
+ // MC particle stack
+ stack = mcEvent->Stack();
+ if (!stack) {
+ AliDebug(AliLog::kError, "Stack not available");
+ return;
+ }
+
+ // get MC vertex
+ genHeader = header->GenEventHeader();
+ if (!genHeader) {
+ AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+ return;
+ }
+ genHeader->PrimaryVertex(vtxMC);
+
+ // multipliticy of all MC primary tracks
+ // in Zv, pt and eta ranges)
+ multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+
+ } // end bUseMC
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+
+ // check event cuts
+ if(isEventOK && isEventTriggered)
+ {
+ //
+ // get IR information
+ //
+ AliESDHeader *esdHeader = 0;
+ esdHeader = esdEvent->GetHeader();
+ if(!esdHeader) return;
+ //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s
+ //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval
+ //
+ Int_t ir1 = 0;
+ Int_t ir2 = 0;
+
+ //
+ Double_t vert[3] = {0};
+ vert[0] = vtxESD->GetXv();
+ vert[1] = vtxESD->GetYv();
+ vert[2] = vtxESD->GetZv();
+ Int_t mult = vtxESD->GetNContributors();
+ Double_t bz = esdEvent->GetMagneticField();
+ Double_t runNumber = esdEvent->GetRunNumber();
+ Double_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ // high pT tracks
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->Charge()==0) continue;
+ if(!esdTrackCuts->AcceptTrack(track)) continue;
+ if(!accCuts->AcceptTrack(track)) continue;
+
+ // downscale low-pT tracks
+ Double_t scalempt= TMath::Min(track->Pt(),10.);
+ Double_t downscaleF = gRandom->Rndm();
+ downscaleF *= fLowPtTrackDownscaligF;
+ if(TMath::Exp(2*scalempt)<downscaleF) continue;
+ //printf("TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
+
+ // Dump to the tree
+ // vertex
+ // TPC constrained tracks
+ // InnerParams constrained tracks
+ // TPC-ITS tracks
+ // ITSout-InnerParams tracks
+ // chi2 distance between TPC constrained and TPC-ITS tracks
+ // chi2 distance between TPC refitted constrained and TPC-ITS tracks
+ // chi2 distance between ITSout and InnerParams tracks
+ // MC information
+
+ Double_t x[3]; track->GetXYZ(x);
+ Double_t b[3]; AliTracker::GetBxByBz(x,b);
+
+ //
+ // Transform TPC inner params to track reference frame
+ //
+ Bool_t isOKtpcInner = kFALSE;
+ AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());
+ if (tpcInner) {
+ // transform to the track reference frame
+ isOKtpcInner = tpcInner->Rotate(track->GetAlpha());
+ isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
+ }
+
+ //
+ // Constrain TPC inner to vertex
+ // clone TPCinner has to be deleted
+ //
+ Bool_t isOKtpcInnerC = kFALSE;
+ AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam()));
+ if (tpcInnerC) {
+ isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b);
+ isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha());
+ isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
+ }
+
+ //
+ // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex
+ // Clone track InnerParams has to be deleted
+ //
+ Bool_t isOKtrackInnerC = kFALSE;
+ AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam()));
+ if (trackInnerC) {
+ isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);
+ isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha());
+ isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());
+ }
+
+ //
+ // calculate chi2 between vi and vj vectors
+ // with covi and covj covariance matrices
+ // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
+ //
+ TMatrixD deltaT(5,1), deltaTtrackC(5,1);
+ TMatrixD delta(1,5), deltatrackC(1,5);
+ TMatrixD covarM(5,5), covarMtrackC(5,5);
+ TMatrixD chi2(1,1);
+ TMatrixD chi2trackC(1,1);
+
+ if(isOKtpcInnerC && isOKtrackInnerC)
+ {
+ for (Int_t ipar=0; ipar<5; ipar++) {
+ deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
+ delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
+
+ deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
+ deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar];
+
+ for (Int_t jpar=0; jpar<5; jpar++) {
+ Int_t index=track->GetIndex(ipar,jpar);
+ covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
+ covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index];
+ }
+ }
+
+ // chi2 distance TPC constrained and TPC+ITS
+ TMatrixD covarMInv = covarM.Invert();
+ TMatrixD mat2 = covarMInv*deltaT;
+ chi2 = delta*mat2;
+ //chi2.Print();
+
+ // chi2 distance TPC refitted constrained and TPC+ITS
+ TMatrixD covarMInvtrackC = covarMtrackC.Invert();
+ TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC;
+ chi2trackC = deltatrackC*mat2trackC;
+ //chi2trackC.Print();
+ }
+
+
+ //
+ // Propagate ITSout to TPC inner wall
+ // and calculate chi2 distance to track (InnerParams)
+ //
+ const Double_t kTPCRadius=85;
+ const Double_t kStep=3;
+
+ // clone track InnerParams has to be deleted
+ Bool_t isOKtrackInnerC2 = kFALSE;
+ AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam()));
+ if (trackInnerC2) {
+ isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE);
+ }
+
+ Bool_t isOKouterITSc = kFALSE;
+ AliExternalTrackParam *outerITSc = NULL;
+ TMatrixD chi2OuterITS(1,1);
+
+ if(esdFriend && esdFriend->TestSkipBit()==kFALSE)
+ {
+ // propagate ITSout to TPC inner wall
+ AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack);
+
+ if(friendTrack)
+ {
+ outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut());
+ if(outerITSc)
+ {
+ isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE);
+ isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha());
+ isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());
+
+ //
+ // calculate chi2 between outerITS and innerParams
+ // cov without z-coordinate at the moment
+ //
+ TMatrixD deltaTouterITS(4,1);
+ TMatrixD deltaouterITS(1,4);
+ TMatrixD covarMouterITS(4,4);
+
+ if(isOKtrackInnerC2 && isOKouterITSc) {
+ Int_t kipar = 0;
+ Int_t kjpar = 0;
+ for (Int_t ipar=0; ipar<5; ipar++) {
+ if(ipar!=1) {
+ deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
+ deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];
+ }
+
+ kjpar=0;
+ for (Int_t jpar=0; jpar<5; jpar++) {
+ Int_t index=outerITSc->GetIndex(ipar,jpar);
+ if(ipar !=1 || jpar!=1) {
+ covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];
+ }
+ if(jpar!=1) kjpar++;
+ }
+ if(ipar!=1) kipar++;
+ }
+
+ // chi2 distance ITSout and InnerParams
+ TMatrixD covarMInvouterITS = covarMouterITS.Invert();
+ TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;
+ chi2OuterITS = deltaouterITS*mat2outerITS;
+ //chi2OuterITS.Print();
+ }
+ }
+ }
+ }
+
+ //
+ // MC info
+ //
+ TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL;
+ TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL;
+ Int_t mech=-1, mechTPC=-1, mechITS=-1;
+ Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE;
+ Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE;
+ Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE;
+ Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE;
+
+ AliTrackReference *refTPCIn = NULL;
+ AliTrackReference *refITS = NULL;
+
+ Bool_t isOKtrackInnerC3 = kFALSE;
+ AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
+
+ if(IsUseMCInfo())
+ {
+ if(!stack) return;
+
+ //
+ // global track
+ //
+ Int_t label = TMath::Abs(track->GetLabel());
+ particle = stack->Particle(label);
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.)
+ {
+ particleMother = GetMother(particle,stack);
+ mech = particle->GetUniqueID();
+ isPrim = stack->IsPhysicalPrimary(label);
+ isFromStrangess = IsFromStrangeness(label,stack);
+ isFromConversion = IsFromConversion(label,stack);
+ isFromMaterial = IsFromMaterial(label,stack);
+ }
+
+ //
+ // TPC track
+ //
+ Int_t labelTPC = TMath::Abs(track->GetTPCLabel());
+ particleTPC = stack->Particle(labelTPC);
+ if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.)
+ {
+ particleMotherTPC = GetMother(particleTPC,stack);
+ mechTPC = particleTPC->GetUniqueID();
+ isPrimTPC = stack->IsPhysicalPrimary(labelTPC);
+ isFromStrangessTPC = IsFromStrangeness(labelTPC,stack);
+ isFromConversionTPC = IsFromConversion(labelTPC,stack);
+ isFromMaterialTPC = IsFromMaterial(labelTPC,stack);
+ }
+
+ //
+ // store first track reference
+ // for TPC track
+ //
+ TParticle *part=0;
+ TClonesArray *trefs=0;
+ Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs);
+
+ if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.)
+ {
+ Int_t nTrackRef = trefs->GetEntries();
+ //printf("nTrackRef %d \n",nTrackRef);
+
+ Int_t countITS = 0;
+ for (Int_t iref = 0; iref < nTrackRef; iref++)
+ {
+ AliTrackReference *ref = (AliTrackReference *)trefs->At(iref);
+
+ // Ref. in the middle ITS
+ if(ref && ref->DetectorId()==AliTrackReference::kITS)
+ {
+ if(!refITS && countITS==2) {
+ refITS = ref;
+ //printf("refITS %p \n",refITS);
+ }
+ countITS++;
+ }
+
+ // TPC
+ if(ref && ref->DetectorId()==AliTrackReference::kTPC)
+ {
+ if(!refTPCIn) {
+ refTPCIn = ref;
+ //printf("refTPCIn %p \n",refTPCIn);
+ //break;
+ }
+ }
+ }
+
+ // transform inner params to TrackRef
+ // reference frame
+ if(refTPCIn && trackInnerC3)
+ {
+ Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X());
+ isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi);
+ isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);
+ }
+ }
+
+ //
+ // ITS track
+ //
+ Int_t labelITS = TMath::Abs(track->GetITSLabel());
+ particleITS = stack->Particle(labelITS);
+ if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)
+ {
+ particleMotherITS = GetMother(particleITS,stack);
+ mechITS = particleITS->GetUniqueID();
+ isPrimITS = stack->IsPhysicalPrimary(labelITS);
+ isFromStrangessITS = IsFromStrangeness(labelITS,stack);
+ isFromConversionITS = IsFromConversion(labelITS,stack);
+ isFromMaterialITS = IsFromMaterial(labelITS,stack);
+ }
+ }
+
+ //
+ Bool_t dumpToTree=kFALSE;
+
+ if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
+ if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
+ if(fUseMCInfo && isOKtrackInnerC3) dumpToTree = kTRUE;
+
+ //
+ if(fTreeSRedirector && dumpToTree)
+ {
+ (*fTreeSRedirector)<<"highPt"<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "Bz="<<bz<<
+ "vtxESD.="<<vtxESD<<
+ "IRtot="<<ir1<<
+ "IRint2="<<ir2<<
+ "mult="<<mult<<
+ "esdTrack.="<<track<<
+ "extTPCInnerC.="<<tpcInnerC<<
+ "extInnerParamC.="<<trackInnerC<<
+ "extInnerParam.="<<trackInnerC2<<
+ "extOuterITS.="<<outerITSc<<
+ "extInnerParamRef.="<<trackInnerC3<<
+ "refTPCIn.="<<refTPCIn<<
+ "refITS.="<<refITS<<
+ "chi2TPCInnerC="<<chi2(0,0)<<
+ "chi2InnerC="<<chi2trackC(0,0)<<
+ "chi2OuterITS="<<chi2OuterITS(0,0)<<
+ "centralityF="<<centralityF<<
+ "particle.="<<particle<<
+ "particleMother.="<<particleMother<<
+ "mech="<<mech<<
+ "isPrim="<<isPrim<<
+ "isFromStrangess="<<isFromStrangess<<
+ "isFromConversion="<<isFromConversion<<
+ "isFromMaterial="<<isFromMaterial<<
+ "particleTPC.="<<particleTPC<<
+ "particleMotherTPC.="<<particleMotherTPC<<
+ "mechTPC="<<mechTPC<<
+ "isPrimTPC="<<isPrimTPC<<
+ "isFromStrangessTPC="<<isFromStrangessTPC<<
+ "isFromConversionTPC="<<isFromConversionTPC<<
+ "isFromMaterialTPC="<<isFromMaterialTPC<<
+ "particleITS.="<<particleITS<<
+ "particleMotherITS.="<<particleMotherITS<<
+ "mechITS="<<mechITS<<
+ "isPrimITS="<<isPrimITS<<
+ "isFromStrangessITS="<<isFromStrangessITS<<
+ "isFromConversionITS="<<isFromConversionITS<<
+ "isFromMaterialITS="<<isFromMaterialITS<<
+ "\n";
+ }
+
+ if(tpcInnerC) delete tpcInnerC;
+ if(trackInnerC) delete trackInnerC;
+ if(trackInnerC2) delete trackInnerC2;
+ if(outerITSc) delete outerITSc;
+ if(trackInnerC3) delete trackInnerC3;
+ }
+ }
+
+}
+
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Fill tree for efficiency studies MC only
+
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ if(!mcEvent) {
+ AliDebug(AliLog::kError, "mcEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AlidNdPtEventCuts *evtCuts = GetEventCuts();
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ AliDebug(AliLog::kError, "cuts not available");
+ return;
+ }
+
+ // trigger selection
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // trigger
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) return;
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) return;
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
+ }
+ }
+
+ // centrality determination
+ Float_t centralityF = -1;
+ AliCentrality *esdCentrality = esdEvent->GetCentrality();
+ centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+ // use MC information
+ AliHeader* header = 0;
+ AliGenEventHeader* genHeader = 0;
+ AliStack* stack = 0;
+ TArrayF vtxMC(3);
+
+ Int_t multMCTrueTracks = 0;
+ if(IsUseMCInfo())
+ {
+ //
+ if(!mcEvent) {
+ AliDebug(AliLog::kError, "mcEvent not available");
+ return;
+ }
+ // get MC event header
+ header = mcEvent->Header();
+ if (!header) {
+ AliDebug(AliLog::kError, "Header not available");
+ return;
+ }
+ // MC particle stack
+ stack = mcEvent->Stack();
+ if (!stack) {
+ AliDebug(AliLog::kError, "Stack not available");
+ return;
+ }
+
+ // get MC vertex
+ genHeader = header->GenEventHeader();
+ if (!genHeader) {
+ AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+ return;
+ }
+ genHeader->PrimaryVertex(vtxMC);
+
+ // multipliticy of all MC primary tracks
+ // in Zv, pt and eta ranges)
+ multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+
+ } // end bUseMC
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+ // check event cuts
+ if(isEventOK && isEventTriggered)
+ {
+ if(IsUseMCInfo())
+ {
+ if(!stack) return;
+
+ //
+ // MC info
+ //
+ TParticle *particle=NULL;
+ TParticle *particleMother=NULL;
+ Int_t mech=-1;
+
+ // reco event info
+ Double_t vert[3] = {0};
+ vert[0] = vtxESD->GetXv();
+ vert[1] = vtxESD->GetYv();
+ vert[2] = vtxESD->GetZv();
+ Int_t mult = vtxESD->GetNContributors();
+ Double_t bz = esdEvent->GetMagneticField();
+ Double_t runNumber = esdEvent->GetRunNumber();
+ Double_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ // loop over MC stack
+ for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
+ {
+ particle = stack->Particle(iMc);
+ if (!particle)
+ continue;
+
+ // only charged particles
+ if(!particle->GetPDG()) continue;
+ Double_t charge = particle->GetPDG()->Charge()/3.;
+ if (TMath::Abs(charge) < 0.001)
+ continue;
+
+ // only primary particles
+ Bool_t prim = stack->IsPhysicalPrimary(iMc);
+ if(!prim) continue;
+
+ // downscale low-pT particles
+ Double_t scalempt= TMath::Min(particle->Pt(),10.);
+ Double_t downscaleF = gRandom->Rndm();
+ downscaleF *= fLowPtTrackDownscaligF;
+ if(TMath::Exp(2*scalempt)<downscaleF) continue;
+
+ // is particle in acceptance
+ if(!accCuts->AcceptTrack(particle)) continue;
+
+ // check if particle reconstructed
+ Bool_t isRec = kFALSE;
+ Int_t trackIndex = -1;
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->Charge()==0) continue;
+ if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track))
+ {
+ Int_t label = TMath::Abs(track->GetLabel());
+ if(label == iMc) {
+ isRec = kTRUE;
+ trackIndex = iTrack;
+ break;
+ }
+ }
+ }
+
+ // Store information in the output tree
+ AliESDtrack *recTrack = NULL;
+ if(trackIndex>-1) {
+ recTrack = esdEvent->GetTrack(trackIndex);
+ } else {
+ recTrack = new AliESDtrack();
+ }
+
+ particleMother = GetMother(particle,stack);
+ mech = particle->GetUniqueID();
+
+ //MC particle track length
+ Double_t tpcTrackLength = 0.;
+ AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc);
+ if(mcParticle) {
+ Int_t counter;
+ tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0);
+ }
+
+
+ //
+ if(fTreeSRedirector) {
+ (*fTreeSRedirector)<<"MCEffTree"<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "Bz="<<bz<<
+ "vtxESD.="<<vtxESD<<
+ "mult="<<mult<<
+ "esdTrack.="<<recTrack<<
+ "isRec="<<isRec<<
+ "tpcTrackLength="<<tpcTrackLength<<
+ "particle.="<<particle<<
+ "particleMother.="<<particleMother<<
+ "mech="<<mech<<
+ "\n";
+ }
+
+ if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
+ }
+ }
+ }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::IsHighDeDxParticle(AliESDtrack * track) {
+ //
+ // check if particle is Z > 1
+ //
+ if (track->GetTPCNcls() < 60) return kFALSE;
+ Double_t mom = track->GetInnerParam()->GetP();
+ if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization
+ Float_t dca[2], bCov[3];
+ track->GetImpactParameters(dca,bCov);
+ //
+
+ Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
+
+ if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE;
+
+ return kFALSE;
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Select real events with V0 (K0s and Lambda) high-pT candidates
+ //
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AlidNdPtEventCuts *evtCuts = GetEventCuts();
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ AliDebug(AliLog::kError, "cuts not available");
+ return;
+ }
+
+ // trigger selection
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // trigger
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) return;
+ //SetPhysicsTriggerSelection(physicsSelection);
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) return;
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
+ }
+ }
+
+ // centrality determination
+ Float_t centralityF = -1;
+ AliCentrality *esdCentrality = esdEvent->GetCentrality();
+ centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
+
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+ // check event cuts
+ if(isEventOK && isEventTriggered) {
+ //
+ // Dump the pt downscaled V0 into the tree
+ //
+ Int_t ntracks = esdEvent->GetNumberOfTracks();
+ Int_t nV0s = esdEvent->GetNumberOfV0s();
+ Int_t run = esdEvent->GetRunNumber();
+ Int_t time= esdEvent->GetTimeStamp();
+ Int_t evNr=esdEvent->GetEventNumberInFile();
+ Double_t bz = esdEvent->GetMagneticField();
+
+
+ for (Int_t iv0=0; iv0<nV0s; iv0++){
+ AliESDv0 * v0 = esdEvent->GetV0(iv0);
+ if (!v0) continue;
+ AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0));
+ AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1));
+ if (!track0) continue;
+ if (!track1) continue;
+ if (track0->GetSign()<0) {
+ track1 = esdEvent->GetTrack(v0->GetIndex(0));
+ track0 = esdEvent->GetTrack(v0->GetIndex(1));
+ }
+ //
+ Bool_t isDownscaled = IsV0Downscaled(v0);
+ if (isDownscaled) continue;
+ AliKFParticle kfparticle; //
+ Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
+ if (type==0) continue;
+
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"V0s"<<
+ "isDownscaled="<<isDownscaled<<
+ "Bz="<<bz<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<run<<
+ "evtTimeStamp="<<time<<
+ "evtNumberInFile="<<evNr<<
+ "type="<<type<<
+ "ntracks="<<ntracks<<
+ "v0.="<<v0<<
+ "kf.="<<&kfparticle<<
+ "track0.="<<track0<<
+ "track1.="<<track1<<
+ "centralityF="<<centralityF<<
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::ProcessdEdx(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Select real events with large TPC dEdx signal
+ //
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AlidNdPtEventCuts *evtCuts = GetEventCuts();
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ AliDebug(AliLog::kError, "cuts not available");
+ return;
+ }
+
+ // get file name
+ TTree *chain = (TChain*)GetInputData(0);
+ if(!chain) {
+ Printf("ERROR: Could not receive input chain");
+ return;
+ }
+ TObjString fileName(chain->GetCurrentFile()->GetName());
+
+ // trigger
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler)
+ {
+ Printf("ERROR: Could not receive input handler");
+ return;
+ }
+
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) return;
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) return;
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
+ }
+ }
+
+ // get reconstructed vertex
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+
+ // check event cuts
+ if(isEventOK && isEventTriggered)
+ {
+ Double_t vert[3] = {0};
+ vert[0] = vtxESD->GetXv();
+ vert[1] = vtxESD->GetYv();
+ vert[2] = vtxESD->GetZv();
+ Int_t mult = vtxESD->GetNContributors();
+ Double_t bz = esdEvent->GetMagneticField();
+ Double_t runNumber = esdEvent->GetRunNumber();
+ Double_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ // large dEdx
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->Charge()==0) continue;
+ if(!esdTrackCuts->AcceptTrack(track)) continue;
+ if(!accCuts->AcceptTrack(track)) continue;
+
+ if(!IsHighDeDxParticle(track)) continue;
+
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"dEdx"<<
+ "fileName.="<<&fileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "Bz="<<bz<<
+ "vtxESD.="<<vtxESD<<
+ "mult="<<mult<<
+ "esdTrack.="<<track<<
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AlidNdPtTrackDumpTask::GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle)
+{
+ //
+ // Create KF particle in case the V0 fullfill selection criteria
+ //
+ // Selection criteria
+ // 0. algorithm cut
+ // 1. track cut
+ // 3. chi2 cut
+ // 4. rough mass cut
+ // 5. Normalized pointing angle cut
+ //
+ const Double_t cutMass=0.2;
+ const Double_t kSigmaDCACut=3;
+ //
+ // 0.) algo cut - accept only on the fly
+ //
+ if (v0->GetOnFlyStatus() ==kFALSE) return 0;
+ //
+ // 1.) track cut
+ //
+ AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0));
+ AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1));
+ /*
+ TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2";
+ TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1";
+ TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100";
+ */
+ if (TMath::Abs(track0->GetTgl())>1) return 0;
+ if (TMath::Abs(track1->GetTgl())>1) return 0;
+ if ((track0->GetTPCClusterInfo(2,1))<100) return 0;
+ if ((track1->GetTPCClusterInfo(2,1))<100) return 0;
+ //if ((track0->GetITSclusters(0))<2) return 0;
+ //if ((track1->GetITSclusters(0))<2) return 0;
+ Float_t pos0[2]={0}, cov0[3]={0};
+ Float_t pos1[2]={0}, cov1[3]={0};
+ track0->GetImpactParameters(pos0,cov0);
+ track0->GetImpactParameters(pos1,cov1);
+ //
+ if (TMath::Abs(pos0[0])<kSigmaDCACut*TMath::Sqrt(cov0[0])) return 0;
+ if (TMath::Abs(pos1[0])<kSigmaDCACut*TMath::Sqrt(cov1[0])) return 0;
+ //
+ //
+ // 3.) Chi2 cut
+ //
+ Double_t chi2KF = v0->GetKFInfo(2,2,2);
+ if (chi2KF>25) return 0;
+ //
+ // 4.) Rough mass cut - 0.200 GeV
+ //
+ static Double_t masses[2]={-1};
+ if (masses[0]<0){
+ masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass();
+ masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass();
+ }
+ Double_t mass00= v0->GetEffMass(0,0);
+ Double_t mass22= v0->GetEffMass(2,2);
+ Double_t mass42= v0->GetEffMass(4,2);
+ Double_t mass24= v0->GetEffMass(2,4);
+ Bool_t massOK=kFALSE;
+ Int_t type=0;
+ Int_t ptype=0;
+ Double_t dmass=1;
+ Int_t p1=0, p2=0;
+ if (TMath::Abs(mass00-0)<cutMass) {
+ massOK=kTRUE; type+=1;
+ if (TMath::Abs(mass00-0)<dmass) {
+ ptype=1;
+ dmass=TMath::Abs(mass00-0);
+ p1=0; p2=0;
+ }
+ }
+ if (TMath::Abs(mass24-masses[1])<cutMass) {
+ massOK=kTRUE; type+=2;
+ if (TMath::Abs(mass24-masses[1])<dmass){
+ dmass = TMath::Abs(mass24-masses[1]);
+ ptype=2;
+ p1=2; p2=4;
+ }
+ }
+ if (TMath::Abs(mass42-masses[1])<cutMass) {
+ massOK=kTRUE; type+=4;
+ if (TMath::Abs(mass42-masses[1])<dmass){
+ dmass = TMath::Abs(mass42-masses[1]);
+ ptype=4;
+ p1=4; p2=2;
+ }
+ }
+ if (TMath::Abs(mass22-masses[0])<cutMass) {
+ massOK=kTRUE; type+=8;
+ if (TMath::Abs(mass22-masses[0])<dmass){
+ dmass = TMath::Abs(mass22-masses[0]);
+ ptype=8;
+ p1=2; p2=2;
+ }
+ }
+ if (type==0) return 0;
+ //
+ const Int_t spdg[5]={kPositron,kMuonPlus,kPiPlus, kKPlus, kProton};
+ const AliExternalTrackParam *paramP = v0->GetParamP();
+ const AliExternalTrackParam *paramN = v0->GetParamN();
+ if (paramP->GetSign()<0){
+ paramP=v0->GetParamP();
+ paramN=v0->GetParamN();
+ }
+ //Double_t *pparam1 = (Double_t*)paramP->GetParameter();
+ //Double_t *pparam2 = (Double_t*)paramN->GetParameter();
+ //
+ AliKFParticle kfp1( *paramP, spdg[p1] );
+ AliKFParticle kfp2( *paramN, -1 *spdg[p2] );
+ AliKFParticle V0KF;
+ (V0KF)+=kfp1;
+ (V0KF)+=kfp2;
+ kfparticle=V0KF;
+ //
+ // Pointing angle
+ //
+ Double_t errPhi = V0KF.GetErrPhi();
+ Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle());
+ if (pointAngle/errPhi>10) return 0;
+ //
+ return ptype;
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::IsV0Downscaled(AliESDv0 *const v0)
+{
+ //
+ // Downscale randomly low pt V0
+ //
+ //return kFALSE;
+ Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt());
+ Double_t scalempt= TMath::Min(maxPt,10.);
+ Double_t downscaleF = gRandom->Rndm();
+ downscaleF *= fLowPtV0DownscaligF;
+
+ //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF);
+ if (TMath::Exp(2*scalempt)<downscaleF) return kTRUE;
+ return kFALSE;
+
+ /*
+ TH1F his1("his1","his1",100,0,10);
+ TH1F his2("his2","his2",100,0,10);
+ {for (Int_t i=0; i<10000; i++){
+ Double_t rnd=gRandom->Exp(1);
+ Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm();
+ his1->Fill(rnd);
+ if (!isDownscaled) his2->Fill(rnd);
+ }}
+
+ */
+
+}
+
+
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])
+{
+ // Constrain TPC inner params constrained
+ //
+ if(!tpcInnerC) return kFALSE;
+ if(!vtx) return kFALSE;
+
+ Double_t dz[2],cov[3];
+ //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
+ //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
+ //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
+ if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
+
+
+ Double_t covar[6]; vtx->GetCovMatrix(covar);
+ Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};
+ Double_t c[3]={covar[2],0.,covar[5]};
+ Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);
+ if (chi2C>kVeryBig) return kFALSE;
+
+ if(!tpcInnerC->Update(p,c)) return kFALSE;
+
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])
+{
+ // Constrain TPC inner params constrained
+ //
+ if(!trackInnerC) return kFALSE;
+ if(!vtx) return kFALSE;
+
+ const Double_t kRadius = 2.8;
+ const Double_t kMaxStep = 1.0;
+
+ Double_t dz[2],cov[3];
+
+ //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();
+ //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE;
+ //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE;
+
+ if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;
+ if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE;
+
+ //
+ Double_t covar[6]; vtx->GetCovMatrix(covar);
+ Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};
+ Double_t c[3]={covar[2],0.,covar[5]};
+ Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);
+ if (chi2C>kVeryBig) return kFALSE;
+
+ if(!trackInnerC->Update(p,c)) return kFALSE;
+
+ return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack)
+{
+ if(!particle) return NULL;
+ if(!stack) return NULL;
+
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0));
+ TParticle* mother = NULL;
+ mother = stack->Particle(motherLabel);
+
+return mother;
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack)
+{
+ Bool_t isFromConversion = kFALSE;
+
+ if(stack) {
+ TParticle* particle = stack->Particle(label);
+
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
+ {
+ Int_t mech = particle->GetUniqueID(); // production mechanism
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);
+
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0));
+ TParticle* mother = stack->Particle(motherLabel);
+ if(mother) {
+ Int_t motherPdg = mother->GetPdgCode();
+
+ if(!isPrim && mech==5 && motherPdg==kGamma) {
+ isFromConversion=kTRUE;
+ }
+ }
+ }
+ }
+
+ return isFromConversion;
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack)
+{
+ Bool_t isFromMaterial = kFALSE;
+
+ if(stack) {
+ TParticle* particle = stack->Particle(label);
+
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
+ {
+ Int_t mech = particle->GetUniqueID(); // production mechanism
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);
+
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0));
+ TParticle* mother = stack->Particle(motherLabel);
+ if(mother) {
+ if(!isPrim && mech==13) {
+ isFromMaterial=kTRUE;
+ }
+ }
+ }
+ }
+
+ return isFromMaterial;
+}
+
+//_____________________________________________________________________________
+Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack)
+{
+ Bool_t isFromStrangeness = kFALSE;
+
+ if(stack) {
+ TParticle* particle = stack->Particle(label);
+
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0)
+ {
+ Int_t mech = particle->GetUniqueID(); // production mechanism
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);
+
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0));
+ TParticle* mother = stack->Particle(motherLabel);
+ if(mother) {
+ Int_t motherPdg = mother->GetPdgCode();
+
+ // K+-, lambda, antilambda, K0s decays
+ if(!isPrim && mech==4 &&
+ (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))
+ {
+ isFromStrangeness = kTRUE;
+ }
+ }
+ }
+ }
+
+ return isFromStrangeness;
+}
+
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::FinishTaskOutput()
+{
+ //
+ // Called one at the end
+ // locally on working node
+ //
+
+ // must be deleted to store trees
+ if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;
+
+ // open temporary file and copy trees to the ouptut container
+
+ TChain* chain = 0;
+ /*
+ TTree* tree1 = 0;
+ TTree* tree2 = 0;
+ TTree* tree3 = 0;
+ TTree* tree4 = 0;
+ TTree* tree5 = 0;
+ TTree* tree6 = 0;
+ */
+ //
+ chain = new TChain("highPt");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree1 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree1) fTree1->Print();
+
+ //
+ chain = new TChain("V0s");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree2 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree2) fTree2->Print();
+
+ //
+ chain = new TChain("dEdx");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree3 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree3) fTree3->Print();
+
+ //
+ chain = new TChain("Laser");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree4 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree4) fTree4->Print();
+
+ //
+ chain = new TChain("MCEffTree");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree5 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree5) fTree5->Print();
+
+ //
+ chain = new TChain("CosmicPairs");
+ if(chain) {
+ chain->Add("jotwinow_Temp_Trees.root");
+ fTree6 = chain->CopyTree("1");
+ delete chain; chain=0;
+ }
+ if(fTree6) fTree6->Print();
+
+
+ OpenFile(1);
+
+ // Post output data.
+ PostData(1, fTree1);
+ PostData(2, fTree2);
+ PostData(3, fTree3);
+ PostData(4, fTree4);
+ PostData(5, fTree5);
+ PostData(6, fTree6);
+}
+
+//_____________________________________________________________________________
+void AlidNdPtTrackDumpTask::Terminate(Option_t *)
+{
+ // Called one at the end
+ /*
+ fOutputSummary = dynamic_cast<TTree*> (GetOutputData(1));
+ if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;
+ TChain* chain = new TChain("highPt");
+ if(!chain) return;
+ chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root");
+ TTree *tree = chain->CopyTree("1");
+ if (chain) { delete chain; chain=0; }
+ if(!tree) return;
+ tree->Print();
+ fOutputSummary = tree;
+
+ if (!fOutputSummary) {
+ Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable %p \n", GetOutputData(1));
+ return;
+ }
+ */
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff