+/**************************************************************************\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
+\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 "TFile.h"\r
+#include "TMatrixD.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 "dNdPt/AlidNdPt.h"\r
+#include "dNdPt/AlidNdPtEventCuts.h"\r
+#include "dNdPt/AlidNdPtAcceptanceCuts.h"\r
+\r
+#include "dNdPt/AlidNdPtTrackDumpTask.h"\r
+\r
+using namespace std;\r
+\r
+ClassImp(AlidNdPtTrackDumpTask)\r
+\r
+//_____________________________________________________________________________\r
+AlidNdPtTrackDumpTask::AlidNdPtTrackDumpTask(const char *name) \r
+ : AliAnalysisTaskSE(name)\r
+ , fESD(0)\r
+ , fMC(0)\r
+ , fESDfriend(0)\r
+ , fOutput(0)\r
+ , fPitList(0)\r
+ , fUseMCInfo(kFALSE)\r
+ , fdNdPtEventCuts(0)\r
+ , fdNdPtAcceptanceCuts(0)\r
+ , fdNdPtRecAcceptanceCuts(0)\r
+ , fEsdTrackCuts(0)\r
+ , fTrigger(AliTriggerAnalysis::kMB1) \r
+ , fAnalysisMode(AlidNdPtHelper::kTPC) \r
+ , fOutputSummary(0)\r
+ , fTreeSRedirector(0)\r
+ , fCentralityEstimator(0)\r
+{\r
+ // Constructor\r
+\r
+ // Define input and output slots here\r
+ DefineOutput(0, TTree::Class());\r
+ //DefineOutput(1, TList::Class());\r
+\r
+}\r
+\r
+//_____________________________________________________________________________\r
+AlidNdPtTrackDumpTask::~AlidNdPtTrackDumpTask()\r
+{\r
+ if(fOutput) delete fOutput; fOutput =0; \r
+ if(fOutputSummary) delete fOutputSummary; fOutputSummary =0; \r
+ if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector =0; \r
+\r
+ if(fdNdPtEventCuts) delete fdNdPtEventCuts; fdNdPtEventCuts=NULL; \r
+ if(fdNdPtAcceptanceCuts) delete fdNdPtAcceptanceCuts; fdNdPtAcceptanceCuts=NULL;\r
+ if(fdNdPtRecAcceptanceCuts) delete fdNdPtRecAcceptanceCuts; fdNdPtRecAcceptanceCuts=NULL; \r
+ if(fEsdTrackCuts) delete fEsdTrackCuts; fEsdTrackCuts=NULL;\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::Notify()\r
+{\r
+ static Int_t count = 0;\r
+ count++;\r
+ TTree *chain = (TChain*)GetInputData(0);\r
+ if(chain)\r
+ {\r
+ Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName());\r
+ }\r
+\r
+return kTRUE;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AlidNdPtTrackDumpTask::UserCreateOutputObjects()\r
+{\r
+ // Create histograms\r
+ // Called once\r
+ fOutput = new TList;\r
+ fOutput->SetOwner();\r
+\r
+ //\r
+ // create output tree\r
+ //\r
+ fTreeSRedirector = new TTreeSRedirector("dNdPtOutliersAnalysisPbPb.root");\r
+\r
+ PostData(0, fOutputSummary);\r
+ //PostData(1, fOutput);\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AlidNdPtTrackDumpTask::UserExec(Option_t *) \r
+{\r
+ //\r
+ // Called for each event\r
+ //\r
+\r
+ // ESD event\r
+ fESD = (AliESDEvent*) (InputEvent());\r
+ if (!fESD) {\r
+ Printf("ERROR: ESD event not available");\r
+ return;\r
+ }\r
+\r
+ // MC event\r
+ if(fUseMCInfo) {\r
+ fMC = MCEvent();\r
+ if (!fMC) {\r
+ Printf("ERROR: MC event not available");\r
+ return;\r
+ }\r
+ }\r
+\r
+ fESDfriend = static_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));\r
+ if(!fESDfriend) {\r
+ Printf("ERROR: ESD friends not available");\r
+ }\r
+\r
+ //\r
+ Process(fESD,fMC,fESDfriend);\r
+\r
+ // Post output data.\r
+ PostData(0, fOutputSummary);\r
+ //PostData(1, fOutput);\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AlidNdPtTrackDumpTask::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend)\r
+{\r
+ //\r
+ // Process real and/or simulated events\r
+ //\r
+ if(!esdEvent) {\r
+ AliDebug(AliLog::kError, "esdEvent not available");\r
+ return;\r
+ }\r
+\r
+ // get selection cuts\r
+ AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
+ AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
+ AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
+\r
+ if(!evtCuts || !accCuts || !esdTrackCuts) {\r
+ AliDebug(AliLog::kError, "cuts not available");\r
+ return;\r
+ }\r
+\r
+ // trigger selection\r
+ Bool_t isEventTriggered = kTRUE;\r
+ AliPhysicsSelection *physicsSelection = NULL;\r
+ AliTriggerAnalysis* triggerAnalysis = NULL;\r
+\r
+ // \r
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
+ if (!inputHandler)\r
+ {\r
+ Printf("ERROR: Could not receive input handler");\r
+ return;\r
+ }\r
+\r
+ if(evtCuts->IsTriggerRequired()) \r
+ {\r
+ // always MB\r
+ isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;\r
+\r
+ physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());\r
+ if(!physicsSelection) return;\r
+ //SetPhysicsTriggerSelection(physicsSelection);\r
+\r
+ if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {\r
+ // set trigger (V0AND)\r
+ triggerAnalysis = physicsSelection->GetTriggerAnalysis();\r
+ if(!triggerAnalysis) return;\r
+ isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
+ }\r
+ }\r
+\r
+ // centrality determination\r
+ Float_t centralityF = -1;\r
+ AliCentrality *esdCentrality = esdEvent->GetCentrality();\r
+ centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
+\r
+ // use MC information\r
+ AliHeader* header = 0;\r
+ AliGenEventHeader* genHeader = 0;\r
+ AliStack* stack = 0;\r
+ TArrayF vtxMC(3);\r
+\r
+ Int_t multMCTrueTracks = 0;\r
+ if(IsUseMCInfo())\r
+ {\r
+ //\r
+ if(!mcEvent) {\r
+ AliDebug(AliLog::kError, "mcEvent not available");\r
+ return;\r
+ }\r
+ // get MC event header\r
+ header = mcEvent->Header();\r
+ if (!header) {\r
+ AliDebug(AliLog::kError, "Header not available");\r
+ return;\r
+ }\r
+ // MC particle stack\r
+ stack = mcEvent->Stack();\r
+ if (!stack) {\r
+ AliDebug(AliLog::kError, "Stack not available");\r
+ return;\r
+ }\r
+\r
+ // get MC vertex\r
+ genHeader = header->GenEventHeader();\r
+ if (!genHeader) {\r
+ AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
+ return;\r
+ }\r
+ genHeader->PrimaryVertex(vtxMC);\r
+\r
+ // multipliticy of all MC primary tracks\r
+ // in Zv, pt and eta ranges)\r
+ multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
+\r
+ } // end bUseMC\r
+\r
+ // get reconstructed vertex \r
+ //const AliESDVertex* vtxESD = 0; \r
+ const AliESDVertex* vtxESD = 0; \r
+ if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
+ vtxESD = esdEvent->GetPrimaryVertexTPC();\r
+ }\r
+ else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
+ vtxESD = esdEvent->GetPrimaryVertexTracks();\r
+ }\r
+ else {\r
+ return;\r
+ }\r
+\r
+ if(!vtxESD) return;\r
+\r
+ Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
+ //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
+ //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
+\r
+ // check event cuts\r
+ if(isEventOK && isEventTriggered)\r
+ {\r
+\r
+ for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)\r
+ {\r
+ AliESDtrack *track = esdEvent->GetTrack(iTrack);\r
+ if(!track) continue;\r
+ if(track->Charge()==0) continue;\r
+ if(!esdTrackCuts->AcceptTrack(track)) continue;\r
+ if(!accCuts->AcceptTrack(track)) continue;\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
+ Bool_t isOK = kFALSE;\r
+\r
+ //\r
+ // Constrain TPC-only tracks (TPCinner) to vertex\r
+ //\r
+ AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam());\r
+ if (!tpcInner) continue;\r
+ // transform to the track reference frame \r
+ isOK = tpcInner->Rotate(track->GetAlpha());\r
+ isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
+ if(!isOK) continue;\r
+\r
+ // clone TPCinner has to be deleted\r
+ AliExternalTrackParam * tpcInnerC = (AliExternalTrackParam *)(track->GetTPCInnerParam()->Clone());\r
+ if (!tpcInnerC) continue;\r
+ \r
+ // constrain TPCinner \r
+ //isOK = ConstrainTPCInner(tpcInnerC,vtxESD,esdEvent->GetMagneticField());\r
+ isOK = ConstrainTPCInner(tpcInnerC,vtxESD,b);\r
+\r
+ // transform to the track reference frame \r
+ isOK = tpcInnerC->Rotate(track->GetAlpha());\r
+ isOK = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
+\r
+ if(!isOK) {\r
+ if(tpcInnerC) delete tpcInnerC;\r
+ continue;\r
+ }\r
+\r
+\r
+ //\r
+ // Constrain TPC refitted tracks at inner TPC wall (InnerParams) \r
+ // to vertex\r
+ //\r
+ // clone track InnerParams has to be deleted\r
+ AliExternalTrackParam * trackInnerC = (AliExternalTrackParam *)(track->GetInnerParam()->Clone());\r
+ if (!trackInnerC) continue;\r
+ \r
+ // constrain track InnerParams \r
+ isOK = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b);\r
+\r
+ // transform to the track reference frame \r
+ isOK = trackInnerC->Rotate(track->GetAlpha());\r
+ isOK = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField());\r
+\r
+ if(!isOK) {\r
+ if(trackInnerC) delete trackInnerC;\r
+ continue;\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
+\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
+ // chi2 distance TPC constrained and TPC+ITS\r
+ TMatrixD covarMInv = covarM.Invert();\r
+ TMatrixD mat2 = covarMInv*deltaT;\r
+ TMatrixD 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
+ TMatrixD chi2trackC = deltatrackC*mat2trackC; \r
+ //chi2trackC.Print();\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
+ AliExternalTrackParam *trackInnerC2 = (AliExternalTrackParam *)(track->GetInnerParam()->Clone());\r
+ if (!trackInnerC2) continue;\r
+ if(!AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
+ {\r
+ if(trackInnerC2) delete trackInnerC2;\r
+ continue;\r
+ }\r
+\r
+ AliExternalTrackParam *outerITSc = new AliExternalTrackParam();\r
+ if(!outerITSc) continue;\r
+\r
+ TMatrixD chi2OuterITS(1,1);\r
+ chi2OuterITS(0,0) = 0;\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
+ if( (outerITSc = (AliExternalTrackParam *)(friendTrack->GetITSOut()->Clone())) ) \r
+ {\r
+ if(AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE))\r
+ {\r
+ // transform ITSout to the track InnerParams reference frame \r
+ isOK = kFALSE;\r
+ isOK = outerITSc->Rotate(trackInnerC2->GetAlpha());\r
+ isOK = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField());\r
+\r
+ if(!isOK) {\r
+ if(outerITSc) delete outerITSc;\r
+ if(trackInnerC2) delete trackInnerC2;\r
+ continue;\r
+ }\r
+ \r
+ //\r
+ // calculate chi2 between outerITS and innerParams\r
+ // cov without z-coordinate at the moment\r
+ //\r
+ TMatrixD deltaTouterITS(4,1);\r
+ TMatrixD deltaouterITS(1,4);\r
+ TMatrixD covarMouterITS(4,4);\r
+\r
+ Int_t kipar = 0;\r
+ Int_t kjpar = 0;\r
+ for (Int_t ipar=0; ipar<5; ipar++) {\r
+ if(ipar!=1) {\r
+ deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
+ deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar];\r
+ }\r
+\r
+ kjpar=0;\r
+ for (Int_t jpar=0; jpar<5; jpar++) {\r
+ Int_t index=outerITSc->GetIndex(ipar,jpar);\r
+ if(ipar !=1 || jpar!=1) {\r
+ covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index];\r
+ }\r
+ if(jpar!=1) kjpar++;\r
+ }\r
+ if(ipar!=1) kipar++;\r
+ }\r
+\r
+ // chi2 distance ITSout and InnerParams\r
+ TMatrixD covarMInvouterITS = covarMouterITS.Invert();\r
+ TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS;\r
+ 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
+ AliExternalTrackParam *trackInnerC3 = (AliExternalTrackParam *)(track->GetInnerParam()->Clone());\r
+ if(!trackInnerC3) continue;\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
+ //printf("ref %p \n",ref);\r
+ //if(ref) printf("ref->DetectorId() %d \n",ref->DetectorId());\r
+ //printf("AliTrackReference::kTPC %d \n",AliTrackReference::kTPC);\r
+\r
+\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
+ isOK = kFALSE;\r
+ isOK = trackInnerC3->Rotate(kRefPhi);\r
+ isOK = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE);\r
+\r
+ if(!isOK){\r
+ if(trackInnerC3) delete trackInnerC3;\r
+ if(refTPCIn) delete refTPCIn;\r
+ }\r
+ }\r
+ }\r
+\r
+ //\r
+ // ITS track\r
+ //\r
+ Int_t labelITS = TMath::Abs(track->GetITSLabel()); \r
+ particleITS = stack->Particle(labelITS);\r
+ if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.)\r
+ {\r
+ particleMotherITS = GetMother(particleITS,stack);\r
+ mechITS = particleITS->GetUniqueID();\r
+ isPrimITS = stack->IsPhysicalPrimary(labelITS);\r
+ isFromStrangessITS = IsFromStrangeness(labelITS,stack);\r
+ isFromConversionITS = IsFromConversion(labelITS,stack);\r
+ isFromMaterialITS = IsFromMaterial(labelITS,stack);\r
+ }\r
+ }\r
+\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
+\r
+ //\r
+ if(!fTreeSRedirector) return;\r
+ (*fTreeSRedirector)<<"dNdPtTree"<<\r
+ "runNumber="<<runNumber<<\r
+ "Bz="<<bz<<\r
+ "vertX="<<vert[0]<<\r
+ "vertY="<<vert[1]<<\r
+ "vertZ="<<vert[2]<<\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
+ if(tpcInnerC) delete tpcInnerC;\r
+ if(trackInnerC) delete trackInnerC;\r
+ if(trackInnerC2) delete trackInnerC2;\r
+ if(outerITSc) delete outerITSc;\r
+\r
+ if(trackInnerC3) delete trackInnerC3;\r
+ }\r
+ }\r
+\r
+ PostData(0, fOutputSummary);\r
+ //PostData(1, fOutput);\r
+}\r
+\r
+\r
+//_____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3])\r
+{\r
+ // Constrain TPC inner params constrained\r
+ //\r
+ if(!tpcInnerC) return kFALSE; \r
+ if(!vtx) return kFALSE;\r
+\r
+ Double_t dz[2],cov[3];\r
+ //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
+ //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
+ //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
+ if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
+\r
+\r
+ Double_t covar[6]; vtx->GetCovMatrix(covar);\r
+ Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]};\r
+ Double_t c[3]={covar[2],0.,covar[5]};\r
+ Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c);\r
+ if (chi2C>kVeryBig) return kFALSE; \r
+\r
+ if(!tpcInnerC->Update(p,c)) return kFALSE;\r
+\r
+ return kTRUE;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3])\r
+{\r
+ // Constrain TPC inner params constrained\r
+ //\r
+ if(!trackInnerC) return kFALSE; \r
+ if(!vtx) return kFALSE;\r
+\r
+ const Double_t kRadius = 2.8; \r
+ const Double_t kMaxStep = 1.0; \r
+\r
+ Double_t dz[2],cov[3];\r
+\r
+ //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex();\r
+ //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; \r
+ //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; \r
+\r
+ if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE;\r
+ if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; \r
+\r
+ //\r
+ Double_t covar[6]; vtx->GetCovMatrix(covar);\r
+ Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]};\r
+ Double_t c[3]={covar[2],0.,covar[5]};\r
+ Double_t chi2C=trackInnerC->GetPredictedChi2(p,c);\r
+ if (chi2C>kVeryBig) return kFALSE; \r
+\r
+ if(!trackInnerC->Update(p,c)) return kFALSE;\r
+\r
+ return kTRUE;\r
+}\r
+\r
+\r
+//_____________________________________________________________________________\r
+TParticle *AlidNdPtTrackDumpTask::GetMother(TParticle *const particle, AliStack *const stack) \r
+{\r
+ if(!particle) return NULL;\r
+ if(!stack) return NULL;\r
+\r
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
+ TParticle* mother = NULL; \r
+ mother = stack->Particle(motherLabel); \r
+\r
+return mother;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::IsFromConversion(const Int_t label, AliStack *const stack) \r
+{\r
+ Bool_t isFromConversion = kFALSE;\r
+\r
+ if(stack) {\r
+ TParticle* particle = stack->Particle(label);\r
+\r
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
+ {\r
+ Int_t mech = particle->GetUniqueID(); // production mechanism \r
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
+\r
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
+ TParticle* mother = stack->Particle(motherLabel); \r
+ if(mother) {\r
+ Int_t motherPdg = mother->GetPdgCode();\r
+\r
+ if(!isPrim && mech==5 && motherPdg==kGamma) { \r
+ isFromConversion=kTRUE; \r
+ }\r
+ }\r
+ } \r
+ } \r
+\r
+ return isFromConversion;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::IsFromMaterial(const Int_t label, AliStack *const stack) \r
+{\r
+ Bool_t isFromMaterial = kFALSE;\r
+\r
+ if(stack) {\r
+ TParticle* particle = stack->Particle(label);\r
+\r
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
+ {\r
+ Int_t mech = particle->GetUniqueID(); // production mechanism \r
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
+\r
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
+ TParticle* mother = stack->Particle(motherLabel); \r
+ if(mother) {\r
+ if(!isPrim && mech==13) { \r
+ isFromMaterial=kTRUE; \r
+ }\r
+ }\r
+ } \r
+ } \r
+\r
+ return isFromMaterial;\r
+}\r
+\r
+//_____________________________________________________________________________\r
+Bool_t AlidNdPtTrackDumpTask::IsFromStrangeness(const Int_t label, AliStack *const stack) \r
+{\r
+ Bool_t isFromStrangeness = kFALSE;\r
+\r
+ if(stack) {\r
+ TParticle* particle = stack->Particle(label);\r
+\r
+ if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) \r
+ {\r
+ Int_t mech = particle->GetUniqueID(); // production mechanism \r
+ Bool_t isPrim = stack->IsPhysicalPrimary(label);\r
+\r
+ Int_t motherLabel = TMath::Abs(particle->GetMother(0)); \r
+ TParticle* mother = stack->Particle(motherLabel); \r
+ if(mother) {\r
+ Int_t motherPdg = mother->GetPdgCode();\r
+\r
+ // K+-, lambda, antilambda, K0s decays\r
+ if(!isPrim && mech==4 && \r
+ (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short))\r
+ {\r
+ isFromStrangeness = kTRUE;\r
+ } \r
+ }\r
+ } \r
+ } \r
+\r
+ return isFromStrangeness;\r
+}\r
+\r
+\r
+//_____________________________________________________________________________\r
+void AlidNdPtTrackDumpTask::FinishTaskOutput() \r
+{\r
+ //\r
+ // Called one at the end \r
+ // locally on working node\r
+ //\r
+\r
+ // Post output data.\r
+ PostData(1, fOutput);\r
+ //PostData(0, fOutputSummary);\r
+}\r
+\r
+//_____________________________________________________________________________\r
+void AlidNdPtTrackDumpTask::Terminate(Option_t *) \r
+{\r
+ // Called one at the end \r
+ \r
+ // check output data\r
+ fOutputSummary = dynamic_cast<TTree*> (GetOutputData(0));\r
+ if(fOutputSummary) delete fOutputSummary; fOutputSummary=0;\r
+ if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0;\r
+\r
+ TChain* chain = new TChain("dNdPtTree");\r
+ if(!chain) return;\r
+ chain->Add("dNdPtOutliersAnalysisPbPb.root");\r
+ TTree *tree = chain->CopyTree("1");\r
+ if (chain) { delete chain; chain=0; }\r
+ if(!tree) return;\r
+ tree->Print();\r
+ fOutputSummary = tree;\r
+\r
+ if (!fOutputSummary) {\r
+ Printf("ERROR: AlidNdPtTrackDumpTask::Terminate(): Output data not avaiable GetOutputData(0)==0x0 ..." );\r
+ return;\r
+ }\r
+ \r
+\r
+\r
+ PostData(0, fOutputSummary);\r
+ //PostData(1, fOutput);\r
+}\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff