-/**************************************************************************\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 "TH3.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 "AliVTrack.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 "AliFilteredTreeEventCuts.h"\r
-#include "AliFilteredTreeAcceptanceCuts.h"\r
-\r
-#include "AliAnalysisTaskFilteredTree.h"\r
-#include "AliKFParticle.h"\r
-#include "AliESDv0.h"\r
-\r
-using namespace std;\r
-\r
-ClassImp(AliAnalysisTaskFilteredTree)\r
-\r
-//_____________________________________________________________________________\r
-AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(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
- , fReducePileUp(kTRUE)\r
- , fFillTree(kTRUE)\r
- , fFilteredTreeEventCuts(0)\r
- , fFilteredTreeAcceptanceCuts(0)\r
- , fFilteredTreeRecAcceptanceCuts(0)\r
- , fEsdTrackCuts(0)\r
- , fTrigger(AliTriggerAnalysis::kMB1) \r
- , fAnalysisMode(kTPCAnalysisMode) \r
- , fTreeSRedirector(0)\r
- , fCentralityEstimator(0)\r
- , fLowPtTrackDownscaligF(0)\r
- , fLowPtV0DownscaligF(0)\r
- , fProcessAll(kFALSE)\r
- , fProcessCosmics(kFALSE)\r
- , fHighPtTree(0)\r
- , fV0Tree(0)\r
- , fdEdxTree(0)\r
- , fLaserTree(0)\r
- , fMCEffTree(0)\r
- , fCosmicPairsTree(0)\r
- , fPtResPhiPtTPC(0)\r
- , fPtResPhiPtTPCc(0)\r
- , fPtResPhiPtTPCITS(0)\r
- , fPtResEtaPtTPC(0)\r
- , fPtResEtaPtTPCc(0)\r
- , fPtResEtaPtTPCITS(0)\r
- , fPtResCentPtTPC(0)\r
- , fPtResCentPtTPCc(0)\r
- , fPtResCentPtTPCITS(0)\r
-{\r
- // Constructor\r
-\r
- // Define input and output slots here\r
- DefineOutput(1, TTree::Class());\r
- DefineOutput(2, TTree::Class());\r
- DefineOutput(3, TTree::Class());\r
- DefineOutput(4, TTree::Class());\r
- DefineOutput(5, TTree::Class());\r
- DefineOutput(6, TTree::Class());\r
- DefineOutput(7, TList::Class());\r
-}\r
-\r
-//_____________________________________________________________________________\r
-AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()\r
-{\r
- //the output trees not to be deleted in case of proof analysis\r
- //Bool_t deleteTrees=kTRUE;\r
- //if ((AliAnalysisManager::GetAnalysisManager()))\r
- //{\r
- // if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == \r
- // AliAnalysisManager::kProofAnalysis)\r
- // deleteTrees=kFALSE;\r
- //}\r
- //if (deleteTrees) delete fTreeSRedirector;\r
-\r
- delete fFilteredTreeEventCuts;\r
- delete fFilteredTreeAcceptanceCuts;\r
- delete fFilteredTreeRecAcceptanceCuts;\r
- delete fEsdTrackCuts;\r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t AliAnalysisTaskFilteredTree::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 AliAnalysisTaskFilteredTree::UserCreateOutputObjects()\r
-{\r
- // Create histograms\r
- // Called once\r
-\r
- //\r
- //get the output file to make sure the trees will be associated to it\r
- OpenFile(1);\r
- fTreeSRedirector = new TTreeSRedirector();\r
- \r
- //\r
- // Create trees\r
- fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree();\r
- fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree();\r
- fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree();\r
- fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree();\r
- fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree();\r
- fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();\r
-\r
-\r
-\r
-\r
- // histogram booking\r
-\r
- Double_t minPt = 0.1; \r
- Double_t maxPt = 100.; \r
- Int_t nbinsPt = 30; \r
-\r
- Double_t logminPt = TMath::Log10(minPt);\r
- Double_t logmaxPt = TMath::Log10(maxPt);\r
- Double_t binwidth = (logmaxPt-logminPt)/nbinsPt;\r
- Double_t *binsPt = new Double_t[nbinsPt+1];\r
- binsPt[0] = minPt;\r
- for (Int_t i=1;i<=nbinsPt;i++) {\r
- binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth);\r
- }\r
-\r
- // 1pT resol cov matrix bins\r
- Double_t min1PtRes = 0.; \r
- Double_t max1PtRes = 0.3; \r
- Int_t nbins1PtRes = 300; \r
- Double_t bins1PtRes[301];\r
- for (Int_t i=0;i<=nbins1PtRes;i++) {\r
- bins1PtRes[i] = min1PtRes + i*(max1PtRes-min1PtRes)/nbins1PtRes;\r
- }\r
-\r
- // phi bins\r
- Double_t minPhi = 0.; \r
- Double_t maxPhi = 6.5; \r
- Int_t nbinsPhi = 100; \r
- Double_t binsPhi[101];\r
- for (Int_t i=0;i<=nbinsPhi;i++) {\r
- binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi;\r
- }\r
-\r
- // eta bins\r
- Double_t minEta = -1.;\r
- Double_t maxEta = 1.;\r
- Int_t nbinsEta = 20;\r
- Double_t binsEta[21];\r
- for (Int_t i=0;i<=nbinsEta;i++) {\r
- binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta;\r
- }\r
-\r
- // mult bins\r
- Double_t minCent = 0.;\r
- Double_t maxCent = 100;\r
- Int_t nbinsCent = 20;\r
- Double_t binsCent[101];\r
- for (Int_t i=0;i<=nbinsCent;i++) {\r
- binsCent[i] = minCent + i*(maxCent-minCent)/nbinsCent;\r
- }\r
- \r
- fPtResPhiPtTPC = new TH3D("fPtResPhiPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);\r
- fPtResPhiPtTPCc = new TH3D("fPtResPhiPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);\r
- fPtResPhiPtTPCITS = new TH3D("fPtResPhiPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);\r
- \r
-fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);\r
- fPtResEtaPtTPCc = new TH3D("fPtResEtaPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);\r
- fPtResEtaPtTPCITS = new TH3D("fPtResEtaPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);\r
- \r
-fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);\r
- fPtResCentPtTPCc = new TH3D("fPtResCentPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);\r
- fPtResCentPtTPCITS = new TH3D("fPtResCentPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);\r
-\r
- \r
- fOutput = new TList; \r
- if(!fOutput) return;\r
- fOutput->SetOwner();\r
-\r
- fOutput->Add(fPtResPhiPtTPC);\r
- fOutput->Add(fPtResPhiPtTPCc);\r
- fOutput->Add(fPtResPhiPtTPCITS);\r
- fOutput->Add(fPtResEtaPtTPC);\r
- fOutput->Add(fPtResEtaPtTPCc);\r
- fOutput->Add(fPtResEtaPtTPCITS);\r
- fOutput->Add(fPtResCentPtTPC);\r
- fOutput->Add(fPtResCentPtTPCc);\r
- fOutput->Add(fPtResCentPtTPCITS);\r
-\r
- // post data to outputs\r
-\r
- PostData(1,fV0Tree);\r
- PostData(2,fHighPtTree);\r
- PostData(3,fdEdxTree);\r
- PostData(4,fLaserTree);\r
- PostData(5,fMCEffTree);\r
- PostData(6,fCosmicPairsTree);\r
-\r
- PostData(7,fOutput);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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 MC info available - use it.\r
- fMC = MCEvent();\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
- }\r
- else {\r
- Process(fESD,fMC,fESDfriend); // only global and TPC tracks\r
- }\r
-\r
- //\r
- ProcessV0(fESD,fMC,fESDfriend);\r
- ProcessLaser(fESD,fMC,fESDfriend);\r
- ProcessdEdx(fESD,fMC,fESDfriend);\r
-\r
- if (fProcessCosmics) { ProcessCosmics(fESD); }\r
- if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);}\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event)\r
-{\r
- //\r
- // Select real events with high-pT tracks \r
- //\r
- if(!event) {\r
- AliDebug(AliLog::kError, "event not available");\r
- return;\r
- }\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
-\r
- // check for cosmic pairs\r
- //\r
- // find cosmic pairs trigger by random trigger\r
- //\r
- //\r
- AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();\r
- AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); \r
- const Double_t kMinPt=0.8;\r
- const Double_t kMinPtMax=0.8;\r
- const Double_t kMinNcl=50;\r
- const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};\r
- Int_t ntracks=event->GetNumberOfTracks(); \r
- // Float_t dcaTPC[2]={0,0};\r
- // Float_t covTPC[3]={0,0,0};\r
-\r
- UInt_t specie = event->GetEventSpecie(); // skip laser events\r
- if (specie==AliRecoParam::kCalib) return;\r
- \r
-\r
-\r
- for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {\r
- AliESDtrack *track0 = event->GetTrack(itrack0);\r
- if (!track0) continue;\r
- if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue;\r
-\r
- if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue;\r
- if (track0->Pt() < kMinPt) continue;\r
- if (track0->GetTPCncls() < kMinNcl) continue;\r
- if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue; \r
- if (track0->GetKinkIndex(0)>0) continue;\r
- const Double_t * par0=track0->GetParameter(); //track param at rhe DCA\r
- //rm primaries\r
- //\r
- //track0->GetImpactParametersTPC(dcaTPC,covTPC);\r
- //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;\r
- //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;\r
- // const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();\r
- for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {\r
- AliESDtrack *track1 = event->GetTrack(itrack1);\r
- if (!track1) continue; \r
- if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue;\r
- if (track1->GetKinkIndex(0)>0) continue;\r
- if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue;\r
- if (track1->Pt() < kMinPt) continue;\r
- if (track1->GetTPCncls()<kMinNcl) continue;\r
- if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;\r
- if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;\r
- //track1->GetImpactParametersTPC(dcaTPC,covTPC);\r
- // if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;\r
- //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;\r
- //\r
- const Double_t* par1=track1->GetParameter(); //track param at rhe DCA\r
- //\r
- Bool_t isPair=kTRUE;\r
- for (Int_t ipar=0; ipar<5; ipar++){\r
- if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off\r
- if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;\r
- }\r
- if (!isPair) continue;\r
- if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;\r
- //delta with correct sign\r
- /*\r
- TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"\r
- TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"\r
- TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"\r
- */\r
- if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign\r
- if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign\r
- if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign\r
- if (!isPair) continue;\r
- TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());\r
- Int_t eventNumber = event->GetEventNumberInFile(); \r
- //Bool_t hasFriend = kFALSE;\r
- //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);\r
- //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);\r
- // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); \r
- //\r
- // \r
- Int_t ntracksSPD = vertexSPD->GetNContributors();\r
- Int_t ntracksTPC = vertexTPC->GetNContributors(); \r
- Int_t runNumber = event->GetRunNumber(); \r
- Int_t timeStamp = event->GetTimeStamp();\r
- ULong64_t triggerMask = event->GetTriggerMask();\r
- Float_t magField = event->GetMagneticField();\r
- TObjString triggerClass = event->GetFiredTriggerClasses().Data();\r
- \r
- //\r
- // Dump to the tree \r
- // vertex\r
- // TPC-ITS tracks\r
- //\r
-\r
- //fCosmicPairsTree->Branch("fileName",&fileName,32000,0);\r
- //fCosmicPairsTree->Branch("runNumber",&runNumber,"runNumber/I");\r
- //fCosmicPairsTree->Branch("timeStamp",&timeStamp,"timeStamp/I");\r
- //fCosmicPairsTree->Branch("eventNumber",&eventNumber,"eventNumber/I");\r
- //fCosmicPairsTree->Branch("triggerMask",&triggerMask,32000,0);\r
- //fCosmicPairsTree->Branch("triggerClass",&triggerClass,32000,0);\r
- //fCosmicPairsTree->Branch("magField",&magField,"magField/F");\r
- //fCosmicPairsTree->Branch("ntracksSPD",&ntracksSPD,"ntracksSPD/I");\r
- //fCosmicPairsTree->Branch("ntracksTPC",&ntracksTPC,"ntracksTPC/I");\r
- //fCosmicPairsTree->Branch("vertexSPD",vertexSPD,32000,0);\r
- //fCosmicPairsTree->Branch("vertexTPC",vertexTPC,32000,0);\r
- //fCosmicPairsTree->Branch("track0",track0,32000,0);\r
- //fCosmicPairsTree->Branch("track1",track1,32000,0);\r
- //\r
- //fCosmicPairsTree->Fill();\r
-\r
- if(!fFillTree) return;\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"CosmicPairs"<<\r
- "fileName.="<<&fileName<< // file name\r
- "runNumber="<<runNumber<< // run number \r
- "evtTimeStamp="<<timeStamp<< // time stamp of event\r
- "evtNumberInFile="<<eventNumber<< // event number \r
- "trigger="<<triggerMask<< // trigger\r
- "triggerClass="<<&triggerClass<< // trigger\r
- "Bz="<<magField<< // magnetic field\r
- //\r
- "multSPD="<<ntracksSPD<<\r
- "multTPC="<<ntracksTPC<<\r
- "vertSPD.="<<vertexSPD<< //primary vertex -SPD\r
- "vertTPC.="<<vertexTPC<< //primary vertex -TPC\r
- "t0.="<<track0<< //track0\r
- "t1.="<<track1<< //track1\r
- "\n"; \r
- }\r
- }\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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
- AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
- AliFilteredTreeAcceptanceCuts *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(mcEvent)\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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
- } // end bUseMC\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == kTPCAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
- vtxESD = (AliESDVertex*)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, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv());\r
- //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
- Int_t ntracks = esdEvent->GetNumberOfTracks();\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
- Float_t bz = esdEvent->GetMagneticField();\r
- Int_t runNumber = esdEvent->GetRunNumber();\r
- Int_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
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
-\r
- //fHighPtTree->Branch("fileName",&fileName,32000,0);\r
- //fHighPtTree->Branch("runNumber",&runNumber,"runNumber/I");\r
- //fHighPtTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I");\r
- //fHighPtTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I");\r
- //fHighPtTree->Branch("triggerClass",&triggerClass,32000,0);\r
- //fHighPtTree->Branch("Bz",&bz,"Bz/F");\r
- //fHighPtTree->Branch("vtxESD",vtxESD,32000,0);\r
- //fHighPtTree->Branch("IRtot",&ir1,"IRtot/I");\r
- //fHighPtTree->Branch("IRint2",&ir2,"IRint2/I");\r
- //fHighPtTree->Branch("mult",&mult,"mult/I");\r
- //fHighPtTree->Branch("esdTrack",track,32000,0);\r
- //fHighPtTree->Branch("centralityF",¢ralityF,"centralityF/F");\r
-\r
- //fHighPtTree->Fill();\r
-\r
- //Double_t vtxX=vtxESD->GetX();\r
- //Double_t vtxY=vtxESD->GetY();\r
- //Double_t vtxZ=vtxESD->GetZ();\r
- if(!fFillTree) return;\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"highPt"<<\r
- "fileName.="<<&fileName<< \r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "triggerClass="<<&triggerClass<< // trigger\r
- "Bz="<<bz<< // magnetic field\r
- "vtxESD.="<<vtxESD<<\r
- "ntracksESD="<<ntracks<< // number of tracks in the ESD\r
- // "vtxESDx="<<vtxX<<\r
- // "vtxESDy="<<vtxY<<\r
- // "vtxESDz="<<vtxZ<<\r
- "IRtot="<<ir1<< // interaction record history info\r
- "IRint2="<<ir2<<\r
- "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex\r
- "esdTrack.="<<track<<\r
- "centralityF="<<centralityF<< \r
- "\n";\r
- }\r
- }\r
- \r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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
- Int_t runNumber = esdEvent->GetRunNumber();\r
- Int_t evtTimeStamp = esdEvent->GetTimeStamp();\r
- Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();\r
- Float_t bz = esdEvent->GetMagneticField();\r
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
-\r
- //fLaserTree->Branch("fileName",&fileName,32000,0);\r
- //fLaserTree->Branch("runNumber",&runNumber,"runNumber/I");\r
- //fLaserTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I");\r
- //fLaserTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I");\r
- //fLaserTree->Branch("triggerClass",&triggerClass,32000,0);\r
- //fLaserTree->Branch("Bz",&bz,"Bz/F");\r
- //fLaserTree->Branch("multTPCtracks",&countLaserTracks,"multTPCtracks/I");\r
-\r
- //fLaserTree->Fill();\r
-\r
- if(!fFillTree) return;\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"Laser"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "triggerClass="<<&triggerClass<< // trigger\r
- "Bz="<<bz<<\r
- "multTPCtracks="<<countLaserTracks<<\r
- "\n";\r
- }\r
- }\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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
- AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
- AliFilteredTreeAcceptanceCuts *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) {AliInfo("no physics selection"); return;}\r
- //SetPhysicsTriggerSelection(physicsSelection);\r
-\r
- if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
- // set trigger (V0AND)\r
- triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
- if(!triggerAnalysis) {AliInfo("no trigger analysis");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
- Int_t mcStackSize=0;\r
-\r
- Int_t multMCTrueTracks = 0;\r
- if(mcEvent)\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
- mcStackSize=stack->GetNtrack();\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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
-\r
- } // end bUseMC\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == kTPCAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();\r
- }\r
- else {\r
- AliInfo("no ESD vertex");\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) {AliInfo("no 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
- Float_t bz = esdEvent->GetMagneticField();\r
- Int_t runNumber = esdEvent->GetRunNumber();\r
- Int_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 *refTPCOut = NULL;\r
- AliTrackReference *refITS = NULL;\r
- AliTrackReference *refTRD = NULL;\r
- AliTrackReference *refTOF = NULL;\r
- AliTrackReference *refEMCAL = NULL;\r
- AliTrackReference *refPHOS = NULL;\r
- Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;\r
-\r
- Bool_t isOKtrackInnerC3 = kFALSE;\r
- AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));\r
-\r
- if(mcEvent) \r
- {\r
- if(!stack) return;\r
-\r
- //\r
- // global track\r
- //\r
- Int_t label = TMath::Abs(track->GetLabel()); \r
- if (label >= mcStackSize) continue;\r
- particle = stack->Particle(label);\r
- if (!particle) continue;\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
- if (labelTPC >= mcStackSize) continue;\r
- particleTPC = stack->Particle(labelTPC);\r
- if (!particleTPC) continue;\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(TMath::Abs(labelTPC), 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->Label()==label && ref->DetectorId()==AliTrackReference::kITS)\r
- {\r
- nrefITS++;\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->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)\r
- {\r
- nrefTPC++;\r
- refTPCOut=ref;\r
- if(!refTPCIn) {\r
- refTPCIn = ref;\r
- //printf("refTPCIn %p \n",refTPCIn);\r
- //break;\r
- }\r
- }\r
- // TRD\r
- if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)\r
- {\r
- nrefTRD++;\r
- if(!refTRD) {\r
- refTRD = ref;\r
- }\r
- }\r
- // TOF\r
- if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)\r
- {\r
- nrefTOF++;\r
- if(!refTOF) {\r
- refTOF = ref;\r
- }\r
- }\r
- // EMCAL\r
- if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)\r
- {\r
- nrefEMCAL++;\r
- if(!refEMCAL) {\r
- refEMCAL = ref;\r
- }\r
- }\r
- // PHOS\r
- // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)\r
-// {\r
-// nrefPHOS++;\r
-// if(!refPHOS) {\r
-// refPHOS = ref;\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
- if (labelITS >= mcStackSize) continue;\r
- particleITS = stack->Particle(labelITS);\r
- if (!particleITS) continue;\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(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;\r
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
- if (fReducePileUp){ \r
- //\r
- // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks\r
- // Done only in case no MC info \r
- //\r
- Float_t dcaTPC[2];\r
- track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);\r
- Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;\r
- Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));\r
- Bool_t keepPileUp=gRandom->Rndm()<0.05;\r
- if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){\r
- dumpToTree=kFALSE;\r
- }\r
- }\r
- /////////////////\r
- //book keeping of created dummy objects (to avoid NULL in trees)\r
- Bool_t newvtxESD=kFALSE;\r
- Bool_t newtrack=kFALSE;\r
- Bool_t newtpcInnerC=kFALSE;\r
- Bool_t newtrackInnerC=kFALSE;\r
- Bool_t newtrackInnerC2=kFALSE;\r
- Bool_t newouterITSc=kFALSE;\r
- Bool_t newtrackInnerC3=kFALSE;\r
- Bool_t newrefTPCIn=kFALSE;\r
- Bool_t newrefITS=kFALSE;\r
- Bool_t newparticle=kFALSE;\r
- Bool_t newparticleMother=kFALSE;\r
- Bool_t newparticleTPC=kFALSE;\r
- Bool_t newparticleMotherTPC=kFALSE;\r
- Bool_t newparticleITS=kFALSE;\r
- Bool_t newparticleMotherITS=kFALSE;\r
- \r
- //check that the vertex is there and that it is OK, \r
- //i.e. no null member arrays, otherwise a problem with merging\r
- //later on.\r
- //this is a very ugly hack!\r
- if (!vtxESD)\r
- {\r
- AliInfo("fixing the ESD vertex for streaming");\r
- vtxESD=new AliESDVertex(); \r
- //vtxESD->SetNContributors(1);\r
- //UShort_t pindices[1]; pindices[0]=0;\r
- //vtxESD->SetIndices(1,pindices);\r
- //vtxESD->SetNContributors(0);\r
- newvtxESD=kTRUE;\r
- }\r
- //\r
- if (!track) {track=new AliESDtrack();newtrack=kTRUE;}\r
- if (!tpcInnerC) {tpcInnerC=new AliExternalTrackParam();newtpcInnerC=kTRUE;}\r
- if (!trackInnerC) {trackInnerC=new AliExternalTrackParam();newtrackInnerC=kTRUE;}\r
- if (!trackInnerC2) {trackInnerC2=new AliExternalTrackParam();newtrackInnerC2=kTRUE;}\r
- if (!outerITSc) {outerITSc=new AliExternalTrackParam();newouterITSc=kTRUE;}\r
- if (!trackInnerC3) {trackInnerC3=new AliExternalTrackParam();newtrackInnerC3=kTRUE;}\r
- if (mcEvent)\r
- {\r
- if (!refTPCIn) {refTPCIn=new AliTrackReference(); newrefTPCIn=kTRUE;}\r
- if (!refITS) {refITS=new AliTrackReference();newrefITS=kTRUE;}\r
- if (!particle) {particle=new TParticle(); newparticle=kTRUE;}\r
- if (!particleMother) {particleMother=new TParticle();newparticleMother=kTRUE;}\r
- if (!particleTPC) {particleTPC=new TParticle();newparticleTPC=kTRUE;}\r
- if (!particleMotherTPC) {particleMotherTPC=new TParticle();newparticleMotherTPC=kTRUE;}\r
- if (!particleITS) {particleITS=new TParticle();newparticleITS=kTRUE;}\r
- if (!particleMotherITS) {particleMotherITS=new TParticle();newparticleMotherITS=kTRUE;}\r
- }\r
- /////////////////////////\r
-\r
- //Double_t vtxX=vtxESD->GetX();\r
- //Double_t vtxY=vtxESD->GetY();\r
- //Double_t vtxZ=vtxESD->GetZ();\r
-\r
- AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone();\r
- AliESDtrack* ptrack=(AliESDtrack*)track->Clone();\r
- AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone();\r
- AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone();\r
- AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone();\r
- AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone();\r
- AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone();\r
- Int_t ntracks = esdEvent->GetNumberOfTracks();\r
- \r
- // fill histograms\r
- FillHistograms(ptrack, ptpcInnerC, centralityF, (Double_t)chi2(0,0));\r
-\r
- if(fTreeSRedirector && dumpToTree && fFillTree) \r
- {\r
-\r
- (*fTreeSRedirector)<<"highPt"<<\r
- "fileName.="<<&fileName<< // name of the chunk file (hopefully full)\r
- "runNumber="<<runNumber<< // runNumber\r
- "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)\r
- "evtNumberInFile="<<evtNumberInFile<< // event number\r
- "triggerClass="<<&triggerClass<< // trigger class as a string\r
- "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)\r
- "vtxESD.="<<pvtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)\r
- //"vtxESDx="<<vtxX<<\r
- //"vtxESDy="<<vtxY<<\r
- //"vtxESDz="<<vtxZ<<\r
- "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1\r
- "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2\r
- "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex\r
- "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)\r
- "esdTrack.="<<ptrack<< // esdTrack as used in the physical analysis\r
- "extTPCInnerC.="<<ptpcInnerC<< // ??? \r
- "extInnerParamC.="<<ptrackInnerC<< // ???\r
- "extInnerParam.="<<ptrackInnerC2<< // ???\r
- "extOuterITS.="<<pouterITSc<< // ???\r
- "extInnerParamRef.="<<ptrackInnerC3<< // ???\r
- "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???\r
- "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined\r
- "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout\r
- "centralityF="<<centralityF;\r
- if (mcEvent)\r
- {\r
- AliTrackReference refDummy;\r
- if (!refITS) refITS = &refDummy;\r
- if (!refTRD) refTRD = &refDummy;\r
- if (!refTOF) refTOF = &refDummy;\r
- if (!refEMCAL) refEMCAL = &refDummy;\r
- if (!refPHOS) refPHOS = &refDummy;\r
- (*fTreeSRedirector)<<"highPt"<< \r
- "nrefITS="<<nrefITS<< // number of track references in the ITS\r
- "nrefTPC="<<nrefTPC<< // number of track references in the TPC\r
- "nrefTRD="<<nrefTRD<< // number of track references in the TRD\r
- "nrefTOF="<<nrefTOF<< // number of track references in the TOF\r
- "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF\r
- "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF\r
- "refTPCIn.="<<refTPCIn<<\r
- "refTPCOut.="<<refTPCOut<<\r
- "refITS.="<<refITS<< \r
- "refTRD.="<<refTRD<< \r
- "refTOF.="<<refTOF<< \r
- "refEMCAL.="<<refEMCAL<< \r
- "refPHOS.="<<refPHOS<< \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
- }\r
- //finish writing the entry\r
- AliInfo("writing tree highPt");\r
- (*fTreeSRedirector)<<"highPt"<<"\n";\r
- }\r
-\r
- delete pvtxESD;\r
- delete ptrack;\r
- delete ptpcInnerC;\r
- delete ptrackInnerC;\r
- delete ptrackInnerC2;\r
- delete pouterITSc;\r
- delete ptrackInnerC3;\r
-\r
- ////////////////////\r
- //delete the dummy objects we might have created.\r
- if (newvtxESD) {delete vtxESD; vtxESD=NULL;}\r
- if (newtrack) {delete track; track=NULL;}\r
- if (newtpcInnerC) {delete tpcInnerC; tpcInnerC=NULL;}\r
- if (newtrackInnerC) {delete trackInnerC; trackInnerC=NULL;}\r
- if (newtrackInnerC2) {delete trackInnerC2; trackInnerC2=NULL;}\r
- if (newouterITSc) {delete outerITSc; outerITSc=NULL;}\r
- if (newtrackInnerC3) {delete trackInnerC3; trackInnerC3=NULL;}\r
- if (newrefTPCIn) {delete refTPCIn; refTPCIn=NULL;}\r
- if (newrefITS) {delete refITS; refITS=NULL;}\r
- if (newparticle) {delete particle; particle=NULL;}\r
- if (newparticleMother) {delete particleMother; particleMother=NULL;}\r
- if (newparticleTPC) {delete particleTPC; particleTPC=NULL;}\r
- if (newparticleMotherTPC) {delete particleMotherTPC; particleMotherTPC=NULL;}\r
- if (newparticleITS) {delete particleITS; particleITS=NULL;}\r
- if (newparticleMotherITS) {delete particleMotherITS; particleMotherITS=NULL;}\r
- ///////////////\r
-\r
- delete tpcInnerC;\r
- delete trackInnerC;\r
- delete trackInnerC2;\r
- delete outerITSc;\r
- delete trackInnerC3;\r
- }\r
- }\r
-\r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Fill tree for efficiency studies MC only\r
- AliInfo("we start!");\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
- AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
- AliFilteredTreeAcceptanceCuts *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
- Int_t mcStackSize=0;\r
- TArrayF vtxMC(3);\r
-\r
- Int_t multMCTrueTracks = 0;\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
- mcStackSize=stack->GetNtrack();\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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
-\r
-\r
- // get reconstructed vertex \r
- //const AliESDVertex* vtxESD = 0; \r
- AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == kTPCAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
- vtxESD = (AliESDVertex*)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
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
-\r
- // check event cuts\r
- if(isEventOK && isEventTriggered)\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
- // loop over MC stack\r
- for (Int_t iMc = 0; iMc < mcStackSize; ++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
- Int_t isESDtrackCut= 0;\r
- Int_t isAccCuts = 0;\r
- AliESDtrack *recTrack = NULL; \r
-\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
- //\r
- Int_t label = TMath::Abs(track->GetLabel());\r
- if (label >= mcStackSize) continue;\r
- if(label == iMc) { \r
- Bool_t isAcc=esdTrackCuts->AcceptTrack(track);\r
- if (isAcc) isESDtrackCut=1;\r
- if (accCuts->AcceptTrack(track)) isAccCuts=1;\r
- isRec = kTRUE;\r
- if (recTrack){\r
- if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track\r
- if (!isAcc) continue;\r
- trackIndex = iTrack;\r
- }\r
- recTrack = esdEvent->GetTrack(trackIndex); \r
- continue;\r
- } \r
- }\r
- \r
- // Store information in the output tree\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 && fFillTree) {\r
- (*fTreeSRedirector)<<"MCEffTree"<<\r
- "fileName.="<<&fileName<<\r
- "triggerClass.="<<&triggerClass<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<< // \r
- "Bz="<<bz<< // magnetic field\r
- "vtxESD.="<<vtxESD<< // vertex info\r
- "mult="<<mult<< // primary vertex 9whatewe found) multiplicity\r
- "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities\r
- //\r
- "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts\r
- "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag\r
- "esdTrack.="<<recTrack<<\r
- "isRec="<<isRec<<\r
- "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection\r
- "particle.="<<particle<< // particle properties\r
- "particleMother.="<<particleMother<<\r
- "mech="<<mech<<\r
- "\n";\r
- }\r
-\r
- if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;\r
- }\r
- }\r
-\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AliAnalysisTaskFilteredTree::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 AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)\r
-{\r
- //\r
- // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates\r
- //\r
- if(!esdEvent) {\r
- AliDebug(AliLog::kError, "esdEvent not available");\r
- return;\r
- }\r
-\r
- // get selection cuts\r
- AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
- AliFilteredTreeAcceptanceCuts *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
- AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == kTPCAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
- vtxESD = (AliESDVertex*)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
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
-\r
- if(!fFillTree) return;\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"V0s"<<\r
- "isDownscaled="<<isDownscaled<<\r
- "triggerClass="<<&triggerClass<< // trigger\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
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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
- AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); \r
- AliFilteredTreeAcceptanceCuts *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
- AliESDVertex* vtxESD = 0; \r
- if(GetAnalysisMode() == kTPCAnalysisMode) {\r
- vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();\r
- }\r
- else if(GetAnalysisMode() == kTPCITSAnalysisMode) {\r
- vtxESD = (AliESDVertex*)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
- TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();\r
-\r
- if(!fFillTree) return;\r
- if(!fTreeSRedirector) return;\r
- (*fTreeSRedirector)<<"dEdx"<<\r
- "fileName.="<<&fileName<<\r
- "runNumber="<<runNumber<<\r
- "evtTimeStamp="<<evtTimeStamp<<\r
- "evtNumberInFile="<<evtNumberInFile<<\r
- "triggerClass="<<&triggerClass<< // trigger\r
- "Bz="<<bz<<\r
- "vtxESD.="<<vtxESD<<\r
- "mult="<<mult<<\r
- "esdTrack.="<<track<<\r
- "\n";\r
- }\r
- }\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Int_t AliAnalysisTaskFilteredTree::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 AliAnalysisTaskFilteredTree::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
- // Special treatment of the gamma conversion pt spectra is softer - \r
- Double_t mass00= v0->GetEffMass(0,0);\r
- const Double_t cutMass=0.2;\r
- if (TMath::Abs(mass00-0)<cutMass){\r
- downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate\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 AliAnalysisTaskFilteredTree::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
-Bool_t AliAnalysisTaskFilteredTree::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
-TParticle *AliAnalysisTaskFilteredTree::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
- Int_t mcStackSize=stack->GetNtrack();\r
- if (motherLabel>=mcStackSize) return NULL;\r
- mother = stack->Particle(motherLabel); \r
-\r
-return mother;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromConversion = kFALSE;\r
-\r
- if(stack) {\r
- Int_t mcStackSize=stack->GetNtrack();\r
- if (label>=mcStackSize) return kFALSE;\r
- TParticle* particle = stack->Particle(label);\r
- if (!particle) return kFALSE;\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
- if (motherLabel>=mcStackSize) return kFALSE;\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 AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromMaterial = kFALSE;\r
-\r
- if(stack) {\r
- Int_t mcStackSize=stack->GetNtrack();\r
- if (label>=mcStackSize) return kFALSE;\r
- TParticle* particle = stack->Particle(label);\r
- if (!particle) return kFALSE;\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
- if (motherLabel>=mcStackSize) return kFALSE;\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 AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
-{\r
- Bool_t isFromStrangeness = kFALSE;\r
-\r
- if(stack) {\r
- Int_t mcStackSize=stack->GetNtrack();\r
- if (label>=mcStackSize) return kFALSE;\r
- TParticle* particle = stack->Particle(label);\r
- if (!particle) return kFALSE;\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
- if (motherLabel>=mcStackSize) return kFALSE;\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 AliAnalysisTaskFilteredTree::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
- ////\r
- //chain = new TChain("highPt");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fHighPtTree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fHighPtTree) fHighPtTree->Print();\r
-\r
- ////\r
- //chain = new TChain("V0s");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fV0Tree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fV0Tree) fV0Tree->Print();\r
-\r
- ////\r
- //chain = new TChain("dEdx");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fdEdxTree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fdEdxTree) fdEdxTree->Print();\r
-\r
- ////\r
- //chain = new TChain("Laser");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fLaserTree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fLaserTree) fLaserTree->Print();\r
-\r
- ////\r
- //chain = new TChain("MCEffTree");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fMCEffTree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fMCEffTree) fMCEffTree->Print();\r
-\r
- ////\r
- //chain = new TChain("CosmicPairs");\r
- //if(chain) { \r
- // chain->Add("jotwinow_Temp_Trees.root");\r
- // fCosmicPairsTree = chain->CopyTree("1");\r
- // delete chain; chain=0; \r
- //}\r
- //if(fCosmicPairsTree) fCosmicPairsTree->Print(); \r
-\r
- Bool_t deleteTrees=kTRUE;\r
- if ((AliAnalysisManager::GetAnalysisManager()))\r
- {\r
- if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == \r
- AliAnalysisManager::kProofAnalysis)\r
- deleteTrees=kFALSE;\r
- }\r
- if (deleteTrees) delete fTreeSRedirector;\r
- fTreeSRedirector=NULL;\r
-\r
- //OpenFile(1);\r
-\r
- // Post output data.\r
- //PostData(1, fHighPtTree);\r
- //PostData(2, fV0Tree);\r
- //PostData(3, fdEdxTree);\r
- //PostData(4, fLaserTree);\r
- //PostData(5, fMCEffTree);\r
- //PostData(6, fCosmicPairsTree);\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::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("highPt");\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: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1));\r
- return;\r
- }\r
- */\r
- \r
-}\r
-\r
-//_____________________________________________________________________________\r
-Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)\r
-{\r
- //\r
- // calculate mc event true track multiplicity\r
- //\r
- if(!mcEvent) return 0;\r
-\r
- AliStack* stack = 0;\r
- Int_t mult = 0;\r
-\r
- // MC particle stack\r
- stack = mcEvent->Stack();\r
- if (!stack) return 0;\r
-\r
- //\r
- //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());\r
- //\r
-\r
- Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);\r
- if(!isEventOK) return 0; \r
-\r
- Int_t nPart = stack->GetNtrack();\r
- for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
- {\r
- TParticle* 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
- // physical primary\r
- Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
- if(!prim) continue;\r
-\r
- // checked accepted without pt cut\r
- //if(accCuts->AcceptTrack(particle)) \r
- if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) \r
- {\r
- mult++;\r
- }\r
- }\r
-\r
-return mult; \r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, const Double_t centralityF, const Double_t chi2TPCInnerC) \r
-{\r
-//\r
-// Fill pT relative resolution histograms for \r
-// TPC only, TPC only constrained to vertex and TPC+ITS tracking\r
-//\r
- if(!ptrack) return; \r
- if(!ptpcInnerC) return; \r
-\r
- const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();\r
- if(!innerParam) return;\r
-\r
- Float_t dxy, dz;\r
- ptrack->GetImpactParameters(dxy,dz);\r
-\r
-// TPC+ITS primary tracks \r
-if( abs(ptrack->Eta())<0.8 && \r
- ptrack->GetTPCClusterInfo(3,1)>120 && \r
- ptrack->IsOn(0x40) && \r
- ptrack->GetTPCclusters(0)>0.0 && \r
- ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && \r
- //abs(innerParam->GetX())>0.0 && \r
- //abs(innerParam->GetY()/innerParam->GetX())<0.14 && \r
- //abs(innerParam->GetTgl())<0.85 && \r
- ptrack->IsOn(0x0004) && \r
- ptrack->GetNcls(0)>0 &&\r
- ptrack->GetITSchi2()>0 && \r
- sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&\r
- sqrt(chi2TPCInnerC)<6 &&\r
- (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&\r
- abs(dz)<2.0 && \r
- abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )\r
- {\r
- fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));\r
- fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));\r
- fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));\r
- }\r
-\r
-// TPC primary tracks \r
-// and TPC constrained primary tracks \r
-\r
- AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam(); \r
- if(!ptpcInner) return;\r
-\r
-\r
- Float_t dxyTPC, dzTPC;\r
- ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);\r
-\r
-if( abs(ptrack->Eta())<0.8 && \r
- ptrack->GetTPCClusterInfo(3,1)>120 && \r
- ptrack->IsOn(0x40)&& \r
- ptrack->GetTPCclusters(0)>0.0 && \r
- ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && \r
- //abs(innerParam->GetX())>0.0 && \r
- //abs(innerParam->GetY()/innerParam->GetX())<0.14 && \r
- //abs(innerParam->GetTgl())<0.85 && \r
- abs(dzTPC)<3.2 && \r
- abs(dxyTPC)<2.4 )\r
- {\r
- // TPC only\r
- fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));\r
- fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));\r
- fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));\r
-\r
- // TPC constrained to vertex \r
- fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));\r
- fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));\r
- fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));\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. *
+ **************************************************************************/
+
+/*
+ Class to process/filter reconstruction information from ESD, ESD friends, MC and provide them for later reprocessing
+ Filtering schema - low pt part is downscaled - to have flat pt specatra of selected topologies (tracks and V0s)
+ Downscaling schema is controlled by downscaling factors
+ Usage:
+ 1.) Filtering on Lego train
+ 2.) expert QA for tracking (resolution efficnecy)
+ 3.) pt reoslution studies using V0s
+ 4.) dEdx calibration using V0s
+ 5.) pt resolution and dEdx studies using cosmic
+ +
+ 6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx)
+
+ Exported trees (with full objects and dereived variables):
+ 1.) "highPt" - filtered trees with esd tracks, derived variables(propagated tracks), optional MC info +optional space points
+ 2.) "V0s" - - filtered trees with selected V0s (rough KF chi2 cut), KF particle and corresponding esd tracks + optionla space points
+ 3.) "Laser" - dump laser tracks with space points if exests
+ 4.) "CosmicTree" - cosmic track candidate (random or triggered) + esdTracks(up/down)+ optional points
+ 5.) "dEdx" - tree with high dEdx tpc tracks
+*/
+
+#include "iostream"
+#include "TSystem.h"
+#include <TPDGCode.h>
+#include <TDatabasePDG.h>
+
+#include "TChain.h"
+#include "TTreeStream.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH3.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 "AliTrackPointArray.h"
+#include "AliSysInfo.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 "AliVTrack.h"
+#include "AliGeomManager.h"
+
+#include "AliCentrality.h"
+#include "AliESDVZERO.h"
+#include "AliMultiplicity.h"
+
+#include "AliESDtrackCuts.h"
+#include "AliMCEventHandler.h"
+#include "AliFilteredTreeEventCuts.h"
+#include "AliFilteredTreeAcceptanceCuts.h"
+
+#include "AliAnalysisTaskFilteredTree.h"
+#include "AliKFParticle.h"
+#include "AliESDv0.h"
+#include "TVectorD.h"
+
+using namespace std;
+
+ClassImp(AliAnalysisTaskFilteredTree)
+
+ //_____________________________________________________________________________
+ AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name)
+ : AliAnalysisTaskSE(name)
+ , fESD(0)
+ , fMC(0)
+ , fESDfriend(0)
+ , fOutput(0)
+ , fPitList(0)
+ , fUseMCInfo(kFALSE)
+ , fUseESDfriends(kFALSE)
+ , fReducePileUp(kTRUE)
+ , fFillTree(kTRUE)
+ , fFilteredTreeEventCuts(0)
+ , fFilteredTreeAcceptanceCuts(0)
+ , fFilteredTreeRecAcceptanceCuts(0)
+ , fEsdTrackCuts(0)
+ , fTrigger(AliTriggerAnalysis::kMB1)
+ , fAnalysisMode(kTPCAnalysisMode)
+ , fTreeSRedirector(0)
+ , fCentralityEstimator(0)
+ , fLowPtTrackDownscaligF(0)
+ , fLowPtV0DownscaligF(0)
+ , fFriendDownscaling(-3.)
+ , fProcessAll(kFALSE)
+ , fProcessCosmics(kFALSE)
+ , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks
+ , fHighPtTree(0)
+ , fV0Tree(0)
+ , fdEdxTree(0)
+ , fLaserTree(0)
+ , fMCEffTree(0)
+ , fCosmicPairsTree(0)
+ , fPtResPhiPtTPC(0)
+ , fPtResPhiPtTPCc(0)
+ , fPtResPhiPtTPCITS(0)
+ , fPtResEtaPtTPC(0)
+ , fPtResEtaPtTPCc(0)
+ , fPtResEtaPtTPCITS(0)
+ , fPtResCentPtTPC(0)
+ , fPtResCentPtTPCc(0)
+ , fPtResCentPtTPCITS(0)
+ , fCurrentFileName("")
+ , fDummyTrack(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());
+ DefineOutput(7, TList::Class());
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree()
+{
+ //
+ // Destructor
+ //
+ delete fFilteredTreeEventCuts;
+ delete fFilteredTreeAcceptanceCuts;
+ delete fFilteredTreeRecAcceptanceCuts;
+ delete fEsdTrackCuts;
+}
+
+//____________________________________________________________________________
+Bool_t AliAnalysisTaskFilteredTree::Notify()
+{
+ //
+ //
+ static Int_t count = 0;
+ count++;
+ TTree *chain = (TChain*)GetInputData(0);
+ if(chain){
+ Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());
+ }
+ fCurrentFileName=chain->GetCurrentFile()->GetName();
+ return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::UserCreateOutputObjects()
+{
+ // Create histograms
+ // Called once
+
+ //
+ //get the output file to make sure the trees will be associated to it
+ OpenFile(1);
+ fTreeSRedirector = new TTreeSRedirector();
+
+ //
+ // Create trees
+ fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree();
+ fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree();
+ fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree();
+ fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree();
+ fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree();
+ fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
+
+ if (!fDummyTrack) {
+ fDummyTrack=new AliESDtrack();
+ }
+
+ // histogram booking
+
+ Double_t minPt = 0.1;
+ Double_t maxPt = 100.;
+ Int_t nbinsPt = 30;
+
+ Double_t logminPt = TMath::Log10(minPt);
+ Double_t logmaxPt = TMath::Log10(maxPt);
+ Double_t binwidth = (logmaxPt-logminPt)/nbinsPt;
+ Double_t *binsPt = new Double_t[nbinsPt+1];
+ binsPt[0] = minPt;
+ for (Int_t i=1;i<=nbinsPt;i++) {
+ binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth);
+ }
+
+ // 1pT resol cov matrix bins
+ Double_t min1PtRes = 0.;
+ Double_t max1PtRes = 0.3;
+ Int_t nbins1PtRes = 300;
+ Double_t bins1PtRes[301];
+ for (Int_t i=0;i<=nbins1PtRes;i++) {
+ bins1PtRes[i] = min1PtRes + i*(max1PtRes-min1PtRes)/nbins1PtRes;
+ }
+
+ // phi bins
+ Double_t minPhi = 0.;
+ Double_t maxPhi = 6.5;
+ Int_t nbinsPhi = 100;
+ Double_t binsPhi[101];
+ for (Int_t i=0;i<=nbinsPhi;i++) {
+ binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi;
+ }
+
+ // eta bins
+ Double_t minEta = -1.;
+ Double_t maxEta = 1.;
+ Int_t nbinsEta = 20;
+ Double_t binsEta[21];
+ for (Int_t i=0;i<=nbinsEta;i++) {
+ binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta;
+ }
+
+ // mult bins
+ Double_t minCent = 0.;
+ Double_t maxCent = 100;
+ Int_t nbinsCent = 20;
+ Double_t binsCent[101];
+ for (Int_t i=0;i<=nbinsCent;i++) {
+ binsCent[i] = minCent + i*(maxCent-minCent)/nbinsCent;
+ }
+
+ fPtResPhiPtTPC = new TH3D("fPtResPhiPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
+ fPtResPhiPtTPCc = new TH3D("fPtResPhiPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
+ fPtResPhiPtTPCITS = new TH3D("fPtResPhiPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes);
+
+ fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
+ fPtResEtaPtTPCc = new TH3D("fPtResEtaPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
+ fPtResEtaPtTPCITS = new TH3D("fPtResEtaPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes);
+
+ fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
+ fPtResCentPtTPCc = new TH3D("fPtResCentPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
+ fPtResCentPtTPCITS = new TH3D("fPtResCentPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes);
+
+
+ fOutput = new TList;
+ if(!fOutput) return;
+ fOutput->SetOwner();
+
+ fOutput->Add(fPtResPhiPtTPC);
+ fOutput->Add(fPtResPhiPtTPCc);
+ fOutput->Add(fPtResPhiPtTPCITS);
+ fOutput->Add(fPtResEtaPtTPC);
+ fOutput->Add(fPtResEtaPtTPCc);
+ fOutput->Add(fPtResEtaPtTPCITS);
+ fOutput->Add(fPtResCentPtTPC);
+ fOutput->Add(fPtResCentPtTPCc);
+ fOutput->Add(fPtResCentPtTPCITS);
+
+ // post data to outputs
+
+ PostData(1,fV0Tree);
+ PostData(2,fHighPtTree);
+ PostData(3,fdEdxTree);
+ PostData(4,fLaserTree);
+ PostData(5,fMCEffTree);
+ PostData(6,fCosmicPairsTree);
+
+ PostData(7,fOutput);
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::UserExec(Option_t *)
+{
+ //
+ // Called for each event
+ //
+
+ // ESD event
+ fESD = (AliESDEvent*) (InputEvent());
+ if (!fESD) {
+ Printf("ERROR: ESD event not available");
+ return;
+ }
+ //if MC info available - use it.
+ fMC = MCEvent();
+ if(fUseESDfriends) {
+ //fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
+ fESDfriend = ESDfriend();
+ if(!fESDfriend) {
+ Printf("ERROR: ESD friends not available");
+ }
+ }
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if (!inputHandler){
+ return;
+ }
+
+ //if set, use the environment variables to set the downscaling factors
+ //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF
+ //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF
+ //AliAnalysisTaskFilteredTree_fFriendDownscaling
+ TString env;
+ env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF");
+ if (!env.IsNull()){
+ fLowPtTrackDownscaligF=env.Atof();
+ AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF));
+ }
+ env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF");
+ if (!env.IsNull()){
+ fLowPtV0DownscaligF=env.Atof();
+ AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF));
+ }
+ env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling");
+ if (!env.IsNull()){
+ fFriendDownscaling=env.Atof();
+ AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling));
+ }
+ //
+ //
+ //
+ 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,fESDfriend); }
+ if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);}
+ if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend);
+ printf("processed event %d\n", Int_t(Entry()));
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend)
+{
+ //
+ // Find cosmic pairs (triggered or random)
+ //
+ //
+ 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();
+ UInt_t specie = event->GetEventSpecie(); // skip laser events
+ if (specie==AliRecoParam::kCalib) return;
+ Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0;
+
+
+ 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();
+ AliESDfriendTrack* friendTrack0=NULL;
+ if (esdFriend &&!esdFriend->TestSkipBit()){
+ if (itrack0<ntracksFriend){
+ friendTrack0 = esdFriend->GetTrack(itrack0);
+ } //this guy can be NULL
+ }
+
+ 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();
+ //
+ //
+ 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();
+
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)event->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+
+
+ AliESDfriendTrack* friendTrack1=NULL;
+ if (esdFriend &&!esdFriend->TestSkipBit()){
+ if (itrack1<ntracksFriend){
+ friendTrack1 = esdFriend->GetTrack(itrack1);
+ } //this guy can be NULL
+ }
+
+ //
+ AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
+ AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
+ if (fFriendDownscaling>=1){ // downscaling number of friend tracks
+ if (gRandom->Rndm()>1./fFriendDownscaling){
+ friendTrackStore0 = 0;
+ friendTrackStore1 = 0;
+ }
+ }
+ if (fFriendDownscaling<=0){
+ if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){
+ TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree();
+ if (tree){
+ Double_t sizeAll=tree->GetZipBytes();
+ TBranch * br= tree->GetBranch("friendTrack0.fPoints");
+ Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
+ br= tree->GetBranch("friendTrack0.fCalibContainer");
+ if (br) sizeFriend+=br->GetZipBytes();
+ if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
+ friendTrackStore0=0;
+ friendTrackStore1=0;
+ }
+ }
+ }
+ }
+ if(!fFillTree) return;
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"CosmicPairs"<<
+ "gid="<<gid<< // global id of track
+ "fileName.="<<&fCurrentFileName<< // file name
+ "runNumber="<<runNumber<< // run number
+ "evtTimeStamp="<<timeStamp<< // time stamp of event
+ "evtNumberInFile="<<eventNumber<< // event number
+ "trigger="<<triggerMask<< // trigger mask
+ "triggerClass="<<&triggerClass<< // trigger class
+ "Bz="<<magField<< // magnetic field
+ //
+ "multSPD="<<ntracksSPD<< // event ultiplicity
+ "multTPC="<<ntracksTPC<< //
+ "vertSPD.="<<vertexSPD<< // primary vertex -SPD
+ "vertTPC.="<<vertexTPC<< // primary vertex -TPC
+ "t0.="<<track0<< // first half of comsic trak
+ "t1.="<<track1<< // second half of cosmic track
+ "friendTrack0.="<<friendTrackStore0<< // friend information first track + points
+ "friendTrack1.="<<friendTrackStore1<< // frined information first track + points
+ "\n";
+ }
+ }
+}
+
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Select real events with high-pT tracks
+ //
+
+ // get selection cuts
+ AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
+ AliFilteredTreeAcceptanceCuts *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();
+ // 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(mcEvent)
+ {
+ // 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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+ } // end bUseMC
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == kTPCAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD();
+
+ if(!vtxESD) return;
+ if(!vtxTPC) return;
+ if(!vtxSPD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ());
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+ Int_t ntracks = esdEvent->GetNumberOfTracks();
+
+ // 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->GetX();
+ //vert[1] = vtxESD->GetY();
+ //vert[2] = vtxESD->GetZ();
+ Int_t mult = vtxESD->GetNContributors();
+ Int_t multSPD = vtxSPD->GetNContributors();
+ Int_t multTPC = vtxTPC->GetNContributors();
+
+ Float_t bz = esdEvent->GetMagneticField();
+ Int_t runNumber = esdEvent->GetRunNumber();
+ Int_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+
+
+ // 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
+ //
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+ if(!fFillTree) return;
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"highPt"<<
+ "gid="<<gid<<
+ "fileName.="<<&fCurrentFileName<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "triggerClass="<<&triggerClass<< // trigger
+ "Bz="<<bz<< // magnetic field
+ "vtxESD.="<<vtxESD<<
+ "ntracksESD="<<ntracks<< // number of tracks in the ESD
+ "IRtot="<<ir1<< // interaction record history info
+ "IRint2="<<ir2<<
+ "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
+ "multSPD="<<multSPD<< // multiplicity of tracks pointing to the SPD primary vertex
+ "multTPC="<<multTPC<< // multiplicity of tracks pointing to the TPC primary vertex
+ "esdTrack.="<<track<<
+ "centralityF="<<centralityF<<
+ "\n";
+ }
+ }
+
+}
+
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::ProcessLaser(AliESDEvent *const esdEvent, AliMCEvent * const /*mcEvent*/, AliESDfriend *const esdFriend)
+{
+ //
+ // Process laser events -> dump tracks and clusters to the special tree
+ //
+ const Double_t kMinPt = 5;
+ if(!fFillTree) return;
+ if(!fTreeSRedirector) return;
+ const AliESDHeader* esdHeader = esdEvent->GetHeader();
+ if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) {
+ Int_t countLaserTracks = 0;
+ Int_t runNumber = esdEvent->GetRunNumber();
+ Int_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+ Float_t bz = esdEvent->GetMagneticField();
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->GetTPCInnerParam()) countLaserTracks++;
+ AliESDfriendTrack* friendTrack=NULL;
+ // suppress beam background and CE random reacks
+ if (track->GetInnerParam()->Pt()<kMinPt) continue;
+ Bool_t skipTrack=gRandom->Rndm()>1/(1+TMath::Abs(fFriendDownscaling));
+ if (skipTrack) continue;
+ if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
+ (*fTreeSRedirector)<<"Laser"<<
+ "gid="<<gid<< // global identifier of event
+ "fileName.="<<&fCurrentFileName<< //
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "triggerClass="<<&triggerClass<< // trigger
+ "Bz="<<bz<< // magnetic field
+ "multTPCtracks="<<countLaserTracks<< // multiplicity of tracks
+ "track.="<<track<< // track parameters
+ "friendTrack.="<<friendTrack<< // friend track information
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::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 ...
+ //
+ // get selection cuts
+ AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
+ AliFilteredTreeAcceptanceCuts *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();
+
+
+ // trigger
+ if(evtCuts->IsTriggerRequired())
+ {
+ // always MB
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
+
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+ if(!physicsSelection) {AliInfo("no physics selection"); return;}
+ //SetPhysicsTriggerSelection(physicsSelection);
+
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
+ // set trigger (V0AND)
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();
+ if(!triggerAnalysis) {AliInfo("no trigger analysis");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 mcStackSize=0;
+
+ Int_t multMCTrueTracks = 0;
+ if(mcEvent)
+ {
+ // 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;
+ }
+ mcStackSize=stack->GetNtrack();
+
+ // 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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+
+ } // end bUseMC
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == kTPCAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ AliInfo("no ESD vertex");
+ return;
+ }
+
+ if(!vtxESD) return;
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+
+ //
+ // Vertex info comparison and track multiplicity
+ //
+ AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
+ AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
+ Int_t contSPD = vertexSPD->GetNContributors();
+ Int_t contTPC = vertexTPC->GetNContributors();
+ TVectorD vertexPosTPC(3), vertexPosSPD(3);
+ vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
+ vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
+ Int_t ntracksTPC=0;
+ Int_t ntracksITS=0;
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
+ if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
+ }
+
+ //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) {AliInfo("no 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->GetX();
+ vert[1] = vtxESD->GetY();
+ vert[2] = vtxESD->GetZ();
+ Int_t mult = vtxESD->GetNContributors();
+ Float_t bz = esdEvent->GetMagneticField();
+ Int_t runNumber = esdEvent->GetRunNumber();
+ Int_t evtTimeStamp = esdEvent->GetTimeStamp();
+ Int_t evtNumberInFile = esdEvent->GetEventNumberInFile();
+
+ Int_t numberOfTracks=esdEvent->GetNumberOfTracks();
+ // high pT tracks
+ for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ AliESDfriendTrack* friendTrack=NULL;
+ Int_t numberOfFriendTracks=0;
+ if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks();
+ if (esdFriend && iTrack<numberOfFriendTracks) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL
+ 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())
+ {
+ // propagate ITSout to TPC inner wall
+ if(friendTrack)
+ {
+
+ outerITSc = NULL;
+ if (friendTrack->GetITSOut()) 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 *refTPCOut = NULL;
+ AliTrackReference *refITS = NULL;
+ AliTrackReference *refTRD = NULL;
+ AliTrackReference *refTOF = NULL;
+ AliTrackReference *refEMCAL = NULL;
+ AliTrackReference *refPHOS = NULL;
+ Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0;
+
+ Bool_t isOKtrackInnerC3 = kFALSE;
+ AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam()));
+ if(mcEvent && stack)
+ {
+ do //artificial loop (once) to make the continue statements jump out of the MC part
+ {
+ multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+ //
+ // global track
+ //
+ Int_t label = TMath::Abs(track->GetLabel());
+ if (label >= mcStackSize) continue;
+ particle = stack->Particle(label);
+ if (!particle) continue;
+ 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());
+ if (labelTPC >= mcStackSize) continue;
+ particleTPC = stack->Particle(labelTPC);
+ if (!particleTPC) continue;
+ 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(TMath::Abs(labelTPC), 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->Label()==label && ref->DetectorId()==AliTrackReference::kITS)
+ {
+ nrefITS++;
+ if(!refITS && countITS==2) {
+ refITS = ref;
+ //printf("refITS %p \n",refITS);
+ }
+ countITS++;
+ }
+
+ // TPC
+ if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC)
+ {
+ nrefTPC++;
+ refTPCOut=ref;
+ if(!refTPCIn) {
+ refTPCIn = ref;
+ //printf("refTPCIn %p \n",refTPCIn);
+ //break;
+ }
+ }
+ // TRD
+ if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD)
+ {
+ nrefTRD++;
+ if(!refTRD) {
+ refTRD = ref;
+ }
+ }
+ // TOF
+ if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF)
+ {
+ nrefTOF++;
+ if(!refTOF) {
+ refTOF = ref;
+ }
+ }
+ // EMCAL
+ if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL)
+ {
+ nrefEMCAL++;
+ if(!refEMCAL) {
+ refEMCAL = ref;
+ }
+ }
+ // PHOS
+ // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS)
+ // {
+ // nrefPHOS++;
+ // if(!refPHOS) {
+ // refPHOS = ref;
+ // }
+ // }
+ }
+
+ // 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());
+ if (labelITS >= mcStackSize) continue;
+ particleITS = stack->Particle(labelITS);
+ if (!particleITS) continue;
+ 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);
+ }
+ }
+ while (0);
+ }
+
+ //
+ Bool_t dumpToTree=kFALSE;
+
+ if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE;
+ //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
+ if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE;
+ if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE;
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+ if (fReducePileUp){
+ //
+ // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks
+ // Done only in case no MC info
+ //
+ Float_t dcaTPC[2];
+ track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]);
+ Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10;
+ Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin));
+ Bool_t keepPileUp=gRandom->Rndm()<0.05;
+ if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){
+ dumpToTree=kFALSE;
+ }
+ }
+
+ //init dummy objects
+ static AliESDVertex dummyvtxESD;
+ //if (!dummyvtxESD)
+ //{
+ // dummyvtxESD=new AliESDVertex();
+ // //dummyvtxESD->SetNContributors(1);
+ // //UShort_t pindices[1]; pindices[0]=0;
+ // //dummyvtxESD->SetIndices(1,pindices);
+ // //dummyvtxESD->SetNContributors(0);
+ //}
+ static AliExternalTrackParam dummyexternaltrackparam;
+ //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam();
+ static AliTrackReference dummytrackreference;
+ //if (!dummytrackreference) dummytrackreference=new AliTrackReference();
+ static TParticle dummyparticle;
+ //if (!dummyparticle) dummyparticle=new TParticle();
+
+ //assign the dummy objects if needed
+ if (!track) {track=fDummyTrack;}
+ AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing
+ if (fFriendDownscaling>=1){ // downscaling number of friend tracks
+ friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0;
+ }
+ if (fFriendDownscaling<=0){
+ if (((*fTreeSRedirector)<<"highPt").GetTree()){
+ TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree();
+ if (tree){
+ Double_t sizeAll=tree->GetZipBytes();
+ TBranch * br= tree->GetBranch("friendTrack.fPoints");
+ Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
+ br= tree->GetBranch("friendTrack.fCalibContainer");
+ if (br) sizeFriend+=br->GetZipBytes();
+ if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0;
+ }
+ }
+ }
+
+
+ // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;}
+ if (!vtxESD) {vtxESD=&dummyvtxESD;}
+ if (mcEvent)
+ {
+ if (!refTPCIn) {refTPCIn=&dummytrackreference;}
+ if (!refITS) {refITS=&dummytrackreference;}
+ if (!particle) {particle=&dummyparticle;}
+ if (!particleMother) {particleMother=&dummyparticle;}
+ if (!particleTPC) {particleTPC=&dummyparticle;}
+ if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;}
+ if (!particleITS) {particleITS=&dummyparticle;}
+ if (!particleMotherITS) {particleMotherITS=&dummyparticle;}
+ }
+
+ Int_t ntracks = esdEvent->GetNumberOfTracks();
+
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+
+ // fill histograms
+ FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0));
+
+ if(fTreeSRedirector && dumpToTree && fFillTree) {
+ (*fTreeSRedirector)<<"highPt"<<
+ "gid="<<gid<<
+ "fileName.="<<&fCurrentFileName<< // name of the chunk file (hopefully full)
+ "runNumber="<<runNumber<< // runNumber
+ "evtTimeStamp="<<evtTimeStamp<< // time stamp of event (in seconds)
+ "evtNumberInFile="<<evtNumberInFile<< // event number
+ "triggerClass="<<&triggerClass<< // trigger class as a string
+ "Bz="<<bz<< // solenoid magnetic field in the z direction (in kGaus)
+ "vtxESD.="<<vtxESD<< // vertexer ESD tracks (can be biased by TPC pileup tracks)
+ "IRtot="<<ir1<< // interaction record (trigger) counters - coutner 1
+ "IRint2="<<ir2<< // interaction record (trigger) coutners - counter 2
+ "mult="<<mult<< // multiplicity of tracks pointing to the primary vertex
+ "ntracks="<<ntracks<< // number of the esd tracks (to take into account the pileup in the TPC)
+ // important variables for the pile-up studies
+ "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
+ "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
+ "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
+ "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
+ "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
+ "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
+ //
+ "esdTrack.="<<track<< // esdTrack as used in the physical analysis
+ // "friendTrack.="<<friendTrack<< // esdFriendTrack associated to the esdTrack
+ "friendTrack.="<<friendTrackStore<< // esdFriendTrack associated to the esdTrack
+ "extTPCInnerC.="<<tpcInnerC<< // ???
+ "extInnerParamC.="<<trackInnerC<< // ???
+ "extInnerParam.="<<trackInnerC2<< // ???
+ "extOuterITS.="<<outerITSc<< // ???
+ "extInnerParamRef.="<<trackInnerC3<< // ???
+ "chi2TPCInnerC="<<chi2(0,0)<< // chi2 of tracks ???
+ "chi2InnerC="<<chi2trackC(0,0)<< // chi2s of tracks TPCinner to the combined
+ "chi2OuterITS="<<chi2OuterITS(0,0)<< // chi2s of tracks TPC at inner wall to the ITSout
+ "centralityF="<<centralityF;
+ if (mcEvent){
+ static AliTrackReference refDummy;
+ if (!refITS) refITS = &refDummy;
+ if (!refTRD) refTRD = &refDummy;
+ if (!refTOF) refTOF = &refDummy;
+ if (!refEMCAL) refEMCAL = &refDummy;
+ if (!refPHOS) refPHOS = &refDummy;
+ (*fTreeSRedirector)<<"highPt"<<
+ "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
+ "nrefITS="<<nrefITS<< // number of track references in the ITS
+ "nrefTPC="<<nrefTPC<< // number of track references in the TPC
+ "nrefTRD="<<nrefTRD<< // number of track references in the TRD
+ "nrefTOF="<<nrefTOF<< // number of track references in the TOF
+ "nrefEMCAL="<<nrefEMCAL<< // number of track references in the TOF
+ "nrefPHOS="<<nrefPHOS<< // number of track references in the TOF
+ "refTPCIn.="<<refTPCIn<<
+ "refTPCOut.="<<refTPCOut<<
+ "refITS.="<<refITS<<
+ "refTRD.="<<refTRD<<
+ "refTOF.="<<refTOF<<
+ "refEMCAL.="<<refEMCAL<<
+ "refPHOS.="<<refPHOS<<
+ "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;
+ }
+ //finish writing the entry
+ AliInfo("writing tree highPt");
+ (*fTreeSRedirector)<<"highPt"<<"\n";
+ }
+ AliSysInfo::AddStamp("filteringTask",iTrack,numberOfTracks,numberOfFriendTracks,(friendTrackStore)?0:1);
+ delete tpcInnerC;
+ delete trackInnerC;
+ delete trackInnerC2;
+ delete outerITSc;
+ delete trackInnerC3;
+ }
+ }
+}
+
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::ProcessMCEff(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const /*esdFriend*/)
+{
+ //
+ // Fill tree for efficiency studies MC only
+ AliInfo("we start!");
+ if(!mcEvent) {
+ AliDebug(AliLog::kError, "mcEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
+ AliFilteredTreeAcceptanceCuts *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();
+
+
+ // 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;
+ Int_t mcStackSize=0;
+ TArrayF vtxMC(3);
+
+ Int_t multMCTrueTracks = 0;
+ //
+ 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;
+ }
+ mcStackSize=stack->GetNtrack();
+
+ // 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 = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
+
+
+ // get reconstructed vertex
+ //const AliESDVertex* vtxESD = 0;
+ AliESDVertex* vtxESD = 0;
+ if(GetAnalysisMode() == kTPCAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks();
+ }
+ else {
+ return;
+ }
+
+ if(!vtxESD) return;
+ //
+ // Vertex info comparison and track multiplicity
+ //
+ AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD();
+ AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC();
+ Int_t contSPD = vertexSPD->GetNContributors();
+ Int_t contTPC = vertexTPC->GetNContributors();
+ TVectorD vertexPosTPC(3), vertexPosSPD(3);
+ vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray());
+ vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray());
+ Int_t ntracksTPC=0;
+ Int_t ntracksITS=0;
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++;
+ if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++;
+ }
+
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD);
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
+
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+
+ // check event cuts
+ if(isEventOK && isEventTriggered)
+ {
+ //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->GetX();
+ vert[1] = vtxESD->GetY();
+ vert[2] = vtxESD->GetZ();
+ 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 < mcStackSize; ++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;
+ Int_t trackLoopIndex = -1;
+ Int_t isESDtrackCut= 0;
+ Int_t isAccCuts = 0;
+ Int_t nRec = 0; // how many times reconstructed
+ Int_t nFakes = 0; // how many times reconstructed as a fake track
+ AliESDtrack *recTrack = NULL;
+
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if(track->Charge()==0) continue;
+ //
+ Int_t label = TMath::Abs(track->GetLabel());
+ if (label >= mcStackSize) continue;
+ if(label == iMc) {
+ Bool_t isAcc=esdTrackCuts->AcceptTrack(track);
+ if (isAcc) isESDtrackCut=1;
+ if (accCuts->AcceptTrack(track)) isAccCuts=1;
+ isRec = kTRUE;
+ trackIndex = iTrack;
+
+ if (recTrack){
+ if (track->GetTPCncls()<recTrack->GetTPCncls()) continue; // in case looper tracks use longer track
+ if (!isAcc) continue;
+ trackLoopIndex = iTrack;
+ }
+ recTrack = esdEvent->GetTrack(trackIndex);
+ nRec++;
+ if(track->GetLabel()<0) nFakes++;
+
+ continue;
+ }
+ }
+
+ // Store information in the output tree
+ if (trackLoopIndex>-1) {
+ recTrack = esdEvent->GetTrack(trackLoopIndex);
+ } else if (trackIndex >-1) {
+ recTrack = esdEvent->GetTrack(trackIndex);
+ } else {
+ recTrack = fDummyTrack;
+ }
+
+ 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 && fFillTree) {
+ (*fTreeSRedirector)<<"MCEffTree"<<
+ "fileName.="<<&fCurrentFileName<<
+ "triggerClass.="<<&triggerClass<<
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<< //
+ "Bz="<<bz<< // magnetic field
+ "vtxESD.="<<vtxESD<< // vertex info
+ //
+ "mult="<<mult<< // primary vertex 9whatewe found) multiplicity
+ "multMCTrueTracks="<<multMCTrueTracks<< // mC track multiplicities
+ // important variables for the pile-up studies
+ "contTPC="<< contTPC<< // number of contributors to the TPC primary vertex candidate
+ "contSPD="<< contSPD<< // number of contributors to the SPD primary vertex candidate
+ "vertexPosTPC.="<<&vertexPosTPC<< // TPC vertex position
+ "vertexPosSPD.="<<&vertexPosSPD<< // SPD vertex position
+ "ntracksTPC="<<ntracksTPC<< // total number of the TPC tracks which were refitted
+ "ntracksITS="<<ntracksITS<< // total number of the ITS tracks which were refitted
+ //
+ //
+ "isAcc0="<<isESDtrackCut<< // track accepted by ESD track cuts
+ "isAcc1="<<isAccCuts<< // track accepted by acceptance cuts flag
+ "esdTrack.="<<recTrack<< // reconstructed track (only the longest from the loopers)
+ "isRec="<<isRec<< // track was reconstructed
+ "tpcTrackLength="<<tpcTrackLength<< // track length in the TPC r projection
+ "particle.="<<particle<< // particle properties
+ "particleMother.="<<particleMother<< // particle mother
+ "mech="<<mech<< // production mechanizm
+ "nRec="<<nRec<< // how many times reconstruted
+ "nFakes="<<nFakes<< // how many times reconstructed as a fake track
+ "\n";
+ }
+
+ //if(trackIndex <0 && recTrack) delete recTrack; recTrack=0;
+ }
+ }
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskFilteredTree::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 AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)
+{
+ //
+ // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates
+ //
+ if(!esdEvent) {
+ AliDebug(AliLog::kError, "esdEvent not available");
+ return;
+ }
+
+ // get selection cuts
+ AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
+ AliFilteredTreeAcceptanceCuts *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();
+
+ // 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() == kTPCAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
+ 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();
+
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+
+
+ 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;
+ AliESDfriendTrack* friendTrack0=NULL;
+ AliESDfriendTrack* friendTrack1=NULL;
+ if (esdFriend) {
+ if (!esdFriend->TestSkipBit()){
+ Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
+ if (v0->GetIndex(0)<ntracksFriend){
+ friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL
+ }
+ if (v0->GetIndex(1)<ntracksFriend){
+ friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL
+ }
+ }
+ }
+ if (track0->GetSign()<0) {
+ track1 = esdEvent->GetTrack(v0->GetIndex(0));
+ track0 = esdEvent->GetTrack(v0->GetIndex(1));
+ }
+
+ //
+ AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing
+ AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing
+ if (fFriendDownscaling>=1){ // downscaling number of friend tracks
+ if (gRandom->Rndm()>1./fFriendDownscaling){
+ friendTrackStore0 = 0;
+ friendTrackStore1 = 0;
+ }
+ }
+ if (fFriendDownscaling<=0){
+ if (((*fTreeSRedirector)<<"V0s").GetTree()){
+ TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree();
+ if (tree){
+ Double_t sizeAll=tree->GetZipBytes();
+ TBranch * br= tree->GetBranch("friendTrack0.fPoints");
+ Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0;
+ br= tree->GetBranch("friendTrack0.fCalibContainer");
+ if (br) sizeFriend+=br->GetZipBytes();
+ if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) {
+ friendTrackStore0=0;
+ friendTrackStore1=0;
+ }
+ }
+ }
+ }
+
+ //
+ Bool_t isDownscaled = IsV0Downscaled(v0);
+ if (isDownscaled) continue;
+ AliKFParticle kfparticle; //
+ Int_t type=GetKFParticle(v0,esdEvent,kfparticle);
+ if (type==0) continue;
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+
+ if(!fFillTree) return;
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"V0s"<<
+ "gid="<<gid<< // global id of event
+ "isDownscaled="<<isDownscaled<< //
+ "triggerClass="<<&triggerClass<< // trigger
+ "Bz="<<bz<< //
+ "fileName.="<<&fCurrentFileName<< // full path - file name with ESD
+ "runNumber="<<run<< //
+ "evtTimeStamp="<<time<< // time stamp of event in secons
+ "evtNumberInFile="<<evNr<< //
+ "type="<<type<< // type of V0-
+ "ntracks="<<ntracks<<
+ "v0.="<<v0<<
+ "kf.="<<&kfparticle<<
+ "track0.="<<track0<< // track
+ "track1.="<<track1<<
+ "friendTrack0.="<<friendTrackStore0<<
+ "friendTrack1.="<<friendTrackStore1<<
+ "centralityF="<<centralityF<<
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::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
+ AliFilteredTreeEventCuts *evtCuts = GetEventCuts();
+ AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts();
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
+
+ if(!evtCuts || !accCuts || !esdTrackCuts) {
+ AliDebug(AliLog::kError, "cuts not available");
+ return;
+ }
+
+
+ // trigger
+ Bool_t isEventTriggered = kTRUE;
+ AliPhysicsSelection *physicsSelection = NULL;
+ AliTriggerAnalysis* triggerAnalysis = NULL;
+
+ //
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+
+ 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() == kTPCAnalysisMode) {
+ vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC();
+ }
+ else if(GetAnalysisMode() == kTPCITSAnalysisMode) {
+ 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->GetX();
+ vert[1] = vtxESD->GetY();
+ vert[2] = vtxESD->GetZ();
+ 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();
+
+ // Global event id calculation using orbitID, bunchCrossingID and periodID
+ ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber();
+ ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber();
+ ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber();
+ ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID);
+
+ // large dEdx
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
+ {
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ AliESDfriendTrack* friendTrack=NULL;
+ if (esdFriend && !esdFriend->TestSkipBit()) {
+ Int_t ntracksFriend = esdFriend->GetNumberOfTracks();
+ if (iTrack<ntracksFriend){
+ friendTrack = esdFriend->GetTrack(iTrack);
+ } //this guy can be NULL
+ }
+
+ if(track->Charge()==0) continue;
+ if(!esdTrackCuts->AcceptTrack(track)) continue;
+ if(!accCuts->AcceptTrack(track)) continue;
+
+ if(!IsHighDeDxParticle(track)) continue;
+ TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data();
+
+ if(!fFillTree) return;
+ if(!fTreeSRedirector) return;
+ (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree
+ "gid="<<gid<< // global id
+ "fileName.="<<&fCurrentFileName<< // file name
+ "runNumber="<<runNumber<<
+ "evtTimeStamp="<<evtTimeStamp<<
+ "evtNumberInFile="<<evtNumberInFile<<
+ "triggerClass="<<&triggerClass<< // trigger
+ "Bz="<<bz<<
+ "vtxESD.="<<vtxESD<< //
+ "mult="<<mult<<
+ "esdTrack.="<<track<<
+ "friendTrack.="<<friendTrack<<
+ "\n";
+ }
+ }
+}
+
+//_____________________________________________________________________________
+Int_t AliAnalysisTaskFilteredTree::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;
+ 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 AliAnalysisTaskFilteredTree::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;
+ //
+ // Special treatment of the gamma conversion pt spectra is softer -
+ Double_t mass00= v0->GetEffMass(0,0);
+ const Double_t cutMass=0.2;
+ if (TMath::Abs(mass00-0)<cutMass){
+ downscaleF/=10.; // 10 times smaller downscaling for the gamma concersion candidate
+ }
+ //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 AliAnalysisTaskFilteredTree::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 AliAnalysisTaskFilteredTree::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 *AliAnalysisTaskFilteredTree::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;
+ Int_t mcStackSize=stack->GetNtrack();
+ if (motherLabel>=mcStackSize) return NULL;
+ mother = stack->Particle(motherLabel);
+
+ return mother;
+}
+
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack)
+{
+ Bool_t isFromConversion = kFALSE;
+
+ if(stack) {
+ Int_t mcStackSize=stack->GetNtrack();
+ if (label>=mcStackSize) return kFALSE;
+ TParticle* particle = stack->Particle(label);
+ if (!particle) return kFALSE;
+
+ 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));
+ if (motherLabel>=mcStackSize) return kFALSE;
+ TParticle* mother = stack->Particle(motherLabel);
+ if(mother) {
+ Int_t motherPdg = mother->GetPdgCode();
+
+ if(!isPrim && mech==5 && motherPdg==kGamma) {
+ isFromConversion=kTRUE;
+ }
+ }
+ }
+ }
+
+ return isFromConversion;
+}
+
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack)
+{
+ Bool_t isFromMaterial = kFALSE;
+
+ if(stack) {
+ Int_t mcStackSize=stack->GetNtrack();
+ if (label>=mcStackSize) return kFALSE;
+ TParticle* particle = stack->Particle(label);
+ if (!particle) return kFALSE;
+
+ 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));
+ if (motherLabel>=mcStackSize) return kFALSE;
+ TParticle* mother = stack->Particle(motherLabel);
+ if(mother) {
+ if(!isPrim && mech==13) {
+ isFromMaterial=kTRUE;
+ }
+ }
+ }
+ }
+
+ return isFromMaterial;
+}
+
+//_____________________________________________________________________________
+Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack)
+{
+ Bool_t isFromStrangeness = kFALSE;
+
+ if(stack) {
+ Int_t mcStackSize=stack->GetNtrack();
+ if (label>=mcStackSize) return kFALSE;
+ TParticle* particle = stack->Particle(label);
+ if (!particle) return kFALSE;
+
+ 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));
+ if (motherLabel>=mcStackSize) return kFALSE;
+ 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 AliAnalysisTaskFilteredTree::FinishTaskOutput()
+{
+ //
+ // Called one at the end
+ // locally on working node
+ //
+ Bool_t deleteTrees=kTRUE;
+ if ((AliAnalysisManager::GetAnalysisManager()))
+ {
+ if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() ==
+ AliAnalysisManager::kProofAnalysis)
+ deleteTrees=kFALSE;
+ }
+ if (deleteTrees) delete fTreeSRedirector;
+ fTreeSRedirector=NULL;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::Terminate(Option_t *)
+{
+ //
+ // Called one at the end
+ //
+}
+
+//_____________________________________________________________________________
+Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts)
+{
+ //
+ // calculate mc event true track multiplicity
+ //
+ if(!mcEvent) return 0;
+
+ AliStack* stack = 0;
+ Int_t mult = 0;
+
+ // MC particle stack
+ stack = mcEvent->Stack();
+ if (!stack) return 0;
+
+ //
+ //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv());
+ //
+
+ Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent);
+ if(!isEventOK) return 0;
+
+ Int_t nPart = stack->GetNtrack();
+ for (Int_t iMc = 0; iMc < nPart; ++iMc)
+ {
+ TParticle* 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;
+
+ // physical primary
+ Bool_t prim = stack->IsPhysicalPrimary(iMc);
+ if(!prim) continue;
+
+ // checked accepted without pt cut
+ //if(accCuts->AcceptTrack(particle))
+ if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() )
+ {
+ mult++;
+ }
+ }
+
+ return mult;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC)
+{
+ //
+ // Fill pT relative resolution histograms for
+ // TPC only, TPC only constrained to vertex and TPC+ITS tracking
+ //
+ if(!ptrack) return;
+ if(!ptpcInnerC) return;
+
+ const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam();
+ if(!innerParam) return;
+
+ Float_t dxy, dz;
+ ptrack->GetImpactParameters(dxy,dz);
+
+ // TPC+ITS primary tracks
+ if( abs(ptrack->Eta())<0.8 &&
+ ptrack->GetTPCClusterInfo(3,1)>120 &&
+ ptrack->IsOn(0x40) &&
+ ptrack->GetTPCclusters(0)>0.0 &&
+ ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
+ //abs(innerParam->GetX())>0.0 &&
+ //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
+ //abs(innerParam->GetTgl())<0.85 &&
+ ptrack->IsOn(0x0004) &&
+ ptrack->GetNcls(0)>0 &&
+ ptrack->GetITSchi2()>0 &&
+ sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 &&
+ sqrt(chi2TPCInnerC)<6 &&
+ (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) &&
+ abs(dz)<2.0 &&
+ abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) )
+ {
+ fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
+ fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
+ fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2()));
+ }
+
+ // TPC primary tracks
+ // and TPC constrained primary tracks
+
+ AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam();
+ if(!ptpcInner) return;
+
+
+ Float_t dxyTPC, dzTPC;
+ ptrack->GetImpactParametersTPC(dxyTPC,dzTPC);
+
+ if( abs(ptrack->Eta())<0.8 &&
+ ptrack->GetTPCClusterInfo(3,1)>120 &&
+ ptrack->IsOn(0x40)&&
+ ptrack->GetTPCclusters(0)>0.0 &&
+ ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 &&
+ //abs(innerParam->GetX())>0.0 &&
+ //abs(innerParam->GetY()/innerParam->GetX())<0.14 &&
+ //abs(innerParam->GetTgl())<0.85 &&
+ abs(dzTPC)<3.2 &&
+ abs(dxyTPC)<2.4 )
+ {
+ // TPC only
+ fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
+ fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
+ fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2()));
+
+ // TPC constrained to vertex
+ fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
+ fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
+ fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2()));
+ }
+}
+
+
+void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){
+ //
+ // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks
+ // marian.ivanov@cern.ch
+ // 0.) Init variables
+ // 1.) GetTrack parameters at TPC inner wall
+ // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria
+ //
+ // Logic to be used in reco:
+ // 1.) Find matching ITSalone->TPCalone
+ // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone)
+ // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack
+ const Double_t radiusMatch=84.; // redius to propagate
+ //
+ const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut
+ const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut
+ const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut
+ const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut
+ const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut
+ const Double_t chi2Cut=100; // chi2 matching cut
+ //
+ if (!esdFriend) return; // not ITS standalone track
+ if (esdFriend->TestSkipBit()) return; // friends tracks not stored
+ Int_t ntracks=esdEvent->GetNumberOfTracks();
+ Float_t bz = esdEvent->GetMagneticField();
+ //
+ // 0.) Get parameters in reference radius TPC Inner wall
+ //
+ //
+ TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius
+ TMatrixD vecMomR0(ntracks,6); //
+ TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track
+ TMatrixD vecMomR1(ntracks,6); //
+ Double_t xyz[3], pxyz[3]; //
+ for (Int_t iTrack=0; iTrack<ntracks; iTrack++){
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);
+ if(!track) continue;
+ if (track->GetInnerParam()){
+ const AliExternalTrackParam *trackTPC=track->GetInnerParam();
+ trackTPC->GetXYZAt(radiusMatch,bz,xyz);
+ trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz);
+ for (Int_t i=0; i<3; i++){
+ vecPosR1(iTrack,i)=xyz[i];
+ vecMomR1(iTrack,i)=pxyz[i];
+ }
+ vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle
+ vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle
+ vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();;
+ vecMomR1(iTrack,5)= trackTPC->GetTgl();;
+ }
+ AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack);
+ if(!friendTrack) continue;
+ if (friendTrack->GetITSOut()){
+ const AliExternalTrackParam *trackITS=friendTrack->GetITSOut();
+ trackITS->GetXYZAt(radiusMatch,bz,xyz);
+ trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz);
+ for (Int_t i=0; i<3; i++){
+ vecPosR0(iTrack,i)=xyz[i];
+ vecMomR0(iTrack,i)=pxyz[i];
+ }
+ vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]);
+ vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]);
+ vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();;
+ vecMomR0(iTrack,5)= trackITS->GetTgl();;
+ }
+ }
+ //
+ // 1.) Find closest matching tracks, between the ITS standalone track
+ // and the all other tracks
+ // a.) caltegory - All
+ // b.) category - without ITS
+ //
+ //
+ Int_t ntracksPropagated=0;
+ AliExternalTrackParam extTrackDummy;
+ AliESDtrack esdTrackDummy;
+ AliExternalTrackParam itsAtTPC;
+ AliExternalTrackParam itsAtITSTPC;
+ for (Int_t iTrack0=0; iTrack0<ntracks; iTrack0++){
+ AliESDtrack *track0 = esdEvent->GetTrack(iTrack0);
+ if(!track0) continue;
+ if (track0->IsOn(AliVTrack::kTPCin)) continue;
+ AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0);
+ if (!friendTrack0) continue;
+ //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue;
+ //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone?
+ ntracksPropagated++;
+ //
+ // 2.) find clostest TPCtrack
+ // a.) all tracks
+ Double_t minChi2All=10000000;
+ Double_t minChi2TPC=10000000;
+ Double_t minChi2TPCITS=10000000;
+ Int_t indexAll=-1;
+ Int_t indexTPC=-1;
+ Int_t indexTPCITS=-1;
+ Int_t ncandidates0=0; // n candidates - rough cut
+ Int_t ncandidates1=0; // n candidates - rough + chi2 cut
+ itsAtTPC=*(friendTrack0->GetITSOut());
+ itsAtITSTPC=*(friendTrack0->GetITSOut());
+ for (Int_t iTrack1=0; iTrack1<ntracks; iTrack1++){
+ AliESDtrack *track1 = esdEvent->GetTrack(iTrack1);
+ if(!track1) continue;
+ if (!track1->IsOn(AliVTrack::kTPCin)) continue;
+ // fast checks
+ //
+ if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue;
+ if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue;
+ if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue;
+ if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue;
+ if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue;
+ ncandidates0++;
+ //
+ const AliExternalTrackParam * param1= track1->GetInnerParam();
+ if (!friendTrack0->GetITSOut()) continue;
+ AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut());
+ if (!outerITS.Rotate(param1->GetAlpha())) continue;
+ if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall
+ Double_t chi2 = outerITS.GetPredictedChi2(param1);
+ if (chi2>chi2Cut) continue;
+ ncandidates1++;
+ if (chi2<minChi2All){
+ minChi2All=chi2;
+ indexAll=iTrack1;
+ }
+ if (chi2<minChi2TPC && track1->IsOn(AliVTrack::kITSin)==0){
+ minChi2TPC=chi2;
+ indexTPC=iTrack1;
+ itsAtTPC=outerITS;
+ }
+ if (chi2<minChi2TPCITS && track1->IsOn(AliVTrack::kITSin)){
+ minChi2TPCITS=chi2;
+ indexTPCITS=iTrack1;
+ itsAtITSTPC=outerITS;
+ }
+ }
+ //
+ AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy;
+ AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy;
+ AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy;
+ (*fTreeSRedirector)<<"itsTPC"<<
+ "indexAll="<<indexAll<< // index of closest track (chi2)
+ "indexTPC="<<indexTPC<< // index of closest TPCalone tracks
+ "indexTPCITS="<<indexTPCITS<< // index of closest cobined tracks
+ "ncandidates0="<<ncandidates0<< // number of candidates
+ "ncandidates1="<<ncandidates1<<
+ //
+ "chi2All="<<minChi2All<< // chi2 of closest tracks
+ "chi2TPC="<<minChi2TPC<<
+ "chi2TPCITS="<<minChi2TPCITS<<
+ //
+ "track0.="<<track0<< // ITS standalone tracks
+ "trackAll.="<<trackAll<< // Closets other track
+ "trackTPC.="<<trackTPC<< // Closest TPC only track
+ "trackTPCITS.="<<trackTPCITS<< // closest combined track
+ //
+ "itsAtTPC.="<<&itsAtTPC<< // ITS track parameters at the TPC alone track frame
+ "itsAtITSTPC.="<<&itsAtITSTPC<< // ITS track parameters at the TPC combeined track frame
+ "\n";
+ }
+}
+
+void AliAnalysisTaskFilteredTree::ProcessTrackMatch(AliESDEvent *const /*esdEvent*/, AliESDfriend *const /*esdFriend*/){
+/*
+ Use TPC standalone, ITS standalone and combined tracks to categorize/resp. recover track information.
+
+ Track categories:
+ -TPC+ITS
+ -TPC only
+ -ITS only
+ Topology categories:
+ -Nice isolated tracks with full information
+ -Overlapped tracks - Refit and sign them
+ -Multiple found (check overlap factor) - Merge and sign
+ -Charge particle (kink) decays - Change of direction - Sign them)
+ Info:
+ -Array of indexes of closest tracks in each individual category
+ -Chi2 distance to the closest tracks at reference radius of TPCin
+ -Overlap factors - fraction of shared clusters or missing region
+ -Chi2 distance at DCA
+ Information matrix:
+ -matrix closest tracks from each categories
+ -characterization - chi2, index,chi2, overlap ratio
+
+ Decissions:
+ 0.) Kink decay or catastophic multiple scaterring
+ (combining all track categories)
+ - small chi2 at the DCA
+ - significantly large deflection angle
+ - Combinatorial algorithm - to decrease CPU time restriction of investigation to tracks with small DCA at refernce radius
+
+ 1.) if (TPConly && !(TPC+ITS) && ITSonly match ) {
+ if (!kink) TPCOnly.addoptITS
+ if (kink) TPConly sign
+ }
+
+ 2.) Overlap tracks - Refit with doUnfold
+
+*/
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff