From 1059f583e92280c291449252efdd040e48b27565 Mon Sep 17 00:00:00 2001 From: jotwinow Date: Fri, 19 Oct 2012 15:41:12 +0000 Subject: [PATCH] Class to fill high pT trees (M. Krzewicki) --- PWGPP/AliAnalysisTaskFilteredTree.cxx | 2120 +++++++++++++++++++++++ PWGPP/AliAnalysisTaskFilteredTree.h | 139 ++ PWGPP/AliFilteredTreeAcceptanceCuts.cxx | 177 ++ PWGPP/AliFilteredTreeAcceptanceCuts.h | 97 ++ PWGPP/AliFilteredTreeEventCuts.cxx | 475 +++++ PWGPP/AliFilteredTreeEventCuts.h | 130 ++ PWGPP/PWGPPLinkDef.h | 4 + PWGPP/macros/AddTaskFilteredTree.C | 161 ++ 8 files changed, 3303 insertions(+) create mode 100644 PWGPP/AliAnalysisTaskFilteredTree.cxx create mode 100644 PWGPP/AliAnalysisTaskFilteredTree.h create mode 100644 PWGPP/AliFilteredTreeAcceptanceCuts.cxx create mode 100644 PWGPP/AliFilteredTreeAcceptanceCuts.h create mode 100644 PWGPP/AliFilteredTreeEventCuts.cxx create mode 100644 PWGPP/AliFilteredTreeEventCuts.h create mode 100644 PWGPP/macros/AddTaskFilteredTree.C diff --git a/PWGPP/AliAnalysisTaskFilteredTree.cxx b/PWGPP/AliAnalysisTaskFilteredTree.cxx new file mode 100644 index 00000000000..2f3f263a531 --- /dev/null +++ b/PWGPP/AliAnalysisTaskFilteredTree.cxx @@ -0,0 +1,2120 @@ +/************************************************************************** +* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +* * +* Author: The ALICE Off-line Project. * +* Contributors are mentioned in the code where appropriate. * +* * +* Permission to use, copy, modify and distribute this software and its * +* documentation strictly for non-commercial purposes is hereby granted * +* without fee, provided that the above copyright notice appears in all * +* copies and that both the copyright notice and this permission notice * +* appear in the supporting documentation. The authors make no claims * +* about the suitability of this software for any purpose. It is * +* provided "as is" without express or implied warranty. * +**************************************************************************/ + +#include "iostream" + +#include +#include + +#include "TChain.h" +#include "TTreeStream.h" +#include "TTree.h" +#include "TH1F.h" +#include "TCanvas.h" +#include "TList.h" +#include "TObjArray.h" +#include "TFile.h" +#include "TMatrixD.h" +#include "TRandom3.h" + +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliInputEventHandler.h" +#include "AliStack.h" +#include "AliTrackReference.h" + +#include "AliPhysicsSelection.h" +#include "AliAnalysisTask.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliESDfriend.h" +#include "AliMCEvent.h" +#include "AliESDInputHandler.h" +#include "AliESDVertex.h" +#include "AliTracker.h" +#include "AliGeomManager.h" + +#include "AliCentrality.h" +#include "AliESDVZERO.h" +#include "AliMultiplicity.h" + +#include "AliESDtrackCuts.h" +#include "AliMCEventHandler.h" +#include "AliFilteredTreeEventCuts.h" +#include "AliFilteredTreeAcceptanceCuts.h" + +#include "AliAnalysisTaskFilteredTree.h" +#include "AliKFParticle.h" +#include "AliESDv0.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) + , fFilteredTreeEventCuts(0) + , fFilteredTreeAcceptanceCuts(0) + , fFilteredTreeRecAcceptanceCuts(0) + , fEsdTrackCuts(0) + , fTrigger(AliTriggerAnalysis::kMB1) + , fAnalysisMode(kTPCAnalysisMode) + , fTreeSRedirector(0) + , fCentralityEstimator(0) + , fLowPtTrackDownscaligF(0) + , fLowPtV0DownscaligF(0) + , fProcessAll(kFALSE) + , fProcessCosmics(kFALSE) + , fHighPtTree(0) + , fV0Tree(0) + , fdEdxTree(0) + , fLaserTree(0) + , fMCEffTree(0) + , fCosmicPairsTree(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()); +} + +//_____________________________________________________________________________ +AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree() +{ + if(fOutput) delete fOutput; fOutput =0; + if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0; + + if(fFilteredTreeEventCuts) delete fFilteredTreeEventCuts; fFilteredTreeEventCuts=NULL; + if(fFilteredTreeAcceptanceCuts) delete fFilteredTreeAcceptanceCuts; fFilteredTreeAcceptanceCuts=NULL; + if(fFilteredTreeRecAcceptanceCuts) delete fFilteredTreeRecAcceptanceCuts; fFilteredTreeRecAcceptanceCuts=NULL; + if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL; +} + +//____________________________________________________________________________ +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()); + } + +return kTRUE; +} + +//_____________________________________________________________________________ +void AliAnalysisTaskFilteredTree::UserCreateOutputObjects() +{ + // Create histograms + // Called once + fOutput = new TList; + fOutput->SetOwner(); + + // + // create temporary file for output tree + fTreeSRedirector = new TTreeSRedirector("jotwinow_Temp_Trees.root"); + + + + // + // Create trees + fHighPtTree = new TTree; + fV0Tree = new TTree; + fdEdxTree = new TTree; + fLaserTree = new TTree; + fMCEffTree = new TTree; + fCosmicPairsTree = new TTree; + + fOutput->Add(fHighPtTree); + fOutput->Add(fV0Tree); + fOutput->Add(fdEdxTree); + fOutput->Add(fLaserTree); + fOutput->Add(fMCEffTree); + fOutput->Add(fCosmicPairsTree); + + PostData(1,fHighPtTree); + PostData(2,fV0Tree); + PostData(3,fdEdxTree); + PostData(4,fLaserTree); + PostData(5,fMCEffTree); + PostData(6,fCosmicPairsTree); + +} + +//_____________________________________________________________________________ +void AliAnalysisTaskFilteredTree::UserExec(Option_t *) +{ + // + // Called for each event + // + + // ESD event + fESD = (AliESDEvent*) (InputEvent()); + if (!fESD) { + Printf("ERROR: ESD event not available"); + return; + } + + // MC event + if(fUseMCInfo) { + fMC = MCEvent(); + if (!fMC) { + Printf("ERROR: MC event not available"); + return; + } + } + + if(fUseESDfriends) { + fESDfriend = static_cast(fESD->FindListObject("AliESDfriend")); + if(!fESDfriend) { + Printf("ERROR: ESD friends not available"); + } + } + + // + if(fProcessAll) { + ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC + } + else { + Process(fESD,fMC,fESDfriend); // only global and TPC tracks + } + + // + ProcessV0(fESD,fMC,fESDfriend); + ProcessLaser(fESD,fMC,fESDfriend); + ProcessdEdx(fESD,fMC,fESDfriend); + + if (fProcessCosmics) { ProcessCosmics(fESD); } + if(IsUseMCInfo()) { ProcessMCEff(fESD,fMC,fESDfriend); } +} + +//_____________________________________________________________________________ +void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event) +{ + // + // Select real events with high-pT tracks + // + if(!event) { + AliDebug(AliLog::kError, "event not available"); + return; + } + + // + AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + // get file name + TTree *chain = (TChain*)GetInputData(0); + if(!chain) { + Printf("ERROR: Could not receive input chain"); + return; + } + TObjString fileName(chain->GetCurrentFile()->GetName()); + + + // check for cosmic pairs + // + // find cosmic pairs trigger by random trigger + // + // + AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); + AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); + const Double_t kMinPt=0.8; + const Double_t kMinPtMax=0.8; + const Double_t kMinNcl=50; + const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; + Int_t ntracks=event->GetNumberOfTracks(); + // Float_t dcaTPC[2]={0,0}; + // Float_t covTPC[3]={0,0,0}; + + UInt_t specie = event->GetEventSpecie(); // skip laser events + if (specie==AliRecoParam::kCalib) return; + + + + for (Int_t itrack0=0;itrack0GetTrack(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())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])GetInnerParam(); + for (Int_t itrack1=itrack0+1;itrack1GetTrack(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()1 && TMath::Max(track1->Pt(), track0->Pt())GetY())GetImpactParametersTPC(dcaTPC,covTPC); + // if (TMath::Abs(dcaTPC[0])GetParameter(); //track param at rhe DCA + // + Bool_t isPair=kTRUE; + for (Int_t ipar=0; ipar<5; ipar++){ + if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off + if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE; + } + if (!isPair) continue; + if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; + //delta with correct sign + /* + TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" + TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" + TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" + */ + if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign + if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign + if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign + if (!isPair) continue; + TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); + Int_t eventNumber = event->GetEventNumberInFile(); + //Bool_t hasFriend = kFALSE; + //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4); + //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS); + // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); + // + // + Int_t ntracksSPD = vertexSPD->GetNContributors(); + Int_t ntracksTPC = vertexTPC->GetNContributors(); + Int_t runNumber = event->GetRunNumber(); + Int_t timeStamp = event->GetTimeStamp(); + ULong64_t triggerMask = event->GetTriggerMask(); + Float_t magField = event->GetMagneticField(); + TObjString triggerClass = event->GetFiredTriggerClasses().Data(); + + // + // Dump to the tree + // vertex + // TPC-ITS tracks + // + if(!fTreeSRedirector) return; + (*fTreeSRedirector)<<"CosmicPairs"<< + "fileName.="<<&fileName<< // file name + "runNumber="<GetInputEventHandler(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + // get file name + TTree *chain = (TChain*)GetInputData(0); + if(!chain) { + Printf("ERROR: Could not receive input chain"); + return; + } + TObjString fileName(chain->GetCurrentFile()->GetName()); + + // trigger + if(evtCuts->IsTriggerRequired()) + { + // always MB + isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; + + physicsSelection = static_cast (inputHandler->GetEventSelection()); + if(!physicsSelection) return; + //SetPhysicsTriggerSelection(physicsSelection); + + if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { + // set trigger (V0AND) + triggerAnalysis = physicsSelection->GetTriggerAnalysis(); + if(!triggerAnalysis) return; + isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); + } + } + + // centrality determination + Float_t centralityF = -1; + AliCentrality *esdCentrality = esdEvent->GetCentrality(); + centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); + + // use MC information + AliHeader* header = 0; + AliGenEventHeader* genHeader = 0; + AliStack* stack = 0; + TArrayF vtxMC(3); + + Int_t multMCTrueTracks = 0; + if(IsUseMCInfo()) + { + // + if(!mcEvent) { + AliDebug(AliLog::kError, "mcEvent not available"); + return; + } + // get MC event header + header = mcEvent->Header(); + if (!header) { + AliDebug(AliLog::kError, "Header not available"); + return; + } + // MC particle stack + stack = mcEvent->Stack(); + if (!stack) { + AliDebug(AliLog::kError, "Stack not available"); + return; + } + + // get MC vertex + genHeader = header->GenEventHeader(); + if (!genHeader) { + AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); + return; + } + genHeader->PrimaryVertex(vtxMC); + + // multipliticy of all MC primary tracks + // in Zv, pt and eta ranges) + multMCTrueTracks = 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; + } + + if(!vtxESD) return; + + Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); + //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZv()); + //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); + + + // check event cuts + if(isEventOK && isEventTriggered) + { + + // + // get IR information + // + AliESDHeader *esdHeader = 0; + esdHeader = esdEvent->GetHeader(); + if(!esdHeader) return; + //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s + //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval + + // Use when Peter commit the changes in the header + Int_t ir1 = 0; + Int_t ir2 = 0; + + // + Double_t vert[3] = {0}; + vert[0] = vtxESD->GetXv(); + vert[1] = vtxESD->GetYv(); + vert[2] = vtxESD->GetZv(); + Int_t mult = vtxESD->GetNContributors(); + Double_t bz = esdEvent->GetMagneticField(); + Double_t runNumber = esdEvent->GetRunNumber(); + Double_t evtTimeStamp = esdEvent->GetTimeStamp(); + Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + + // high pT tracks + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + { + AliESDtrack *track = esdEvent->GetTrack(iTrack); + if(!track) continue; + if(track->Charge()==0) continue; + if(!esdTrackCuts->AcceptTrack(track)) continue; + if(!accCuts->AcceptTrack(track)) continue; + + // downscale low-pT tracks + Double_t scalempt= TMath::Min(track->Pt(),10.); + Double_t downscaleF = gRandom->Rndm(); + downscaleF *= fLowPtTrackDownscaligF; + if(TMath::Exp(2*scalempt)GetTPCInnerParam()); + if (!tpcInner) continue; + // transform to the track reference frame + Bool_t isOK = kFALSE; + isOK = tpcInner->Rotate(track->GetAlpha()); + isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); + if(!isOK) continue; + + // Dump to the tree + // vertex + // TPC-ITS tracks + // + if(!fTreeSRedirector) return; + (*fTreeSRedirector)<<"highPt"<< + "fileName.="<<&fileName<< + "runNumber="<GetCurrentFile()->GetName()); + + // laser events + const AliESDHeader* esdHeader = esdEvent->GetHeader(); + if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) + { + Int_t countLaserTracks = 0; + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + { + AliESDtrack *track = esdEvent->GetTrack(iTrack); + if(!track) continue; + + if(track->GetTPCInnerParam()) countLaserTracks++; + } + + if(countLaserTracks > 100) + { + Double_t runNumber = esdEvent->GetRunNumber(); + Double_t evtTimeStamp = esdEvent->GetTimeStamp(); + Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + Double_t bz = esdEvent->GetMagneticField(); + + if(!fTreeSRedirector) return; + (*fTreeSRedirector)<<"Laser"<< + "fileName.="<<&fileName<< + "runNumber="<GetInputEventHandler(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + // get file name + TTree *chain = (TChain*)GetInputData(0); + if(!chain) { + Printf("ERROR: Could not receive input chain"); + return; + } + TObjString fileName(chain->GetCurrentFile()->GetName()); + + // trigger + if(evtCuts->IsTriggerRequired()) + { + // always MB + isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; + + physicsSelection = static_cast (inputHandler->GetEventSelection()); + if(!physicsSelection) return; + //SetPhysicsTriggerSelection(physicsSelection); + + if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { + // set trigger (V0AND) + triggerAnalysis = physicsSelection->GetTriggerAnalysis(); + if(!triggerAnalysis) return; + isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); + } + } + + // centrality determination + Float_t centralityF = -1; + AliCentrality *esdCentrality = esdEvent->GetCentrality(); + centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); + + // use MC information + AliHeader* header = 0; + AliGenEventHeader* genHeader = 0; + AliStack* stack = 0; + TArrayF vtxMC(3); + + Int_t multMCTrueTracks = 0; + if(IsUseMCInfo()) + { + // + if(!mcEvent) { + AliDebug(AliLog::kError, "mcEvent not available"); + return; + } + // get MC event header + header = mcEvent->Header(); + if (!header) { + AliDebug(AliLog::kError, "Header not available"); + return; + } + // MC particle stack + stack = mcEvent->Stack(); + if (!stack) { + AliDebug(AliLog::kError, "Stack not available"); + return; + } + + // get MC vertex + genHeader = header->GenEventHeader(); + if (!genHeader) { + AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); + return; + } + genHeader->PrimaryVertex(vtxMC); + + // multipliticy of all MC primary tracks + // in Zv, pt and eta ranges) + multMCTrueTracks = 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; + } + + if(!vtxESD) return; + + Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); + //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); + //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); + + + // check event cuts + if(isEventOK && isEventTriggered) + { + // + // get IR information + // + AliESDHeader *esdHeader = 0; + esdHeader = esdEvent->GetHeader(); + if(!esdHeader) return; + //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s + //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval + // + Int_t ir1 = 0; + Int_t ir2 = 0; + + // + Double_t vert[3] = {0}; + vert[0] = vtxESD->GetXv(); + vert[1] = vtxESD->GetYv(); + vert[2] = vtxESD->GetZv(); + Int_t mult = vtxESD->GetNContributors(); + Double_t bz = esdEvent->GetMagneticField(); + Double_t runNumber = esdEvent->GetRunNumber(); + Double_t evtTimeStamp = esdEvent->GetTimeStamp(); + Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + + // high pT tracks + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + { + AliESDtrack *track = esdEvent->GetTrack(iTrack); + if(!track) continue; + if(track->Charge()==0) continue; + if(!esdTrackCuts->AcceptTrack(track)) continue; + if(!accCuts->AcceptTrack(track)) continue; + + // downscale low-pT tracks + Double_t scalempt= TMath::Min(track->Pt(),10.); + Double_t downscaleF = gRandom->Rndm(); + downscaleF *= fLowPtTrackDownscaligF; + if(TMath::Exp(2*scalempt)GetXYZ(x); + Double_t b[3]; AliTracker::GetBxByBz(x,b); + + // + // Transform TPC inner params to track reference frame + // + Bool_t isOKtpcInner = kFALSE; + AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam()); + if (tpcInner) { + // transform to the track reference frame + isOKtpcInner = tpcInner->Rotate(track->GetAlpha()); + isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); + } + + // + // Constrain TPC inner to vertex + // clone TPCinner has to be deleted + // + Bool_t isOKtpcInnerC = kFALSE; + AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam())); + if (tpcInnerC) { + isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b); + isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha()); + isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); + } + + // + // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex + // Clone track InnerParams has to be deleted + // + Bool_t isOKtrackInnerC = kFALSE; + AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam())); + if (trackInnerC) { + isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b); + isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha()); + isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); + } + + // + // calculate chi2 between vi and vj vectors + // with covi and covj covariance matrices + // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj) + // + TMatrixD deltaT(5,1), deltaTtrackC(5,1); + TMatrixD delta(1,5), deltatrackC(1,5); + TMatrixD covarM(5,5), covarMtrackC(5,5); + TMatrixD chi2(1,1); + TMatrixD chi2trackC(1,1); + + if(isOKtpcInnerC && isOKtrackInnerC) + { + for (Int_t ipar=0; ipar<5; ipar++) { + deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; + delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; + + deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; + deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; + + for (Int_t jpar=0; jpar<5; jpar++) { + Int_t index=track->GetIndex(ipar,jpar); + covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index]; + covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index]; + } + } + + // chi2 distance TPC constrained and TPC+ITS + TMatrixD covarMInv = covarM.Invert(); + TMatrixD mat2 = covarMInv*deltaT; + chi2 = delta*mat2; + //chi2.Print(); + + // chi2 distance TPC refitted constrained and TPC+ITS + TMatrixD covarMInvtrackC = covarMtrackC.Invert(); + TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC; + chi2trackC = deltatrackC*mat2trackC; + //chi2trackC.Print(); + } + + + // + // Propagate ITSout to TPC inner wall + // and calculate chi2 distance to track (InnerParams) + // + const Double_t kTPCRadius=85; + const Double_t kStep=3; + + // clone track InnerParams has to be deleted + Bool_t isOKtrackInnerC2 = kFALSE; + AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam())); + if (trackInnerC2) { + isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE); + } + + Bool_t isOKouterITSc = kFALSE; + AliExternalTrackParam *outerITSc = NULL; + TMatrixD chi2OuterITS(1,1); + + if(esdFriend && esdFriend->TestSkipBit()==kFALSE) + { + // propagate ITSout to TPC inner wall + AliESDfriendTrack *friendTrack = esdFriend->GetTrack(iTrack); + + if(friendTrack) + { + outerITSc = new AliExternalTrackParam(*friendTrack->GetITSOut()); + if(outerITSc) + { + isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE); + isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha()); + isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField()); + + // + // calculate chi2 between outerITS and innerParams + // cov without z-coordinate at the moment + // + TMatrixD deltaTouterITS(4,1); + TMatrixD deltaouterITS(1,4); + TMatrixD covarMouterITS(4,4); + + if(isOKtrackInnerC2 && isOKouterITSc) { + Int_t kipar = 0; + Int_t kjpar = 0; + for (Int_t ipar=0; ipar<5; ipar++) { + if(ipar!=1) { + deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar]; + deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar]; + } + + kjpar=0; + for (Int_t jpar=0; jpar<5; jpar++) { + Int_t index=outerITSc->GetIndex(ipar,jpar); + if(ipar !=1 || jpar!=1) { + covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index]; + } + if(jpar!=1) kjpar++; + } + if(ipar!=1) kipar++; + } + + // chi2 distance ITSout and InnerParams + TMatrixD covarMInvouterITS = covarMouterITS.Invert(); + TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS; + chi2OuterITS = deltaouterITS*mat2outerITS; + //chi2OuterITS.Print(); + } + } + } + } + + // + // MC info + // + TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL; + TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL; + Int_t mech=-1, mechTPC=-1, mechITS=-1; + Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE; + Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE; + Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE; + Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE; + + AliTrackReference *refTPCIn = NULL; + AliTrackReference *refITS = NULL; + + Bool_t isOKtrackInnerC3 = kFALSE; + AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam())); + + if(IsUseMCInfo()) + { + if(!stack) return; + + // + // global track + // + Int_t label = TMath::Abs(track->GetLabel()); + particle = stack->Particle(label); + if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.) + { + particleMother = GetMother(particle,stack); + mech = particle->GetUniqueID(); + isPrim = stack->IsPhysicalPrimary(label); + isFromStrangess = IsFromStrangeness(label,stack); + isFromConversion = IsFromConversion(label,stack); + isFromMaterial = IsFromMaterial(label,stack); + } + + // + // TPC track + // + Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); + particleTPC = stack->Particle(labelTPC); + if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.) + { + particleMotherTPC = GetMother(particleTPC,stack); + mechTPC = particleTPC->GetUniqueID(); + isPrimTPC = stack->IsPhysicalPrimary(labelTPC); + isFromStrangessTPC = IsFromStrangeness(labelTPC,stack); + isFromConversionTPC = IsFromConversion(labelTPC,stack); + isFromMaterialTPC = IsFromMaterial(labelTPC,stack); + } + + // + // store first track reference + // for TPC track + // + TParticle *part=0; + TClonesArray *trefs=0; + Int_t status = mcEvent->GetParticleAndTR(track->GetTPCLabel(), part, trefs); + + if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) + { + Int_t nTrackRef = trefs->GetEntries(); + //printf("nTrackRef %d \n",nTrackRef); + + Int_t countITS = 0; + for (Int_t iref = 0; iref < nTrackRef; iref++) + { + AliTrackReference *ref = (AliTrackReference *)trefs->At(iref); + + // Ref. in the middle ITS + if(ref && ref->DetectorId()==AliTrackReference::kITS) + { + if(!refITS && countITS==2) { + refITS = ref; + //printf("refITS %p \n",refITS); + } + countITS++; + } + + // TPC + if(ref && ref->DetectorId()==AliTrackReference::kTPC) + { + if(!refTPCIn) { + refTPCIn = ref; + //printf("refTPCIn %p \n",refTPCIn); + //break; + } + } + } + + // transform inner params to TrackRef + // reference frame + if(refTPCIn && trackInnerC3) + { + Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X()); + isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi); + isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE); + } + } + + // + // ITS track + // + Int_t labelITS = TMath::Abs(track->GetITSLabel()); + particleITS = stack->Particle(labelITS); + if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.) + { + particleMotherITS = GetMother(particleITS,stack); + mechITS = particleITS->GetUniqueID(); + isPrimITS = stack->IsPhysicalPrimary(labelITS); + isFromStrangessITS = IsFromStrangeness(labelITS,stack); + isFromConversionITS = IsFromConversion(labelITS,stack); + isFromMaterialITS = IsFromMaterial(labelITS,stack); + } + } + + // + Bool_t dumpToTree=kFALSE; + + if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE; + if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; + if(fUseMCInfo && isOKtrackInnerC3) dumpToTree = kTRUE; + + // + if(fTreeSRedirector && dumpToTree) + { + (*fTreeSRedirector)<<"highPt"<< + "fileName.="<<&fileName<< + "runNumber="<GetInputEventHandler(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + // get file name + TTree *chain = (TChain*)GetInputData(0); + if(!chain) { + Printf("ERROR: Could not receive input chain"); + return; + } + TObjString fileName(chain->GetCurrentFile()->GetName()); + + // trigger + if(evtCuts->IsTriggerRequired()) + { + // always MB + isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; + + physicsSelection = static_cast (inputHandler->GetEventSelection()); + if(!physicsSelection) return; + + if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { + // set trigger (V0AND) + triggerAnalysis = physicsSelection->GetTriggerAnalysis(); + if(!triggerAnalysis) return; + isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); + } + } + + // centrality determination + Float_t centralityF = -1; + AliCentrality *esdCentrality = esdEvent->GetCentrality(); + centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); + + // use MC information + AliHeader* header = 0; + AliGenEventHeader* genHeader = 0; + AliStack* stack = 0; + TArrayF vtxMC(3); + + Int_t multMCTrueTracks = 0; + if(IsUseMCInfo()) + { + // + if(!mcEvent) { + AliDebug(AliLog::kError, "mcEvent not available"); + return; + } + // get MC event header + header = mcEvent->Header(); + if (!header) { + AliDebug(AliLog::kError, "Header not available"); + return; + } + // MC particle stack + stack = mcEvent->Stack(); + if (!stack) { + AliDebug(AliLog::kError, "Stack not available"); + return; + } + + // get MC vertex + genHeader = header->GenEventHeader(); + if (!genHeader) { + AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); + return; + } + genHeader->PrimaryVertex(vtxMC); + + // multipliticy of all MC primary tracks + // in Zv, pt and eta ranges) + multMCTrueTracks = 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; + } + + if(!vtxESD) return; + + Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); + //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); + //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); + + // check event cuts + if(isEventOK && isEventTriggered) + { + if(IsUseMCInfo()) + { + if(!stack) return; + + // + // MC info + // + TParticle *particle=NULL; + TParticle *particleMother=NULL; + Int_t mech=-1; + + // reco event info + Double_t vert[3] = {0}; + vert[0] = vtxESD->GetXv(); + vert[1] = vtxESD->GetYv(); + vert[2] = vtxESD->GetZv(); + Int_t mult = vtxESD->GetNContributors(); + Double_t bz = esdEvent->GetMagneticField(); + Double_t runNumber = esdEvent->GetRunNumber(); + Double_t evtTimeStamp = esdEvent->GetTimeStamp(); + Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + + // loop over MC stack + for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) + { + particle = stack->Particle(iMc); + if (!particle) + continue; + + // only charged particles + if(!particle->GetPDG()) continue; + Double_t charge = particle->GetPDG()->Charge()/3.; + if (TMath::Abs(charge) < 0.001) + continue; + + // only primary particles + Bool_t prim = stack->IsPhysicalPrimary(iMc); + if(!prim) continue; + + // downscale low-pT particles + Double_t scalempt= TMath::Min(particle->Pt(),10.); + Double_t downscaleF = gRandom->Rndm(); + downscaleF *= fLowPtTrackDownscaligF; + if(TMath::Exp(2*scalempt)AcceptTrack(particle)) continue; + + // check if particle reconstructed + Bool_t isRec = kFALSE; + Int_t trackIndex = -1; + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + { + + AliESDtrack *track = esdEvent->GetTrack(iTrack); + if(!track) continue; + if(track->Charge()==0) continue; + if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) + { + Int_t label = TMath::Abs(track->GetLabel()); + if(label == iMc) { + isRec = kTRUE; + trackIndex = iTrack; + break; + } + } + } + + // Store information in the output tree + AliESDtrack *recTrack = NULL; + if(trackIndex>-1) { + recTrack = esdEvent->GetTrack(trackIndex); + } else { + recTrack = new AliESDtrack(); + } + + particleMother = GetMother(particle,stack); + mech = particle->GetUniqueID(); + + //MC particle track length + Double_t tpcTrackLength = 0.; + AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc); + if(mcParticle) { + Int_t counter; + tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0); + } + + + // + if(fTreeSRedirector) { + (*fTreeSRedirector)<<"MCEffTree"<< + "fileName.="<<&fileName<< + "runNumber="< 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) 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(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + // get file name + TTree *chain = (TChain*)GetInputData(0); + if(!chain) { + Printf("ERROR: Could not receive input chain"); + return; + } + TObjString fileName(chain->GetCurrentFile()->GetName()); + + // trigger + if(evtCuts->IsTriggerRequired()) + { + // always MB + isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; + + physicsSelection = static_cast (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(); + + + for (Int_t iv0=0; iv0GetV0(iv0); + if (!v0) continue; + AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0)); + AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1)); + if (!track0) continue; + if (!track1) continue; + if (track0->GetSign()<0) { + track1 = esdEvent->GetTrack(v0->GetIndex(0)); + track0 = esdEvent->GetTrack(v0->GetIndex(1)); + } + // + Bool_t isDownscaled = IsV0Downscaled(v0); + if (isDownscaled) continue; + AliKFParticle kfparticle; // + Int_t type=GetKFParticle(v0,esdEvent,kfparticle); + if (type==0) continue; + + if(!fTreeSRedirector) return; + (*fTreeSRedirector)<<"V0s"<< + "isDownscaled="<GetCurrentFile()->GetName()); + + // trigger + Bool_t isEventTriggered = kTRUE; + AliPhysicsSelection *physicsSelection = NULL; + AliTriggerAnalysis* triggerAnalysis = NULL; + + // + AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); + if (!inputHandler) + { + Printf("ERROR: Could not receive input handler"); + return; + } + + if(evtCuts->IsTriggerRequired()) + { + // always MB + isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; + + physicsSelection = static_cast (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->GetXv(); + vert[1] = vtxESD->GetYv(); + vert[2] = vtxESD->GetZv(); + Int_t mult = vtxESD->GetNContributors(); + Double_t bz = esdEvent->GetMagneticField(); + Double_t runNumber = esdEvent->GetRunNumber(); + Double_t evtTimeStamp = esdEvent->GetTimeStamp(); + Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); + + // large dEdx + for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) + { + AliESDtrack *track = esdEvent->GetTrack(iTrack); + if(!track) continue; + if(track->Charge()==0) continue; + if(!esdTrackCuts->AcceptTrack(track)) continue; + if(!accCuts->AcceptTrack(track)) continue; + + if(!IsHighDeDxParticle(track)) continue; + + if(!fTreeSRedirector) return; + (*fTreeSRedirector)<<"dEdx"<< + "fileName.="<<&fileName<< + "runNumber="<GetOnFlyStatus() ==kFALSE) return 0; + // + // 1.) track cut + // + AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0)); + AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1)); + /* + TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2"; + TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1"; + TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100"; + */ + if (TMath::Abs(track0->GetTgl())>1) return 0; + if (TMath::Abs(track1->GetTgl())>1) return 0; + if ((track0->GetTPCClusterInfo(2,1))<100) return 0; + if ((track1->GetTPCClusterInfo(2,1))<100) return 0; + //if ((track0->GetITSclusters(0))<2) return 0; + //if ((track1->GetITSclusters(0))<2) return 0; + Float_t pos0[2]={0}, cov0[3]={0}; + Float_t pos1[2]={0}, cov1[3]={0}; + track0->GetImpactParameters(pos0,cov0); + track0->GetImpactParameters(pos1,cov1); + // + if (TMath::Abs(pos0[0])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)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; + + //printf("V0 TMath::Exp(2*scalempt) %e, downscaleF %e \n",TMath::Exp(2*scalempt), downscaleF); + if (TMath::Exp(2*scalempt)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; + mother = stack->Particle(motherLabel); + +return mother; +} + +//_____________________________________________________________________________ +Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(const Int_t label, AliStack *const stack) +{ + Bool_t isFromConversion = kFALSE; + + if(stack) { + TParticle* particle = stack->Particle(label); + + if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) + { + Int_t mech = particle->GetUniqueID(); // production mechanism + Bool_t isPrim = stack->IsPhysicalPrimary(label); + + Int_t motherLabel = TMath::Abs(particle->GetMother(0)); + TParticle* mother = stack->Particle(motherLabel); + if(mother) { + Int_t motherPdg = mother->GetPdgCode(); + + if(!isPrim && mech==5 && motherPdg==kGamma) { + isFromConversion=kTRUE; + } + } + } + } + + return isFromConversion; +} + +//_____________________________________________________________________________ +Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(const Int_t label, AliStack *const stack) +{ + Bool_t isFromMaterial = kFALSE; + + if(stack) { + TParticle* particle = stack->Particle(label); + + if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) + { + Int_t mech = particle->GetUniqueID(); // production mechanism + Bool_t isPrim = stack->IsPhysicalPrimary(label); + + Int_t motherLabel = TMath::Abs(particle->GetMother(0)); + TParticle* mother = stack->Particle(motherLabel); + if(mother) { + if(!isPrim && mech==13) { + isFromMaterial=kTRUE; + } + } + } + } + + return isFromMaterial; +} + +//_____________________________________________________________________________ +Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(const Int_t label, AliStack *const stack) +{ + Bool_t isFromStrangeness = kFALSE; + + if(stack) { + TParticle* particle = stack->Particle(label); + + if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) + { + Int_t mech = particle->GetUniqueID(); // production mechanism + Bool_t isPrim = stack->IsPhysicalPrimary(label); + + Int_t motherLabel = TMath::Abs(particle->GetMother(0)); + TParticle* mother = stack->Particle(motherLabel); + if(mother) { + Int_t motherPdg = mother->GetPdgCode(); + + // K+-, lambda, antilambda, K0s decays + if(!isPrim && mech==4 && + (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short)) + { + isFromStrangeness = kTRUE; + } + } + } + } + + return isFromStrangeness; +} + + +//_____________________________________________________________________________ +void AliAnalysisTaskFilteredTree::FinishTaskOutput() +{ + // + // Called one at the end + // locally on working node + // + + // must be deleted to store trees + if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0; + + // open temporary file and copy trees to the ouptut container + + TChain* chain = 0; + // + chain = new TChain("highPt"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fHighPtTree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fHighPtTree) fHighPtTree->Print(); + + // + chain = new TChain("V0s"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fV0Tree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fV0Tree) fV0Tree->Print(); + + // + chain = new TChain("dEdx"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fdEdxTree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fdEdxTree) fdEdxTree->Print(); + + // + chain = new TChain("Laser"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fLaserTree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fLaserTree) fLaserTree->Print(); + + // + chain = new TChain("MCEffTree"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fMCEffTree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fMCEffTree) fMCEffTree->Print(); + + // + chain = new TChain("CosmicPairs"); + if(chain) { + chain->Add("jotwinow_Temp_Trees.root"); + fCosmicPairsTree = chain->CopyTree("1"); + delete chain; chain=0; + } + if(fCosmicPairsTree) fCosmicPairsTree->Print(); + + + OpenFile(1); + + // Post output data. + PostData(1, fHighPtTree); + PostData(2, fV0Tree); + PostData(3, fdEdxTree); + PostData(4, fLaserTree); + PostData(5, fMCEffTree); + PostData(6, fCosmicPairsTree); +} + +//_____________________________________________________________________________ +void AliAnalysisTaskFilteredTree::Terminate(Option_t *) +{ + // Called one at the end + /* + fOutputSummary = dynamic_cast (GetOutputData(1)); + if(fOutputSummary) delete fOutputSummary; fOutputSummary=0; + TChain* chain = new TChain("highPt"); + if(!chain) return; + chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root"); + TTree *tree = chain->CopyTree("1"); + if (chain) { delete chain; chain=0; } + if(!tree) return; + tree->Print(); + fOutputSummary = tree; + + if (!fOutputSummary) { + Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1)); + return; + } + */ +} + +//_____________________________________________________________________________ +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; +} + + diff --git a/PWGPP/AliAnalysisTaskFilteredTree.h b/PWGPP/AliAnalysisTaskFilteredTree.h new file mode 100644 index 00000000000..bd6c5b15ce2 --- /dev/null +++ b/PWGPP/AliAnalysisTaskFilteredTree.h @@ -0,0 +1,139 @@ +#ifndef ALIDNDPTTRACKDUMPTASK_H +#define ALIDNDPTTRACKDUMPTASK_H + +//------------------------------------------------------------------------------ +// Task to dump track information +// TPC constrained and TPC+ITS combined +// for outliers analysis. +// +// Author: J.Otwinowski 19/06/2011 +//------------------------------------------------------------------------------ + +class AliESDEvent; +class AliMCEvent; +class AliFilteredTreeEventCuts; +class AliFilteredTreeAcceptanceCuts; +class AliESDtrackCuts; +class AliMagFMaps; +class AliESDEvent; +class AliMCEvent; +class AliKFParticle; +class AliESDv0; +class TList; +class TObjArray; +class TTree; +class TTreeSRedirector; + +#include "AliTriggerAnalysis.h" +#include "AliAnalysisTaskSE.h" + +class AliAnalysisTaskFilteredTree : public AliAnalysisTaskSE { + public: + + enum EAnalysisMode { kInvalidAnalysisMode=-1, + kTPCITSAnalysisMode=0, + kTPCAnalysisMode=1 }; + + AliAnalysisTaskFilteredTree(const char *name = "AliAnalysisTaskFilteredTree"); + virtual ~AliAnalysisTaskFilteredTree(); + + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *); + virtual Bool_t Notify(); + virtual void FinishTaskOutput(); + void SetUseMCInfo(Bool_t info) { fUseMCInfo = info; } + Bool_t IsUseMCInfo() const { return fUseMCInfo; } + void SetUseESDfriends(Bool_t friends) { fUseESDfriends = friends; } + Bool_t IsUseESDfriends() const { return fUseESDfriends; } + + // Process events + void ProcessAll(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void Process(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void ProcessV0(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void ProcessdEdx(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void ProcessLaser(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void ProcessMCEff(AliESDEvent *const esdEvent=0, AliMCEvent *const mcEvent=0, AliESDfriend *const esdFriend=0); + void ProcessCosmics(AliESDEvent *const esdEvent=0); + + void SetEventCuts(AliFilteredTreeEventCuts* const cuts) { fFilteredTreeEventCuts = cuts; } + void SetAcceptanceCuts(AliFilteredTreeAcceptanceCuts* const cuts) { fFilteredTreeAcceptanceCuts = cuts; } + void SetRecAcceptanceCuts(AliFilteredTreeAcceptanceCuts* const cuts) { fFilteredTreeRecAcceptanceCuts = cuts; } + void SetTrackCuts(AliESDtrackCuts* const cuts) { fEsdTrackCuts = cuts; } + void SetTrigger(const AliTriggerAnalysis::Trigger trigger) { fTrigger = trigger; } + void SetAnalysisMode(const EAnalysisMode mode) { fAnalysisMode = mode; } + + AliFilteredTreeEventCuts* GetEventCuts() const { return fFilteredTreeEventCuts; } + AliFilteredTreeAcceptanceCuts* GetAcceptanceCuts() const { return fFilteredTreeAcceptanceCuts; } + AliFilteredTreeAcceptanceCuts* GetRecAcceptanceCuts() const { return fFilteredTreeRecAcceptanceCuts; } + AliESDtrackCuts* GetTrackCuts() const { return fEsdTrackCuts; } + AliTriggerAnalysis::Trigger GetTrigger() const { return fTrigger; } + EAnalysisMode GetAnalysisMode() const { return fAnalysisMode; } + + TString GetCentralityEstimator() const {return fCentralityEstimator; } + void SetCentralityEstimator(TString centEst="V0M") { fCentralityEstimator = centEst; } + + Bool_t IsFromConversion(const Int_t label, AliStack *const stack); + Bool_t IsFromMaterial(const Int_t label, AliStack *const stack); + Bool_t IsFromStrangeness(const Int_t label, AliStack *const stack); + TParticle *GetMother(TParticle *const particle, AliStack *const stack); + + Bool_t ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3]); + Bool_t ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3]); + + // v0s selection + Int_t GetKFParticle(AliESDv0 *const v0, AliESDEvent * const event, AliKFParticle & kfparticle); + Bool_t IsV0Downscaled(AliESDv0 *const v0); + Bool_t IsHighDeDxParticle(AliESDtrack * const track); + + void SetLowPtTrackDownscaligF(Double_t fact) { fLowPtTrackDownscaligF = fact; } + void SetLowPtV0DownscaligF(Double_t fact) { fLowPtV0DownscaligF = fact; } + + void SetProcessCosmics(Bool_t flag) { fProcessCosmics = flag; } + Bool_t GetProcessCosmics() { return fProcessCosmics; } + + void SetProcessAll(Bool_t proc) { fProcessAll = proc; } + static Int_t GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts); + + private: + + AliESDEvent *fESD; //! ESD event + AliMCEvent *fMC; //! MC event + AliESDfriend *fESDfriend; //! ESDfriend event + TList* fOutput; //! list send on output slot 0 + TIterator *fPitList; //! iterator over the output objetcs + + Bool_t fUseMCInfo; // use MC information + Bool_t fUseESDfriends; // use esd friends + + AliFilteredTreeEventCuts *fFilteredTreeEventCuts; // event cuts + AliFilteredTreeAcceptanceCuts *fFilteredTreeAcceptanceCuts; // acceptance cuts + AliFilteredTreeAcceptanceCuts *fFilteredTreeRecAcceptanceCuts; // additional recontruction acceptance cuts (not used for MC truth) + AliESDtrackCuts *fEsdTrackCuts; // esd track cuts + AliTriggerAnalysis::Trigger fTrigger; // trigger settings + EAnalysisMode fAnalysisMode; // analysis mode TPC only, TPC + ITS + + TTreeSRedirector* fTreeSRedirector; //! temp tree to dump output + + TString fCentralityEstimator; // use centrality can be "VOM" (default), "FMD", "TRK", "TKL", "CL0", "CL1", "V0MvsFMD", "TKLvsV0M", "ZEMvsZDC" + + Double_t fLowPtTrackDownscaligF; // low pT track downscaling factor + Double_t fLowPtV0DownscaligF; // low pT V0 downscaling factor + Double_t fProcessAll; // Calculate all track properties including MC + + Bool_t fProcessCosmics; // look for cosmic pairs from random trigger + + TTree* fHighPtTree; //! list send on output slot 0 + TTree* fV0Tree; //! list send on output slot 0 + TTree* fdEdxTree; //! list send on output slot 0 + TTree* fLaserTree; //! list send on output slot 0 + TTree* fMCEffTree; //! list send on output slot 0 + TTree* fCosmicPairsTree; //! list send on output slot 0 + + AliAnalysisTaskFilteredTree(const AliAnalysisTaskFilteredTree&); // not implemented + AliAnalysisTaskFilteredTree& operator=(const AliAnalysisTaskFilteredTree&); // not implemented + + ClassDef(AliAnalysisTaskFilteredTree, 1); // example of analysis +}; + +#endif diff --git a/PWGPP/AliFilteredTreeAcceptanceCuts.cxx b/PWGPP/AliFilteredTreeAcceptanceCuts.cxx new file mode 100644 index 00000000000..b4f70bf52a1 --- /dev/null +++ b/PWGPP/AliFilteredTreeAcceptanceCuts.cxx @@ -0,0 +1,177 @@ +/************************************************************************** +* 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. * +**************************************************************************/ +// last change: 2011-04-04 by M.Knichel + +#include +#include + +#include "AliLog.h" +#include "AliESDtrack.h" +#include "AliExternalTrackParam.h" +#include "TParticle.h" + +#include "AliFilteredTreeAcceptanceCuts.h" + +using namespace std; + +ClassImp(AliFilteredTreeAcceptanceCuts) + +//_____________________________________________________________________________ +AliFilteredTreeAcceptanceCuts::AliFilteredTreeAcceptanceCuts(const Char_t* name,const Char_t *title) : +AliAnalysisCuts(name, title) +, fMinEta(0) +, fMaxEta(0) +, fMinPhi(0) +, fMaxPhi(0) +, fMinPt(0) +, fMaxPt(0) +, fExcludeMinEta(0) +, fExcludeMaxEta(0) +, fExcludeMinPhi(0) +, fExcludeMaxPhi(0) +, fExcludeMinEta2(0) +, fExcludeMaxEta2(0) +, fExcludeMinPhi2(0) +, fExcludeMaxPhi2(0) +, fCheckRange(kFALSE) +, fMaxDCAr(0) +, fMaxDCAz(0) +{ + // default constructor + + // init data members with defaults + Init(); +} + +//_____________________________________________________________________________ +AliFilteredTreeAcceptanceCuts::~AliFilteredTreeAcceptanceCuts() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliFilteredTreeAcceptanceCuts::Init() +{ + // set default values + SetEtaRange(); + SetPhiRange(); + SetPtRange(); + SetMaxDCAr(); + SetMaxDCAz(); +} + +//_____________________________________________________________________________ +Bool_t AliFilteredTreeAcceptanceCuts::AcceptTrack(AliESDtrack *track) +{ + // check acceptance cuts for AliESDtrack + if(!track) return kFALSE; + + Float_t eta = track->Eta(); + Float_t phi = track->Phi(); + Float_t pt = track->Pt(); + + if(eta < fMinEta) return kFALSE; + if(eta > fMaxEta) return kFALSE; + if(phi < fMinPhi) return kFALSE; + if(phi > fMaxPhi) return kFALSE; + if(pt < fMinPt) return kFALSE; + if(pt > fMaxPt) return kFALSE; + +return kTRUE; +} + +Bool_t AliFilteredTreeAcceptanceCuts::AcceptTrackLocalTPC(AliESDtrack *track) +{ + // check acceptance cuts for AliESDtrack + if(!track) return kFALSE; + const AliExternalTrackParam *innerParam = track->GetInnerParam(); + if(!innerParam) return kFALSE; + + Float_t eta = track->Eta(); + Float_t phi = TMath::ATan2(innerParam->Py(),innerParam->Px()); + + if (fCheckRange) { + if ((eta > fExcludeMinEta) && (eta < fExcludeMaxEta) && (phi > fExcludeMinPhi) && (phi < fExcludeMaxPhi)) { return kFALSE; } + if ((eta > fExcludeMinEta2) && (eta < fExcludeMaxEta2) && (phi > fExcludeMinPhi2) && (phi < fExcludeMaxPhi2)) { return kFALSE; } + } + +return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliFilteredTreeAcceptanceCuts::AcceptTrack(AliExternalTrackParam *track) +{ + // check acceptance cuts for AliESDtrack + if(!track) return kFALSE; + + Float_t eta = track->Eta(); + Float_t phi = track->Phi(); + Float_t pt = track->Pt(); + + if(eta < fMinEta) return kFALSE; + if(eta > fMaxEta) return kFALSE; + if(phi < fMinPhi) return kFALSE; + if(phi > fMaxPhi) return kFALSE; + if(pt < fMinPt) return kFALSE; + if(pt > fMaxPt) return kFALSE; + +return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliFilteredTreeAcceptanceCuts::AcceptTrack(TParticle *particle) +{ + // check acceptance cuts for TParticle + if(!particle) return kFALSE; + + Float_t eta = particle->Eta(); + Float_t phi = particle->Phi(); + Float_t pt = particle->Pt(); + + if(eta < fMinEta) return kFALSE; + if(eta > fMaxEta) return kFALSE; + if(phi < fMinPhi) return kFALSE; + if(phi > fMaxPhi) return kFALSE; + if(pt < fMinPt) return kFALSE; + if(pt > fMaxPt) return kFALSE; + +return kTRUE; +} + +//_____________________________________________________________________________ +Long64_t AliFilteredTreeAcceptanceCuts::Merge(TCollection* list) +{ + // Merge list of objects (needed by PROOF) + if (!list) + return 0; + + if (list->IsEmpty()) + return 1; + + TIterator* iter = list->MakeIterator(); + TObject* obj = 0; + + Int_t count=0; + while((obj = iter->Next()) != 0) + { + AliFilteredTreeAcceptanceCuts* entry = dynamic_cast(obj); + if (entry == 0) + continue; + + count++; + } + +return count; +} diff --git a/PWGPP/AliFilteredTreeAcceptanceCuts.h b/PWGPP/AliFilteredTreeAcceptanceCuts.h new file mode 100644 index 00000000000..87c27308d1e --- /dev/null +++ b/PWGPP/AliFilteredTreeAcceptanceCuts.h @@ -0,0 +1,97 @@ +#ifndef ALIFILTEREDTREEACCEPTANCECUTS_H +#define ALIFILTEREDTREEACCEPTANCECUTS_H + +//------------------------------------------------------------------------------ +// Class to keep selection cuts for MC tracks. +// +// Author: J.Otwinowski 03/11/2008 +// last change: 2011-04-04 by M.Knichel +//------------------------------------------------------------------------------ + +class TParticle; +class AliESDtrack; +class AliExternalTrackParam; + +#include "AliAnalysisCuts.h" + +class AliFilteredTreeAcceptanceCuts : public AliAnalysisCuts +{ +public: + AliFilteredTreeAcceptanceCuts(const Char_t* name ="AliFilteredTreeAcceptanceCuts", const Char_t *title =""); + virtual ~AliFilteredTreeAcceptanceCuts(); + + // setters + void SetEtaRange(const Float_t min=-1e99, const Float_t max=1e99) { fMinEta=min; fMaxEta=max; } + void SetPhiRange(const Float_t min=-1e99, const Float_t max=1e99) { fMinPhi=min; fMaxPhi=max;} + void SetPtRange(const Float_t min=-1e99, const Float_t max=1e99) { fMinPt=min; fMaxPt=max;} + void SetExcludeEtaPhiRange(const Float_t etaMin, const Float_t etaMax, const Float_t phiMin, const Float_t phiMax) + { fExcludeMinEta = etaMin; fExcludeMaxEta = etaMax; fExcludeMinPhi = phiMin; fExcludeMaxPhi = phiMax; fCheckRange=kTRUE; } + void SetExcludeEtaPhiRange2(const Float_t etaMin, const Float_t etaMax, const Float_t phiMin, const Float_t phiMax) + { fExcludeMinEta2 = etaMin; fExcludeMaxEta2 = etaMax; fExcludeMinPhi2 = phiMin; fExcludeMaxPhi2 = phiMax; fCheckRange=kTRUE; } + + void SetMaxDCAr(const Float_t max=1e99) { fMaxDCAr=max;} + void SetMaxDCAz(const Float_t max=1e99) { fMaxDCAz=max;} + + // getters + Float_t GetMinEta() const {return fMinEta;} + Float_t GetMaxEta() const {return fMaxEta;} + Float_t GetMinPhi() const {return fMinPhi;} + Float_t GetMaxPhi() const {return fMaxPhi;} + Float_t GetMinPt() const {return fMinPt;} + Float_t GetMaxPt() const {return fMaxPt;} + + Bool_t GetCheckRange() const { return fCheckRange; } + Float_t GetExcludeMinEta() const { return fExcludeMinEta; } + Float_t GetExcludeMaxEta() const { return fExcludeMaxEta; } + Float_t GetExcludeMinPhi() const { return fExcludeMinPhi; } + Float_t GetExcludeMaxPhi() const { return fExcludeMaxPhi; } + + Float_t GetMaxDCAr() const {return fMaxDCAr;} + Float_t GetMaxDCAz() const {return fMaxDCAz;} + + // cuts init function + void Init(); + + // check MC tracks + virtual Bool_t IsSelected(TObject *) {return kTRUE;} + virtual Bool_t IsSelected(TList *) {return kTRUE;} + + // + Bool_t AcceptTrack(AliESDtrack *track); + Bool_t AcceptTrackLocalTPC(AliESDtrack *track); + Bool_t AcceptTrack(AliExternalTrackParam *track); + Bool_t AcceptTrack(TParticle *particle); + + // Merge output objects (needed by PROOF) + virtual Long64_t Merge(TCollection* list); + +private: + Float_t fMinEta; // min pseudorapidity + Float_t fMaxEta; // max pseudorapidity + Float_t fMinPhi; // min azimuthal angle (rad) + Float_t fMaxPhi; // max azimuthal angle (rad) + Float_t fMinPt; // min pt + Float_t fMaxPt; // max pt + + Float_t fExcludeMinEta; + Float_t fExcludeMaxEta; + Float_t fExcludeMinPhi; + Float_t fExcludeMaxPhi; + Float_t fExcludeMinEta2; + Float_t fExcludeMaxEta2; + Float_t fExcludeMinPhi2; + Float_t fExcludeMaxPhi2; + Bool_t fCheckRange; + + // max DCAr and DCAz with respect + // to nominal vertex position + Float_t fMaxDCAr; // min DCAr + Float_t fMaxDCAz; // max DCAz + + AliFilteredTreeAcceptanceCuts(const AliFilteredTreeAcceptanceCuts&); // not implemented + AliFilteredTreeAcceptanceCuts& operator=(const AliFilteredTreeAcceptanceCuts&); // not implemented + + ClassDef(AliFilteredTreeAcceptanceCuts, 1) +}; + +#endif // diff --git a/PWGPP/AliFilteredTreeEventCuts.cxx b/PWGPP/AliFilteredTreeEventCuts.cxx new file mode 100644 index 00000000000..04c8caf9ec6 --- /dev/null +++ b/PWGPP/AliFilteredTreeEventCuts.cxx @@ -0,0 +1,475 @@ +/************************************************************************** +* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +* * +* Author: The ALICE Off-line Project. * +* Contributors are mentioned in the code where appropriate. * +* * +* Permission to use, copy, modify and distribute this software and its * +* documentation strictly for non-commercial purposes is hereby granted * +* without fee, provided that the above copyright notice appears in all * +* copies and that both the copyright notice and this permission notice * +* appear in the supporting documentation. The authors make no claims * +* about the suitability of this software for any purpose. It is * +* provided "as is" without express or implied warranty. * +**************************************************************************/ + +#include +#include + +#include "AliLog.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliMCEvent.h" +#include "AliHeader.h" +#include "AliGenEventHeader.h" +#include "AliGenPythiaEventHeader.h" +#include "AliGenCocktailEventHeader.h" +#include "AliGenDPMjetEventHeader.h" +#include "AliStack.h" + + +#include "AliFilteredTreeEventCuts.h" + +using namespace std; + +ClassImp(AliFilteredTreeEventCuts) + +Int_t AliFilteredTreeEventCuts::fgLastProcessType = -1; + +//_____________________________________________________________________________ +AliFilteredTreeEventCuts::AliFilteredTreeEventCuts(const Char_t* name,const Char_t *title) : +AliAnalysisCuts(name, title) +, fTriggerRequired(kTRUE) +, fRecVertexRequired(kTRUE) +, fEventProcessType(kInvalidProcess) +, fMinNContributors(0) +, fMaxNContributors(0) +, fMaxR(0) +, fMinZv(0) +, fMaxZv(0) +, fMeanXv(0) +, fMeanYv(0) +, fMeanZv(0) +, fSigmaMeanXv(0) +, fSigmaMeanYv(0) +, fSigmaMeanZv(0) +, fRedoTPCVertex(kTRUE) +, fUseBeamSpotConstraint(kTRUE) +, fEventSelectedRequired(kFALSE) +{ + // default constructor + + // init data members with defaults + Init(); +} + +//_____________________________________________________________________________ +AliFilteredTreeEventCuts::~AliFilteredTreeEventCuts() +{ + // destructor +} + +//_____________________________________________________________________________ +void AliFilteredTreeEventCuts::Init() +{ + // set default values + SetTriggerRequired(); + SetRecVertexRequired(); + SetEventProcessType(); + SetNContributorsRange(); + SetMaxR(); + SetZvRange(); + SetMeanXYZv(); + SetSigmaMeanXYZv(); + SetRedoTPCVertex(); + SetUseBeamSpotConstraint(); +} + +//_____________________________________________________________________________ +Bool_t AliFilteredTreeEventCuts::AcceptEvent(AliESDEvent *esdEvent,AliMCEvent *mcEvent, const AliESDVertex *vtx) +{ + // Check event selection cuts + Bool_t retValue=kTRUE; + + if(!esdEvent) return kFALSE; + if(!IsRecVertexRequired()) return kTRUE; + if(!vtx) return kFALSE; + if(!vtx->GetStatus()) return kFALSE; + + if(mcEvent) { + // check MC event conditions + AliHeader* header = mcEvent->Header(); + if(!header) return kFALSE; + + // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive) + if(fEventProcessType == kInvalidProcess) { + retValue=kTRUE; + } + else if(fEventProcessType == kSD || fEventProcessType == kDD) { + MCProcessType processType = GetEventProcessType(header); + if(processType == kND) retValue=kFALSE; + else retValue=kTRUE; + } + else if(fEventProcessType == GetEventProcessType(header)) { + retValue=kTRUE; + } + else + retValue=kFALSE; + } + + if(vtx->GetZv() < fMinZv) return kFALSE; + if(vtx->GetZv() > fMaxZv) return kFALSE; + +return retValue; +} + +//_____________________________________________________________________________ +Bool_t AliFilteredTreeEventCuts::AcceptMCEvent(AliMCEvent *mcEvent) +{ + // Check event selection cuts + if(!mcEvent) return kFALSE; + + Bool_t retValue=kTRUE; + + // check MC event conditions + AliHeader* header = mcEvent->Header(); + if(!header) return kFALSE; + + AliGenEventHeader* genHeader = header->GenEventHeader(); + if (!genHeader) { + AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); + return kFALSE; + } + TArrayF vtxMC(3); + genHeader->PrimaryVertex(vtxMC); + + // select event type (ND-non diffractive, SD-single diffractive, DD-double diffractive) + if(fEventProcessType == kInvalidProcess) { + retValue=kTRUE; + } else { + if(fEventProcessType == GetEventProcessType(header)) retValue=kTRUE; + else retValue=kFALSE; + } + + /* + Float_t R = TMath::Sqrt(vtxMC[0]*vtxMC[0]+vtxMC[1]*vtxMC[1]); + if(R > fMaxR) return kFALSE; + */ + + if(vtxMC[2] < fMinZv) return kFALSE; + if(vtxMC[2] > fMaxZv) return kFALSE; + +return retValue; +} + +//_____________________________________________________________________________ +Long64_t AliFilteredTreeEventCuts::Merge(TCollection* list) +{ + // Merge list of objects (needed by PROOF) + if (!list) + return 0; + + if (list->IsEmpty()) + return 1; + + TIterator* iter = list->MakeIterator(); + TObject* obj = 0; + + Int_t count=0; + while((obj = iter->Next()) != 0) + { + AliFilteredTreeEventCuts* entry = dynamic_cast(obj); + if (entry == 0) + continue; + + count++; + } + +return count; +} + +//_____________________________________________________________________________ +AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) { + // + // get the process type of the event. + // + + + // Check for simple headers first + + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(aHeader->GenEventHeader()); + if (pythiaGenHeader) { + return GetPythiaEventProcessType(pythiaGenHeader,adebug); + } + + AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast(aHeader->GenEventHeader()); + if (dpmJetGenHeader) { + return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); + } + + + // check for cocktail + + AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast(aHeader->GenEventHeader()); + if (!genCocktailHeader) { + printf("AliFilteredTreeEventCuts::GetProcessType : Unknown header type (not Pythia or Cocktail). \n"); + return kInvalidProcess; + } + + TList* headerList = genCocktailHeader->GetHeaders(); + if (!headerList) { + return kInvalidProcess; + } + + for (Int_t i=0; iGetEntries(); i++) { + + pythiaGenHeader = dynamic_cast(headerList->At(i)); + if (pythiaGenHeader) { + return GetPythiaEventProcessType(pythiaGenHeader,adebug); + } + + dpmJetGenHeader = dynamic_cast(headerList->At(i)); + if (dpmJetGenHeader) { + return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); + } + } + return kInvalidProcess; +} + +//____________________________________________________________________ +AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment) +{ + // + // get process type + // + // diffTreatment: + // kMCFlags: use MC flags + // kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND + // kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND + // kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND + // + + MCProcessType mcProcessType = GetEventProcessType(header); + + if (diffTreatment == kMCFlags) + return mcProcessType; + + if (!esd) + { + Printf("ERROR: AliFilteredTreeEventCuts::GetEventProcessType: diffTreatment != kMCFlags and esd == 0"); + return kInvalidProcess; + } + + Float_t cms = esd->GetESDRun()->GetBeamEnergy(); + if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV()) + cms *= 2; + //Printf("cms = %f", cms); + + if (diffTreatment == kUA5Cuts && mcProcessType == kSD) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) + return kSD; + } + else if (diffTreatment == kE710Cuts && mcProcessType == kSD) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05)) + return kSD; + } + else if (diffTreatment == kALICEHadronLevel) + { + if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) + return kSD; + } + + if (mcProcessType == kSD) + return kND; + + return mcProcessType; +} + + +AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { + + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(aHeader); + + if (!pythiaGenHeader) { + printf("AliFilteredTreeEventCuts::GetProcessType : Unknown gen Header type). \n"); + return kInvalidProcess; + } + + + Int_t pythiaType = pythiaGenHeader->ProcessType(); + fgLastProcessType = pythiaType; + MCProcessType globalType = kInvalidProcess; + + + if (adebug) { + printf("AliFilteredTreeEventCuts::GetProcessType : Pythia process type found: %d \n",pythiaType); + } + + + if(pythiaType==92||pythiaType==93){ + globalType = kSD; + } + else if(pythiaType==94){ + globalType = kDD; + } + //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??} + else { + globalType = kND; + } + return globalType; +} + + +AliFilteredTreeEventCuts::MCProcessType AliFilteredTreeEventCuts::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { + // + // get the process type of the event. + // + + // can only read pythia headers, either directly or from cocktalil header + AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast(aHeader); + + if (!dpmJetGenHeader) { + printf("AliFilteredTreeEventCuts::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n"); + return kInvalidProcess; + } + + Int_t dpmJetType = dpmJetGenHeader->ProcessType(); + fgLastProcessType = dpmJetType; + MCProcessType globalType = kInvalidProcess; + + + if (adebug) { + printf("AliFilteredTreeEventCuts::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType); + } + + + if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction + globalType = kND; + } + else if (dpmJetType==5 || dpmJetType==6) { + globalType = kSD; + } + else if (dpmJetType==7) { + globalType = kDD; + } + return globalType; +} + +//____________________________________________________________________ +Bool_t AliFilteredTreeEventCuts::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax) +{ + // + // return if process is single diffractive on hadron level + // + // xiMax and xiMin cut on M^2/s + // + // Based on code from Martin Poghoysan + // + + TParticle* part1 = 0; + TParticle* part2 = 0; + + Double_t smallestY = 1e10; + Double_t largestY = -1e10; + + for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++) + { + TParticle* part = stack->Particle(iParticle); + if (!part) + continue; + + Int_t pdg = TMath::Abs(part->GetPdgCode()); + + Int_t child1 = part->GetFirstDaughter(); + Int_t ist = part->GetStatusCode(); + + Int_t mfl = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg)))); + if (child1 > -1 || ist != 1) + mfl = 0; // select final state charm and beauty + if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4)) + continue; + Int_t imother = part->GetFirstMother(); + if (imother>0) + { + TParticle *partM = stack->Particle(imother); + Int_t pdgM=TMath::Abs(partM->GetPdgCode()); + if (pdgM==111 || pdgM==3124 || pdgM==3212) + continue; + } + + Double_t y = 0; + + // fix for problem with getting mass of particle 3124 + if (pdg != 3124) + y = Rapidity(part->Pt(), part->Pz(), part->GetMass()); + else + y = Rapidity(part->Pt(), part->Pz(), 1.5195); + + if (y < smallestY) + { + smallestY = y; + part1 = part; + } + + if (y > largestY) + { + largestY = y; + part2 = part; + } + } + + if (part1 == 0 || part2 == 0) + return kFALSE; + + Int_t pdg1 = part1->GetPdgCode(); + Int_t pdg2 = part2->GetPdgCode(); + + Double_t pt1 = part1->Pt(); + Double_t pt2 = part2->Pt(); + Double_t pz1 = part1->Pz(); + Double_t pz2 = part2->Pz(); + + Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938)); + Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938)); + + Int_t arm = -99999; + if (pdg1 == 2212 && pdg2 == 2212) + { + if (y1 > y2) + arm = 0; + else + arm = 1; + } + else if (pdg1 == 2212) + arm = 0; + else if (pdg2 == 2212) + arm = 1; + + Double_t M02s = 1. - 2 * part1->Energy() / cms; + Double_t M12s = 1. - 2 * part2->Energy() / cms; + + if (arm == 0 && M02s > xiMin && M02s < xiMax) + return kTRUE; + else if (arm == 1 && M12s > xiMin && M12s < xiMax) + return kTRUE; + + return kFALSE; +} + +//____________________________________________________________________ +Double_t AliFilteredTreeEventCuts::Rapidity(Double_t pt, Double_t pz, Double_t m) +{ + // + // calculates rapidity keeping the sign in case E == pz + // + + Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m); + if (energy != TMath::Abs(pz)) + return 0.5*TMath::Log((energy+pz)/(energy-pz)); + + Printf("W- mt=0"); + return TMath::Sign(1.e30,pz); +} + diff --git a/PWGPP/AliFilteredTreeEventCuts.h b/PWGPP/AliFilteredTreeEventCuts.h new file mode 100644 index 00000000000..a5a97a4ad45 --- /dev/null +++ b/PWGPP/AliFilteredTreeEventCuts.h @@ -0,0 +1,130 @@ +#ifndef ALIFILTEREDTREEEVENTCUTS_H +#define ALIFILTEREDTREEEVENTCUTS_H + +//------------------------------------------------------------------------------ +// Class to keep event selection cuts for dNdPt analysis. +// +// Author: J.Otwinowski 01/11/2008 +//------------------------------------------------------------------------------ + +#include "AliAnalysisCuts.h" + +class AliESDEvent; +class AliESDVertex; +class AliMCEvent; +class AliHeader; +class AliGenEventHeader; +class AliStack; + +class AliFilteredTreeEventCuts : public AliAnalysisCuts +{ +public: + AliFilteredTreeEventCuts(const Char_t* name ="AliFilteredTreeEventCuts", const Char_t *title =""); + virtual ~AliFilteredTreeEventCuts(); + + //TODO: copied from AliPWG0Helper, find a more central place for these + enum AnalysisMode { kInvalid = -1, kSPD = 0x1, kTPC = 0x2, kTPCITS = 0x4, kFieldOn = 0x8, kSPDOnlyL0 = 0x10, kTPCSPD = 0x20}; + enum MCProcessType { kInvalidProcess = -1, kND = 0x1, kDD = 0x2, kSD = 0x4, kOnePart = 0x8 }; + enum DiffTreatment { kMCFlags = 0, kUA5Cuts = 1, kE710Cuts, kALICEHadronLevel }; + + // setters + void SetTriggerRequired(const Bool_t bFlag=kTRUE) {fTriggerRequired=bFlag;} + void SetRecVertexRequired(const Bool_t bFlag=kTRUE) {fRecVertexRequired=bFlag;} + void SetEventProcessType(MCProcessType type=kInvalidProcess) {fEventProcessType=type;} + void SetNContributorsRange(const Float_t min=0.,const Float_t max=1e99) {fMinNContributors=min; fMaxNContributors=max;} + void SetMaxR(const Float_t max=1e99) {fMaxR=max;} + void SetZvRange(const Float_t min=-1e99, const Float_t max=1e99) {fMinZv=min; fMaxZv=max;} + + void SetMeanXYZv(const Float_t xv=0.0, const Float_t yv=0.0, const Float_t zv=0.0) { + fMeanXv = xv; fMeanYv = yv; fMeanZv = zv; + } + + void SetSigmaMeanXYZv(const Float_t sxv=1.0, const Float_t syv=1.0, const Float_t szv=10.0) { + fSigmaMeanXv = sxv; fSigmaMeanYv = syv; fSigmaMeanZv = szv; + } + + + void SetRedoTPCVertex(const Bool_t redo = kTRUE) {fRedoTPCVertex = redo;} + void SetUseBeamSpotConstraint(const Bool_t useConstr = kTRUE) {fUseBeamSpotConstraint = useConstr;} + void SetEventSelectedRequired(const Bool_t evtSel = kTRUE) {fEventSelectedRequired = evtSel;} + + + // getters + Bool_t IsEventSelectedRequired() const {return fEventSelectedRequired;} + Bool_t IsTriggerRequired() const {return fTriggerRequired;} + Bool_t IsRecVertexRequired() const {return fRecVertexRequired;} + Int_t GetEventProcessType() const {return fEventProcessType;} + Float_t GetMinNContributors() const {return fMinNContributors;} + Float_t GetMaxNContributors() const {return fMaxNContributors;} + Float_t GetMaxR() const {return fMaxR;} + Float_t GetMinZv() const {return fMinZv;} + Float_t GetMaxZv() const {return fMaxZv;} + + Float_t GetMeanXv() const {return fMeanXv;} + Float_t GetMeanYv() const {return fMeanYv;} + Float_t GetMeanZv() const {return fMeanZv;} + + Float_t GetSigmaMeanXv() const {return fSigmaMeanXv;} + Float_t GetSigmaMeanYv() const {return fSigmaMeanYv;} + Float_t GetSigmaMeanZv() const {return fSigmaMeanZv;} + + Bool_t IsRedoTPCVertex() const {return fRedoTPCVertex;} + Bool_t IsUseBeamSpotConstraint() const {return fUseBeamSpotConstraint;} + + + // cuts init function + void Init(); + + // check MC tracks + Bool_t IsSelected(TObject *) {return kTRUE;} + Bool_t IsSelected(TList *) {return kTRUE;} + + // accept event + Bool_t AcceptEvent(AliESDEvent *event=0, AliMCEvent *mcEvent=0, const AliESDVertex *vtx=0); + Bool_t AcceptMCEvent(AliMCEvent *mcEvent=0); + + // Merge output objects (needed by PROOF) + virtual Long64_t Merge(TCollection* list); + + //statics copied from AliPWG0Helper + static MCProcessType GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, DiffTreatment diffTreatment); + static MCProcessType GetEventProcessType(AliHeader* aHeader, Bool_t adebug = kFALSE); + static MCProcessType GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE); + static MCProcessType GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug = kFALSE); + static Bool_t IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax); + static Double_t Rapidity(Double_t pt, Double_t pz, Double_t m); +protected: + static Int_t fgLastProcessType; // stores the raw value of the last process type extracted + // + +private: + Bool_t fTriggerRequired; // trigger required + Bool_t fRecVertexRequired; // reconstructed event vertex required + Int_t fEventProcessType; // select MC event process type (ND, SD, DD) + Float_t fMinNContributors; // min. number of contributing vertex tracks + Float_t fMaxNContributors; // max. number of contributing vertex tracks + Float_t fMaxR; // max. vertex radii (R = sqrt(Xv^2+Yv^2) + Float_t fMinZv; // min. Zv vertex + Float_t fMaxZv; // max. Zv vertex + + // interaction spot constraint + Float_t fMeanXv; // mean Xv position + Float_t fMeanYv; // mean Yv position + Float_t fMeanZv; // mean Zv position + + Float_t fSigmaMeanXv; // sigma mean Xv position + Float_t fSigmaMeanYv; // sigma mean Yv position + Float_t fSigmaMeanZv; // sigma mean Zv position + + Bool_t fRedoTPCVertex; // redo vertex + Bool_t fUseBeamSpotConstraint; // use beam spot contraints + + Bool_t fEventSelectedRequired; // event with at least one track (pT>0.5 GeV, |eta|<0.8) required + + AliFilteredTreeEventCuts(const AliFilteredTreeEventCuts&); // not implemented + AliFilteredTreeEventCuts& operator=(const AliFilteredTreeEventCuts&); // not implemented + + ClassDef(AliFilteredTreeEventCuts, 1) +}; + +#endif // ALIFILTEREDTREEEVENTCUTS_H diff --git a/PWGPP/PWGPPLinkDef.h b/PWGPP/PWGPPLinkDef.h index a574d4b622a..aecf2d24800 100644 --- a/PWGPP/PWGPPLinkDef.h +++ b/PWGPP/PWGPPLinkDef.h @@ -79,6 +79,10 @@ #pragma link C++ class AliAnaVZEROQA+; #pragma link C++ class AliAnaVZEROPbPb+; +#pragma link C++ class AliAnalysisTaskFilteredTree+; +#pragma link C++ class AliFilteredTreeEventCuts+; +#pragma link C++ class AliFilteredTreeAcceptanceCuts+; + // TRD performance classes #pragma link C++ class AliTenderSupplyTRD+; diff --git a/PWGPP/macros/AddTaskFilteredTree.C b/PWGPP/macros/AddTaskFilteredTree.C new file mode 100644 index 00000000000..fc314eebed5 --- /dev/null +++ b/PWGPP/macros/AddTaskFilteredTree.C @@ -0,0 +1,161 @@ +void AddTaskFilteredTree() +{ + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libTENDER"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWGUDbase"); + gSystem->Load("libTPCcalib"); + gSystem->Load("libPWGPP"); + gSystem->Load("libPWGLFspectra"); + + + gRandom->SetSeed(0); + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + + if (!mgr) { + Error("AddTaskFilteredTree", "No analysis manager found."); + return 0; + } + + // Switch off all AliInfo (too much output!!!) + AliLog::SetGlobalLogLevel(AliLog::kError); + mgr->SetDebugLevel(0); + + + + // + // Create physics trigger selection class + // + //AliPhysicsSelection *physTrigSel = new AliPhysicsSelection(); + + // + // Create event cuts + // + Float_t zvWindow = 30. ; + + AliFilteredTreeEventCuts *evtCuts = new AliFilteredTreeEventCuts("AliFilteredTreeEventCuts","Event cuts"); + evtCuts->SetZvRange(-zvWindow,zvWindow); + evtCuts->SetMeanXYZv(0.0,0.0,0.0); + evtCuts->SetSigmaMeanXYZv(1.0,1.0,10.0); + evtCuts->SetTriggerRequired(kFALSE); + //evtCuts->SetTriggerRequired(kTRUE); + + // + // Create geom. acceptance cuts + // + Float_t etaWindow = 1.0 ; + Float_t ptMin = 0.15 ; + + AliFilteredTreeAcceptanceCuts *accCuts = new AliFilteredTreeAcceptanceCuts("AliFilteredTreeAcceptanceCuts","Geom. acceptance cuts"); + accCuts->SetEtaRange(-etaWindow,etaWindow); + accCuts->SetPtRange(ptMin,1.e10); + accCuts->SetMaxDCAr(3.0); + accCuts->SetMaxDCAz(30.0); + + // + // Create standard esd track cuts + // + + AliESDtrackCuts* esdTrackCuts = CreateCuts(); + if (!esdTrackCuts) { + printf("ERROR: esdTrackCuts could not be created\n"); + return; + } else { + esdTrackCuts->SetHistogramsOn(kTRUE); + } + + Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0); + + // + // Create task + // + AliAnalysisTaskFilteredTree *task = new AliAnalysisTaskFilteredTree("AliAnalysisTaskFilteredTree"); + task->SetUseMCInfo(hasMC); + //task->SetLowPtTrackDownscaligF(1.e4); + //task->SetLowPtV0DownscaligF(1.e2); + task->SetLowPtTrackDownscaligF(1.e7); + task->SetLowPtV0DownscaligF(1.e4); + task->SetProcessAll(kTRUE); + task->SetProcessCosmics(kTRUE); + //task->SetProcessAll(kFALSE); + + // trigger + //task->SelectCollisionCandidates(AliVEvent::kMB); + + // + // set analysis options from the Helper here !!! + // + // AlidNdPtHelper::OutputObject outputObject = AlidNdPtHelper::kCutAnalysisPbPb; + // AlidNdPtHelper::ParticleMode particleMode = AlidNdPtHelper::kAllPart ; + //AlidNdPtHelper::AnalysisMode analysisMode = AlidNdPtHelper::kTPCITS; + + task->SetUseMCInfo(hasMC); + task->SetEventCuts(evtCuts); + task->SetAcceptanceCuts(accCuts); + task->SetTrackCuts(esdTrackCuts); + task->SetAnalysisMode(AliAnalysisTaskFilteredTree::kTPCITSAnalysisMode); + task->SetCentralityEstimator("V0M"); + + // Add task + mgr->AddTask(task); + + // Create containers for input + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + mgr->ConnectInput(task, 0, cinput); + + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("filtered1", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 1, coutput1); + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("filtered2", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 2, coutput2); + AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("filtered3", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 3, coutput3); + AliAnalysisDataContainer *coutput4 = mgr->CreateContainer("filtered4", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 4, coutput4); + AliAnalysisDataContainer *coutput5 = mgr->CreateContainer("filtered5", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 5, coutput5); + AliAnalysisDataContainer *coutput6 = mgr->CreateContainer("filtered6", TTree::Class(), AliAnalysisManager::kOutputContainer, "FilterEvents_Trees.root"); + mgr->ConnectOutput(task, 6, coutput6); +} + +AliESDtrackCuts* CreateCuts(Bool_t fieldOn = kTRUE, Bool_t hists = kTRUE) +{ + AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts"); + if(!esdTrackCuts) return 0; + + if (hists) + esdTrackCuts->DefineHistograms(1); + + Double_t cov1, cov2, cov3, cov4, cov5; + Double_t nSigma; + Double_t maxDCAtoVertex, maxDCAtoVertexXY, maxDCAtoVertexZ; + Double_t minNClustersTPC; + Double_t maxChi2PerClusterTPC; + Double_t minPt, maxPt; + + // + // TPC + // + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + // + // ITS + // + //esdTrackCuts->SetRequireITSRefit(kTRUE); + //esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + // + + TString tag = "TPC+ITS refit and KinkRejection required - for cut studies"; + + // cuts for data without field + if (!fieldOn) + { + cov5 = 1e10; + tag += " without field"; + } + + Printf("Created track cuts for: %s", tag.Data()); + + return esdTrackCuts; +} -- 2.39.3