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