/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /* Class to process/filter reconstruction information from ESD, ESD friends, MC and provide them for later reprocessing Filtering schema - low pt part is downscaled - to have flat pt specatra of selected topologies (tracks and V0s) Downscaling schema is controlled by downscaling factors Usage: 1.) Filtering on Lego train 2.) expert QA for tracking (resolution efficnecy) 3.) pt reoslution studies using V0s 4.) dEdx calibration using V0s 5.) pt resolution and dEdx studies using cosmic + 6.) Info used for later raw data OFFLINE triggering (highPt, V0, laser, cosmic, high dEdx) Exported trees (with full objects and dereived variables): 1.) "highPt" - filtered trees with esd tracks, derived variables(propagated tracks), optional MC info +optional space points 2.) "V0s" - - filtered trees with selected V0s (rough KF chi2 cut), KF particle and corresponding esd tracks + optionla space points 3.) "Laser" - dump laser tracks with space points if exests 4.) "CosmicTree" - cosmic track candidate (random or triggered) + esdTracks(up/down)+ optional points 5.) "dEdx" - tree with high dEdx tpc tracks */ #include "iostream" #include "TSystem.h" #include #include #include "TChain.h" #include "TTreeStream.h" #include "TTree.h" #include "TH1F.h" #include "TH3.h" #include "TCanvas.h" #include "TList.h" #include "TObjArray.h" #include "TFile.h" #include "TMatrixD.h" #include "TRandom3.h" #include "AliHeader.h" #include "AliGenEventHeader.h" #include "AliInputEventHandler.h" #include "AliStack.h" #include "AliTrackReference.h" #include "AliTrackPointArray.h" #include "AliSysInfo.h" #include "AliPhysicsSelection.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliESDEvent.h" #include "AliESDfriend.h" #include "AliMCEvent.h" #include "AliESDInputHandler.h" #include "AliESDVertex.h" #include "AliTracker.h" #include "AliVTrack.h" #include "AliGeomManager.h" #include "AliCentrality.h" #include "AliESDVZERO.h" #include "AliMultiplicity.h" #include "AliESDtrackCuts.h" #include "AliMCEventHandler.h" #include "AliFilteredTreeEventCuts.h" #include "AliFilteredTreeAcceptanceCuts.h" #include "AliAnalysisTaskFilteredTree.h" #include "AliKFParticle.h" #include "AliESDv0.h" #include "TVectorD.h" using namespace std; ClassImp(AliAnalysisTaskFilteredTree) //_____________________________________________________________________________ AliAnalysisTaskFilteredTree::AliAnalysisTaskFilteredTree(const char *name) : AliAnalysisTaskSE(name) , fESD(0) , fMC(0) , fESDfriend(0) , fOutput(0) , fPitList(0) , fUseMCInfo(kFALSE) , fUseESDfriends(kFALSE) , fReducePileUp(kTRUE) , fFillTree(kTRUE) , fFilteredTreeEventCuts(0) , fFilteredTreeAcceptanceCuts(0) , fFilteredTreeRecAcceptanceCuts(0) , fEsdTrackCuts(0) , fTrigger(AliTriggerAnalysis::kMB1) , fAnalysisMode(kTPCAnalysisMode) , fTreeSRedirector(0) , fCentralityEstimator(0) , fLowPtTrackDownscaligF(0) , fLowPtV0DownscaligF(0) , fFriendDownscaling(-3.) , fProcessAll(kFALSE) , fProcessCosmics(kFALSE) , fProcessITSTPCmatchOut(kFALSE) // swittch to process ITS/TPC standalone tracks , fHighPtTree(0) , fV0Tree(0) , fdEdxTree(0) , fLaserTree(0) , fMCEffTree(0) , fCosmicPairsTree(0) , fPtResPhiPtTPC(0) , fPtResPhiPtTPCc(0) , fPtResPhiPtTPCITS(0) , fPtResEtaPtTPC(0) , fPtResEtaPtTPCc(0) , fPtResEtaPtTPCITS(0) , fPtResCentPtTPC(0) , fPtResCentPtTPCc(0) , fPtResCentPtTPCITS(0) , fCurrentFileName("") , fDummyTrack(0) { // Constructor // Define input and output slots here DefineOutput(1, TTree::Class()); DefineOutput(2, TTree::Class()); DefineOutput(3, TTree::Class()); DefineOutput(4, TTree::Class()); DefineOutput(5, TTree::Class()); DefineOutput(6, TTree::Class()); DefineOutput(7, TList::Class()); } //_____________________________________________________________________________ AliAnalysisTaskFilteredTree::~AliAnalysisTaskFilteredTree() { // // Destructor // delete fFilteredTreeEventCuts; delete fFilteredTreeAcceptanceCuts; delete fFilteredTreeRecAcceptanceCuts; delete fEsdTrackCuts; } //____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::Notify() { // // static Int_t count = 0; count++; TTree *chain = (TChain*)GetInputData(0); if(chain){ Printf("Processing %d. file: %s", count, chain->GetCurrentFile()->GetName()); } fCurrentFileName=chain->GetCurrentFile()->GetName(); return kTRUE; } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::UserCreateOutputObjects() { // Create histograms // Called once // //get the output file to make sure the trees will be associated to it OpenFile(1); fTreeSRedirector = new TTreeSRedirector(); // // Create trees fV0Tree = ((*fTreeSRedirector)<<"V0s").GetTree(); fHighPtTree = ((*fTreeSRedirector)<<"highPt").GetTree(); fdEdxTree = ((*fTreeSRedirector)<<"dEdx").GetTree(); fLaserTree = ((*fTreeSRedirector)<<"Laser").GetTree(); fMCEffTree = ((*fTreeSRedirector)<<"MCEffTree").GetTree(); fCosmicPairsTree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree(); if (!fDummyTrack) { fDummyTrack=new AliESDtrack(); } // histogram booking Double_t minPt = 0.1; Double_t maxPt = 100.; Int_t nbinsPt = 30; Double_t logminPt = TMath::Log10(minPt); Double_t logmaxPt = TMath::Log10(maxPt); Double_t binwidth = (logmaxPt-logminPt)/nbinsPt; Double_t *binsPt = new Double_t[nbinsPt+1]; binsPt[0] = minPt; for (Int_t i=1;i<=nbinsPt;i++) { binsPt[i] = minPt + TMath::Power(10,logminPt+i*binwidth); } // 1pT resol cov matrix bins Double_t min1PtRes = 0.; Double_t max1PtRes = 0.3; Int_t nbins1PtRes = 300; Double_t bins1PtRes[301]; for (Int_t i=0;i<=nbins1PtRes;i++) { bins1PtRes[i] = min1PtRes + i*(max1PtRes-min1PtRes)/nbins1PtRes; } // phi bins Double_t minPhi = 0.; Double_t maxPhi = 6.5; Int_t nbinsPhi = 100; Double_t binsPhi[101]; for (Int_t i=0;i<=nbinsPhi;i++) { binsPhi[i] = minPhi + i*(maxPhi-minPhi)/nbinsPhi; } // eta bins Double_t minEta = -1.; Double_t maxEta = 1.; Int_t nbinsEta = 20; Double_t binsEta[21]; for (Int_t i=0;i<=nbinsEta;i++) { binsEta[i] = minEta + i*(maxEta-minEta)/nbinsEta; } // mult bins Double_t minCent = 0.; Double_t maxCent = 100; Int_t nbinsCent = 20; Double_t binsCent[101]; for (Int_t i=0;i<=nbinsCent;i++) { binsCent[i] = minCent + i*(maxCent-minCent)/nbinsCent; } fPtResPhiPtTPC = new TH3D("fPtResPhiPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes); fPtResPhiPtTPCc = new TH3D("fPtResPhiPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes); fPtResPhiPtTPCITS = new TH3D("fPtResPhiPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsPhi,binsPhi,nbins1PtRes,bins1PtRes); fPtResEtaPtTPC = new TH3D("fPtResEtaPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes); fPtResEtaPtTPCc = new TH3D("fPtResEtaPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes); fPtResEtaPtTPCITS = new TH3D("fPtResEtaPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsEta,binsEta,nbins1PtRes,bins1PtRes); fPtResCentPtTPC = new TH3D("fPtResCentPtTPC","pt rel. resolution from cov. matrix TPC tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes); fPtResCentPtTPCc = new TH3D("fPtResCentPtTPCc","pt rel. resolution from cov. matrix TPC constrained tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes); fPtResCentPtTPCITS = new TH3D("fPtResCentPtTPCITS","pt rel. resolution from cov. matrix TPC+ITS tracks",nbinsPt,binsPt,nbinsCent,binsCent,nbins1PtRes,bins1PtRes); fOutput = new TList; if(!fOutput) return; fOutput->SetOwner(); fOutput->Add(fPtResPhiPtTPC); fOutput->Add(fPtResPhiPtTPCc); fOutput->Add(fPtResPhiPtTPCITS); fOutput->Add(fPtResEtaPtTPC); fOutput->Add(fPtResEtaPtTPCc); fOutput->Add(fPtResEtaPtTPCITS); fOutput->Add(fPtResCentPtTPC); fOutput->Add(fPtResCentPtTPCc); fOutput->Add(fPtResCentPtTPCITS); // post data to outputs PostData(1,fV0Tree); PostData(2,fHighPtTree); PostData(3,fdEdxTree); PostData(4,fLaserTree); PostData(5,fMCEffTree); PostData(6,fCosmicPairsTree); PostData(7,fOutput); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::UserExec(Option_t *) { // // Called for each event // // ESD event fESD = (AliESDEvent*) (InputEvent()); if (!fESD) { Printf("ERROR: ESD event not available"); return; } //if MC info available - use it. fMC = MCEvent(); if(fUseESDfriends) { //fESDfriend = dynamic_cast(fESD->FindListObject("AliESDfriend")); fESDfriend = ESDfriend(); if(!fESDfriend) { Printf("ERROR: ESD friends not available"); } } AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); if (!inputHandler){ return; } //if set, use the environment variables to set the downscaling factors //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF //AliAnalysisTaskFilteredTree_fFriendDownscaling TString env; env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF"); if (!env.IsNull()){ fLowPtTrackDownscaligF=env.Atof(); AliInfo(Form("fLowPtTrackDownscaligF=%f",fLowPtTrackDownscaligF)); } env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF"); if (!env.IsNull()){ fLowPtV0DownscaligF=env.Atof(); AliInfo(Form("fLowPtV0DownscaligF=%f",fLowPtTrackDownscaligF)); } env = gSystem->Getenv("AliAnalysisTaskFilteredTree_fFriendDownscaling"); if (!env.IsNull()){ fFriendDownscaling=env.Atof(); AliInfo(Form(" fFriendDownscaling=%f",fFriendDownscaling)); } // // // if(fProcessAll) { ProcessAll(fESD,fMC,fESDfriend); // all track stages and MC } else { Process(fESD,fMC,fESDfriend); // only global and TPC tracks } // ProcessV0(fESD,fMC,fESDfriend); ProcessLaser(fESD,fMC,fESDfriend); ProcessdEdx(fESD,fMC,fESDfriend); if (fProcessCosmics) { ProcessCosmics(fESD,fESDfriend); } if(fMC) { ProcessMCEff(fESD,fMC,fESDfriend);} if (fProcessITSTPCmatchOut) ProcessITSTPCmatchOut(fESD, fESDfriend); printf("processed event %d\n", Int_t(Entry())); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend) { // // Find cosmic pairs (triggered or random) // // AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); const Double_t kMinPt=0.8; const Double_t kMinPtMax=0.8; const Double_t kMinNcl=50; const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; Int_t ntracks=event->GetNumberOfTracks(); UInt_t specie = event->GetEventSpecie(); // skip laser events if (specie==AliRecoParam::kCalib) return; Int_t ntracksFriend = esdFriend ? esdFriend->GetNumberOfTracks() : 0; for (Int_t itrack0=0;itrack0GetTrack(itrack0); if (!track0) continue; if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue; if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue; if (track0->Pt() < kMinPt) continue; if (track0->GetTPCncls() < kMinNcl) continue; if (TMath::Abs(track0->GetY())GetKinkIndex(0)>0) continue; const Double_t * par0=track0->GetParameter(); //track param at rhe DCA //rm primaries // //track0->GetImpactParametersTPC(dcaTPC,covTPC); //if (TMath::Abs(dcaTPC[0])GetInnerParam(); AliESDfriendTrack* friendTrack0=NULL; if (esdFriend &&!esdFriend->TestSkipBit()){ if (itrack0GetTrack(itrack0); } //this guy can be NULL } for (Int_t itrack1=itrack0+1;itrack1GetTrack(itrack1); if (!track1) continue; if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue; if (track1->GetKinkIndex(0)>0) continue; if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue; if (track1->Pt() < kMinPt) continue; if (track1->GetTPCncls()1 && TMath::Max(track1->Pt(), track0->Pt())GetY())GetImpactParametersTPC(dcaTPC,covTPC); // if (TMath::Abs(dcaTPC[0])GetParameter(); //track param at rhe DCA // Bool_t isPair=kTRUE; for (Int_t ipar=0; ipar<5; ipar++){ if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE; } if (!isPair) continue; if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; //delta with correct sign /* TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" */ if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign if (!isPair) continue; TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); Int_t eventNumber = event->GetEventNumberInFile(); // // Int_t ntracksSPD = vertexSPD->GetNContributors(); Int_t ntracksTPC = vertexTPC->GetNContributors(); Int_t runNumber = event->GetRunNumber(); Int_t timeStamp = event->GetTimeStamp(); ULong64_t triggerMask = event->GetTriggerMask(); Float_t magField = event->GetMagneticField(); TObjString triggerClass = event->GetFiredTriggerClasses().Data(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)event->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)event->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)event->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); AliESDfriendTrack* friendTrack1=NULL; if (esdFriend &&!esdFriend->TestSkipBit()){ if (itrack1GetTrack(itrack1); } //this guy can be NULL } // AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing if (fFriendDownscaling>=1){ // downscaling number of friend tracks if (gRandom->Rndm()>1./fFriendDownscaling){ friendTrackStore0 = 0; friendTrackStore1 = 0; } } if (fFriendDownscaling<=0){ if (((*fTreeSRedirector)<<"CosmicPairs").GetTree()){ TTree * tree = ((*fTreeSRedirector)<<"CosmicPairs").GetTree(); if (tree){ Double_t sizeAll=tree->GetZipBytes(); TBranch * br= tree->GetBranch("friendTrack0.fPoints"); Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; br= tree->GetBranch("friendTrack0.fCalibContainer"); if (br) sizeFriend+=br->GetZipBytes(); if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) { friendTrackStore0=0; friendTrackStore1=0; } } } } if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"CosmicPairs"<< "gid="<GetInputEventHandler(); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; //SetPhysicsTriggerSelection(physicsSelection); if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // use MC information AliHeader* header = 0; AliGenEventHeader* genHeader = 0; AliStack* stack = 0; TArrayF vtxMC(3); Int_t multMCTrueTracks = 0; if(mcEvent) { // get MC event header header = mcEvent->Header(); if (!header) { AliDebug(AliLog::kError, "Header not available"); return; } // MC particle stack stack = mcEvent->Stack(); if (!stack) { AliDebug(AliLog::kError, "Stack not available"); return; } // get MC vertex genHeader = header->GenEventHeader(); if (!genHeader) { AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); return; } genHeader->PrimaryVertex(vtxMC); // multipliticy of all MC primary tracks // in Zv, pt and eta ranges) multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); } // end bUseMC // get reconstructed vertex //const AliESDVertex* vtxESD = 0; AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } AliESDVertex* vtxTPC = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); AliESDVertex* vtxSPD = (AliESDVertex*)esdEvent->GetPrimaryVertexSPD(); if(!vtxESD) return; if(!vtxTPC) return; if(!vtxSPD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d, status %d, vz %f \n",isEventOK, isEventTriggered, vtxESD->GetStatus(), vtxESD->GetZ()); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); Int_t ntracks = esdEvent->GetNumberOfTracks(); // check event cuts if(isEventOK && isEventTriggered) { // // get IR information // AliESDHeader *esdHeader = 0; esdHeader = esdEvent->GetHeader(); if(!esdHeader) return; //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval // Use when Peter commit the changes in the header Int_t ir1 = 0; Int_t ir2 = 0; // //Double_t vert[3] = {0}; //vert[0] = vtxESD->GetX(); //vert[1] = vtxESD->GetY(); //vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Int_t multSPD = vtxSPD->GetNContributors(); Int_t multTPC = vtxTPC->GetNContributors(); Float_t bz = esdEvent->GetMagneticField(); Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); // high pT tracks for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; if(track->Charge()==0) continue; if(!esdTrackCuts->AcceptTrack(track)) continue; if(!accCuts->AcceptTrack(track)) continue; // downscale low-pT tracks Double_t scalempt= TMath::Min(track->Pt(),10.); Double_t downscaleF = gRandom->Rndm(); downscaleF *= fLowPtTrackDownscaligF; if(TMath::Exp(2*scalempt)GetTPCInnerParam()); if (!tpcInner) continue; // transform to the track reference frame Bool_t isOK = kFALSE; isOK = tpcInner->Rotate(track->GetAlpha()); isOK = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); if(!isOK) continue; // Dump to the tree // vertex // TPC-ITS tracks // TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"highPt"<< "gid="< dump tracks and clusters to the special tree // const Double_t kMinPt = 5; if(!fFillTree) return; if(!fTreeSRedirector) return; const AliESDHeader* esdHeader = esdEvent->GetHeader(); if(esdHeader && esdHeader->GetEventSpecie()==AliRecoParam::kCalib) { Int_t countLaserTracks = 0; Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); Float_t bz = esdEvent->GetMagneticField(); TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){ AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; if(track->GetTPCInnerParam()) countLaserTracks++; AliESDfriendTrack* friendTrack=NULL; // suppress beam background and CE random reacks if (track->GetInnerParam()->Pt()Rndm()>1/(1+TMath::Abs(fFriendDownscaling)); if (skipTrack) continue; if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL (*fTreeSRedirector)<<"Laser"<< "gid="<GetInputEventHandler(); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) {AliInfo("no physics selection"); return;} //SetPhysicsTriggerSelection(physicsSelection); if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) {AliInfo("no trigger analysis");return;} isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // use MC information AliHeader* header = 0; AliGenEventHeader* genHeader = 0; AliStack* stack = 0; TArrayF vtxMC(3); Int_t mcStackSize=0; Int_t multMCTrueTracks = 0; if(mcEvent) { // get MC event header header = mcEvent->Header(); if (!header) { AliDebug(AliLog::kError, "Header not available"); return; } // MC particle stack stack = mcEvent->Stack(); if (!stack) { AliDebug(AliLog::kError, "Stack not available"); return; } mcStackSize=stack->GetNtrack(); // get MC vertex genHeader = header->GenEventHeader(); if (!genHeader) { AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); return; } genHeader->PrimaryVertex(vtxMC); // multipliticy of all MC primary tracks // in Zv, pt and eta ranges) multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); } // end bUseMC // get reconstructed vertex //const AliESDVertex* vtxESD = 0; AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { AliInfo("no ESD vertex"); return; } if(!vtxESD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); // // Vertex info comparison and track multiplicity // AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD(); AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC(); Int_t contSPD = vertexSPD->GetNContributors(); Int_t contTPC = vertexTPC->GetNContributors(); TVectorD vertexPosTPC(3), vertexPosSPD(3); vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray()); vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray()); Int_t ntracksTPC=0; Int_t ntracksITS=0; for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){ AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++; if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++; } //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); // check event cuts if(isEventOK && isEventTriggered) { // // get IR information // AliESDHeader *esdHeader = 0; esdHeader = esdEvent->GetHeader(); if(!esdHeader) {AliInfo("no esdHeader");return;} //Int_t ir1 = esdHeader->GetTriggerIREntries(); //all ir-s //Int_t ir2 = esdHeader->GetTriggerIREntries(-1,1); // int2 set, 180 ms time interval // Int_t ir1 = 0; Int_t ir2 = 0; // Double_t vert[3] = {0}; vert[0] = vtxESD->GetX(); vert[1] = vtxESD->GetY(); vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Float_t bz = esdEvent->GetMagneticField(); Int_t runNumber = esdEvent->GetRunNumber(); Int_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); Int_t numberOfTracks=esdEvent->GetNumberOfTracks(); // high pT tracks for (Int_t iTrack = 0; iTrack < numberOfTracks; iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); AliESDfriendTrack* friendTrack=NULL; Int_t numberOfFriendTracks=0; if (esdFriend) numberOfFriendTracks=esdFriend->GetNumberOfTracks(); if (esdFriend && iTrackTestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL if(!track) continue; if(track->Charge()==0) continue; if(!esdTrackCuts->AcceptTrack(track)) continue; if(!accCuts->AcceptTrack(track)) continue; // downscale low-pT tracks Double_t scalempt= TMath::Min(track->Pt(),10.); Double_t downscaleF = gRandom->Rndm(); downscaleF *= fLowPtTrackDownscaligF; if(TMath::Exp(2*scalempt)GetXYZ(x); Double_t b[3]; AliTracker::GetBxByBz(x,b); // // Transform TPC inner params to track reference frame // Bool_t isOKtpcInner = kFALSE; AliExternalTrackParam * tpcInner = (AliExternalTrackParam *)(track->GetTPCInnerParam()); if (tpcInner) { // transform to the track reference frame isOKtpcInner = tpcInner->Rotate(track->GetAlpha()); isOKtpcInner = tpcInner->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); } // // Constrain TPC inner to vertex // clone TPCinner has to be deleted // Bool_t isOKtpcInnerC = kFALSE; AliExternalTrackParam * tpcInnerC = new AliExternalTrackParam(*(track->GetTPCInnerParam())); if (tpcInnerC) { isOKtpcInnerC = ConstrainTPCInner(tpcInnerC,vtxESD,b); isOKtpcInnerC = tpcInnerC->Rotate(track->GetAlpha()); isOKtpcInnerC = tpcInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); } // // Constrain TPC refitted tracks at inner TPC wall (InnerParams) to vertex // Clone track InnerParams has to be deleted // Bool_t isOKtrackInnerC = kFALSE; AliExternalTrackParam * trackInnerC = new AliExternalTrackParam(*(track->GetInnerParam())); if (trackInnerC) { isOKtrackInnerC = ConstrainTrackInner(trackInnerC,vtxESD,track->GetMass(),b); isOKtrackInnerC = trackInnerC->Rotate(track->GetAlpha()); isOKtrackInnerC = trackInnerC->PropagateTo(track->GetX(),esdEvent->GetMagneticField()); } // // calculate chi2 between vi and vj vectors // with covi and covj covariance matrices // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj) // TMatrixD deltaT(5,1), deltaTtrackC(5,1); TMatrixD delta(1,5), deltatrackC(1,5); TMatrixD covarM(5,5), covarMtrackC(5,5); TMatrixD chi2(1,1); TMatrixD chi2trackC(1,1); if(isOKtpcInnerC && isOKtrackInnerC) { for (Int_t ipar=0; ipar<5; ipar++) { deltaT(ipar,0)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; delta(0,ipar)=tpcInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; deltaTtrackC(ipar,0)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; deltatrackC(0,ipar)=trackInnerC->GetParameter()[ipar]-track->GetParameter()[ipar]; for (Int_t jpar=0; jpar<5; jpar++) { Int_t index=track->GetIndex(ipar,jpar); covarM(ipar,jpar)=track->GetCovariance()[index]+tpcInnerC->GetCovariance()[index]; covarMtrackC(ipar,jpar)=track->GetCovariance()[index]+trackInnerC->GetCovariance()[index]; } } // chi2 distance TPC constrained and TPC+ITS TMatrixD covarMInv = covarM.Invert(); TMatrixD mat2 = covarMInv*deltaT; chi2 = delta*mat2; //chi2.Print(); // chi2 distance TPC refitted constrained and TPC+ITS TMatrixD covarMInvtrackC = covarMtrackC.Invert(); TMatrixD mat2trackC = covarMInvtrackC*deltaTtrackC; chi2trackC = deltatrackC*mat2trackC; //chi2trackC.Print(); } // // Propagate ITSout to TPC inner wall // and calculate chi2 distance to track (InnerParams) // const Double_t kTPCRadius=85; const Double_t kStep=3; // clone track InnerParams has to be deleted Bool_t isOKtrackInnerC2 = kFALSE; AliExternalTrackParam *trackInnerC2 = new AliExternalTrackParam(*(track->GetInnerParam())); if (trackInnerC2) { isOKtrackInnerC2 = AliTracker::PropagateTrackToBxByBz(trackInnerC2,kTPCRadius,track->GetMass(),kStep,kFALSE); } Bool_t isOKouterITSc = kFALSE; AliExternalTrackParam *outerITSc = NULL; TMatrixD chi2OuterITS(1,1); if(esdFriend && !esdFriend->TestSkipBit()) { // propagate ITSout to TPC inner wall if(friendTrack) { outerITSc = NULL; if (friendTrack->GetITSOut()) outerITSc = new AliExternalTrackParam(*(friendTrack->GetITSOut())); if(outerITSc) { isOKouterITSc = AliTracker::PropagateTrackToBxByBz(outerITSc,kTPCRadius,track->GetMass(),kStep,kFALSE); isOKouterITSc = outerITSc->Rotate(trackInnerC2->GetAlpha()); isOKouterITSc = outerITSc->PropagateTo(trackInnerC2->GetX(),esdEvent->GetMagneticField()); // // calculate chi2 between outerITS and innerParams // cov without z-coordinate at the moment // TMatrixD deltaTouterITS(4,1); TMatrixD deltaouterITS(1,4); TMatrixD covarMouterITS(4,4); if(isOKtrackInnerC2 && isOKouterITSc) { Int_t kipar = 0; Int_t kjpar = 0; for (Int_t ipar=0; ipar<5; ipar++) { if(ipar!=1) { deltaTouterITS(kipar,0)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar]; deltaouterITS(0,kipar)=outerITSc->GetParameter()[ipar]-trackInnerC2->GetParameter()[ipar]; } kjpar=0; for (Int_t jpar=0; jpar<5; jpar++) { Int_t index=outerITSc->GetIndex(ipar,jpar); if(ipar !=1 || jpar!=1) { covarMouterITS(kipar,kjpar)=outerITSc->GetCovariance()[index]+trackInnerC2->GetCovariance()[index]; } if(jpar!=1) kjpar++; } if(ipar!=1) kipar++; } // chi2 distance ITSout and InnerParams TMatrixD covarMInvouterITS = covarMouterITS.Invert(); TMatrixD mat2outerITS = covarMInvouterITS*deltaTouterITS; chi2OuterITS = deltaouterITS*mat2outerITS; //chi2OuterITS.Print(); } } } } // // MC info // TParticle *particle=NULL, *particleTPC=NULL, *particleITS=NULL; TParticle *particleMother=NULL, *particleMotherTPC=NULL, *particleMotherITS=NULL; Int_t mech=-1, mechTPC=-1, mechITS=-1; Bool_t isPrim=kFALSE, isPrimTPC=kFALSE, isPrimITS=kFALSE; Bool_t isFromStrangess=kFALSE, isFromStrangessTPC=kFALSE, isFromStrangessITS=kFALSE; Bool_t isFromConversion=kFALSE, isFromConversionTPC=kFALSE, isFromConversionITS=kFALSE; Bool_t isFromMaterial=kFALSE, isFromMaterialTPC=kFALSE, isFromMaterialITS=kFALSE; AliTrackReference *refTPCIn = NULL; AliTrackReference *refTPCOut = NULL; AliTrackReference *refITS = NULL; AliTrackReference *refTRD = NULL; AliTrackReference *refTOF = NULL; AliTrackReference *refEMCAL = NULL; AliTrackReference *refPHOS = NULL; Int_t nrefTPC=0, nrefTRD=0, nrefTOF=0, nrefITS=0, nrefEMCAL=0, nrefPHOS=0; Bool_t isOKtrackInnerC3 = kFALSE; AliExternalTrackParam *trackInnerC3 = new AliExternalTrackParam(*(track->GetInnerParam())); if(mcEvent && stack) { do //artificial loop (once) to make the continue statements jump out of the MC part { multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); // // global track // Int_t label = TMath::Abs(track->GetLabel()); if (label >= mcStackSize) continue; particle = stack->Particle(label); if (!particle) continue; if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0.) { particleMother = GetMother(particle,stack); mech = particle->GetUniqueID(); isPrim = stack->IsPhysicalPrimary(label); isFromStrangess = IsFromStrangeness(label,stack); isFromConversion = IsFromConversion(label,stack); isFromMaterial = IsFromMaterial(label,stack); } // // TPC track // Int_t labelTPC = TMath::Abs(track->GetTPCLabel()); if (labelTPC >= mcStackSize) continue; particleTPC = stack->Particle(labelTPC); if (!particleTPC) continue; if(particleTPC && particleTPC->GetPDG() && particleTPC->GetPDG()->Charge()!=0.) { particleMotherTPC = GetMother(particleTPC,stack); mechTPC = particleTPC->GetUniqueID(); isPrimTPC = stack->IsPhysicalPrimary(labelTPC); isFromStrangessTPC = IsFromStrangeness(labelTPC,stack); isFromConversionTPC = IsFromConversion(labelTPC,stack); isFromMaterialTPC = IsFromMaterial(labelTPC,stack); } // // store first track reference // for TPC track // TParticle *part=0; TClonesArray *trefs=0; Int_t status = mcEvent->GetParticleAndTR(TMath::Abs(labelTPC), part, trefs); if(status>0 && part && trefs && part->GetPDG() && part->GetPDG()->Charge()!=0.) { Int_t nTrackRef = trefs->GetEntries(); //printf("nTrackRef %d \n",nTrackRef); Int_t countITS = 0; for (Int_t iref = 0; iref < nTrackRef; iref++) { AliTrackReference *ref = (AliTrackReference *)trefs->At(iref); // Ref. in the middle ITS if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kITS) { nrefITS++; if(!refITS && countITS==2) { refITS = ref; //printf("refITS %p \n",refITS); } countITS++; } // TPC if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTPC) { nrefTPC++; refTPCOut=ref; if(!refTPCIn) { refTPCIn = ref; //printf("refTPCIn %p \n",refTPCIn); //break; } } // TRD if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTRD) { nrefTRD++; if(!refTRD) { refTRD = ref; } } // TOF if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kTOF) { nrefTOF++; if(!refTOF) { refTOF = ref; } } // EMCAL if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kEMCAL) { nrefEMCAL++; if(!refEMCAL) { refEMCAL = ref; } } // PHOS // if(ref && ref->Label()==label && ref->DetectorId()==AliTrackReference::kPHOS) // { // nrefPHOS++; // if(!refPHOS) { // refPHOS = ref; // } // } } // transform inner params to TrackRef // reference frame if(refTPCIn && trackInnerC3) { Double_t kRefPhi = TMath::ATan2(refTPCIn->Y(),refTPCIn->X()); isOKtrackInnerC3 = trackInnerC3->Rotate(kRefPhi); isOKtrackInnerC3 = AliTracker::PropagateTrackToBxByBz(trackInnerC3,refTPCIn->R(),track->GetMass(),kStep,kFALSE); } } // // ITS track // Int_t labelITS = TMath::Abs(track->GetITSLabel()); if (labelITS >= mcStackSize) continue; particleITS = stack->Particle(labelITS); if (!particleITS) continue; if(particleITS && particleITS->GetPDG() && particleITS->GetPDG()->Charge()!=0.) { particleMotherITS = GetMother(particleITS,stack); mechITS = particleITS->GetUniqueID(); isPrimITS = stack->IsPhysicalPrimary(labelITS); isFromStrangessITS = IsFromStrangeness(labelITS,stack); isFromConversionITS = IsFromConversion(labelITS,stack); isFromMaterialITS = IsFromMaterial(labelITS,stack); } } while (0); } // Bool_t dumpToTree=kFALSE; if(isOKtpcInnerC && isOKtrackInnerC) dumpToTree = kTRUE; //if(fUseESDfriends && isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; if(isOKtrackInnerC2 && isOKouterITSc) dumpToTree = kTRUE; if(mcEvent && isOKtrackInnerC3) dumpToTree = kTRUE; TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); if (fReducePileUp){ // // 18.03 - Reduce pile-up chunks, done outside of the ESDTrackCuts for 2012/2013 data pile-up about 95 % of tracks // Done only in case no MC info // Float_t dcaTPC[2]; track->GetImpactParametersTPC(dcaTPC[0],dcaTPC[1]); Bool_t isRoughPrimary = TMath::Abs(dcaTPC[1])<10; Bool_t hasOuter=(track->IsOn(AliVTrack::kITSin))||(track->IsOn(AliVTrack::kTOFout))||(track->IsOn(AliVTrack::kTRDin)); Bool_t keepPileUp=gRandom->Rndm()<0.05; if ( (!hasOuter) && (!isRoughPrimary) && (!keepPileUp)){ dumpToTree=kFALSE; } } //init dummy objects static AliESDVertex dummyvtxESD; //if (!dummyvtxESD) //{ // dummyvtxESD=new AliESDVertex(); // //dummyvtxESD->SetNContributors(1); // //UShort_t pindices[1]; pindices[0]=0; // //dummyvtxESD->SetIndices(1,pindices); // //dummyvtxESD->SetNContributors(0); //} static AliExternalTrackParam dummyexternaltrackparam; //if (!dummyexternaltrackparam) dummyexternaltrackparam=new AliExternalTrackParam(); static AliTrackReference dummytrackreference; //if (!dummytrackreference) dummytrackreference=new AliTrackReference(); static TParticle dummyparticle; //if (!dummyparticle) dummyparticle=new TParticle(); //assign the dummy objects if needed if (!track) {track=fDummyTrack;} AliESDfriendTrack *friendTrackStore=friendTrack; // store friend track for later processing if (fFriendDownscaling>=1){ // downscaling number of friend tracks friendTrackStore = (gRandom->Rndm()<1./fFriendDownscaling)? friendTrack:0; } if (fFriendDownscaling<=0){ if (((*fTreeSRedirector)<<"highPt").GetTree()){ TTree * tree = ((*fTreeSRedirector)<<"highPt").GetTree(); if (tree){ Double_t sizeAll=tree->GetZipBytes(); TBranch * br= tree->GetBranch("friendTrack.fPoints"); Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; br= tree->GetBranch("friendTrack.fCalibContainer"); if (br) sizeFriend+=br->GetZipBytes(); if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) friendTrackStore=0; } } } // if (!friendTrackStore && fFriendDownscaling<=1) {friendTrack=fDummyFriendTrack;} if (!vtxESD) {vtxESD=&dummyvtxESD;} if (mcEvent) { if (!refTPCIn) {refTPCIn=&dummytrackreference;} if (!refITS) {refITS=&dummytrackreference;} if (!particle) {particle=&dummyparticle;} if (!particleMother) {particleMother=&dummyparticle;} if (!particleTPC) {particleTPC=&dummyparticle;} if (!particleMotherTPC) {particleMotherTPC=&dummyparticle;} if (!particleITS) {particleITS=&dummyparticle;} if (!particleMotherITS) {particleMotherITS=&dummyparticle;} } Int_t ntracks = esdEvent->GetNumberOfTracks(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); // fill histograms FillHistograms(track, tpcInnerC, centralityF, (Double_t)chi2(0,0)); if(fTreeSRedirector && dumpToTree && fFillTree) { (*fTreeSRedirector)<<"highPt"<< "gid="<GetInputEventHandler(); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // use MC information AliHeader* header = 0; AliGenEventHeader* genHeader = 0; AliStack* stack = 0; Int_t mcStackSize=0; TArrayF vtxMC(3); Int_t multMCTrueTracks = 0; // if(!mcEvent) { AliDebug(AliLog::kError, "mcEvent not available"); return; } // get MC event header header = mcEvent->Header(); if (!header) { AliDebug(AliLog::kError, "Header not available"); return; } // MC particle stack stack = mcEvent->Stack(); if (!stack) { AliDebug(AliLog::kError, "Stack not available"); return; } mcStackSize=stack->GetNtrack(); // get MC vertex genHeader = header->GenEventHeader(); if (!genHeader) { AliDebug(AliLog::kError, "Could not retrieve genHeader from Header"); return; } genHeader->PrimaryVertex(vtxMC); // multipliticy of all MC primary tracks // in Zv, pt and eta ranges) multMCTrueTracks = GetMCTrueTrackMult(mcEvent,evtCuts,accCuts); // get reconstructed vertex //const AliESDVertex* vtxESD = 0; AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } if(!vtxESD) return; // // Vertex info comparison and track multiplicity // AliESDVertex *vertexSPD = (AliESDVertex *)esdEvent->GetPrimaryVertexSPD(); AliESDVertex *vertexTPC = (AliESDVertex *)esdEvent->GetPrimaryVertexTPC(); Int_t contSPD = vertexSPD->GetNContributors(); Int_t contTPC = vertexTPC->GetNContributors(); TVectorD vertexPosTPC(3), vertexPosSPD(3); vertexSPD->GetXYZ(vertexPosSPD.GetMatrixArray()); vertexTPC->GetXYZ(vertexPosTPC.GetMatrixArray()); Int_t ntracksTPC=0; Int_t ntracksITS=0; for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++){ AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; if (track->IsOn(AliVTrack::kTPCrefit)) ntracksTPC++; if (track->IsOn(AliVTrack::kITSrefit)) ntracksITS++; } Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); // check event cuts if(isEventOK && isEventTriggered) { //if(!stack) return; // // MC info // TParticle *particle=NULL; TParticle *particleMother=NULL; Int_t mech=-1; // reco event info Double_t vert[3] = {0}; vert[0] = vtxESD->GetX(); vert[1] = vtxESD->GetY(); vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); Double_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); // loop over MC stack for (Int_t iMc = 0; iMc < mcStackSize; ++iMc) { particle = stack->Particle(iMc); if (!particle) continue; // only charged particles if(!particle->GetPDG()) continue; Double_t charge = particle->GetPDG()->Charge()/3.; if (TMath::Abs(charge) < 0.001) continue; // only primary particles Bool_t prim = stack->IsPhysicalPrimary(iMc); if(!prim) continue; // downscale low-pT particles Double_t scalempt= TMath::Min(particle->Pt(),10.); Double_t downscaleF = gRandom->Rndm(); downscaleF *= fLowPtTrackDownscaligF; if(TMath::Exp(2*scalempt)AcceptTrack(particle)) continue; // check if particle reconstructed Bool_t isRec = kFALSE; Int_t trackIndex = -1; Int_t trackLoopIndex = -1; Int_t isESDtrackCut= 0; Int_t isAccCuts = 0; Int_t nRec = 0; // how many times reconstructed Int_t nFakes = 0; // how many times reconstructed as a fake track AliESDtrack *recTrack = NULL; for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; if(track->Charge()==0) continue; // Int_t label = TMath::Abs(track->GetLabel()); if (label >= mcStackSize) continue; if(label == iMc) { Bool_t isAcc=esdTrackCuts->AcceptTrack(track); if (isAcc) isESDtrackCut=1; if (accCuts->AcceptTrack(track)) isAccCuts=1; isRec = kTRUE; trackIndex = iTrack; if (recTrack){ if (track->GetTPCncls()GetTPCncls()) continue; // in case looper tracks use longer track if (!isAcc) continue; trackLoopIndex = iTrack; } recTrack = esdEvent->GetTrack(trackIndex); nRec++; if(track->GetLabel()<0) nFakes++; continue; } } // Store information in the output tree if (trackLoopIndex>-1) { recTrack = esdEvent->GetTrack(trackLoopIndex); } else if (trackIndex >-1) { recTrack = esdEvent->GetTrack(trackIndex); } else { recTrack = fDummyTrack; } particleMother = GetMother(particle,stack); mech = particle->GetUniqueID(); //MC particle track length Double_t tpcTrackLength = 0.; AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(iMc); if(mcParticle) { Int_t counter; tpcTrackLength = mcParticle->GetTPCTrackLength(bz,0.05,counter,3.0); } // if(fTreeSRedirector && fFillTree) { (*fTreeSRedirector)<<"MCEffTree"<< "fileName.="<<&fCurrentFileName<< "triggerClass.="<<&triggerClass<< "runNumber="< 1 // if (track->GetTPCNcls() < 60) return kFALSE; Double_t mom = track->GetInnerParam()->GetP(); if (mom < 0.2) return kFALSE; // protection against unexpected behavior of Aleph parameterization Float_t dca[2], bCov[3]; track->GetImpactParameters(dca,bCov); // Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541); if (track->GetTPCsignal() > triggerDeDx && track->GetTPCsignal()<1000 && TMath::Abs(dca[0])<3.) return kTRUE; return kFALSE; } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::ProcessV0(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent, AliESDfriend *const esdFriend) { // // Select real events with V0 (K0s and Lambda and Gamma) high-pT candidates // if(!esdEvent) { AliDebug(AliLog::kError, "esdEvent not available"); return; } // get selection cuts AliFilteredTreeEventCuts *evtCuts = GetEventCuts(); AliFilteredTreeAcceptanceCuts *accCuts = GetAcceptanceCuts(); AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); if(!evtCuts || !accCuts || !esdTrackCuts) { AliDebug(AliLog::kError, "cuts not available"); return; } // trigger selection Bool_t isEventTriggered = kTRUE; AliPhysicsSelection *physicsSelection = NULL; AliTriggerAnalysis* triggerAnalysis = NULL; // AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; //SetPhysicsTriggerSelection(physicsSelection); if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // get reconstructed vertex //const AliESDVertex* vtxESD = 0; AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } if(!vtxESD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); // check event cuts if(isEventOK && isEventTriggered) { // // Dump the pt downscaled V0 into the tree // Int_t ntracks = esdEvent->GetNumberOfTracks(); Int_t nV0s = esdEvent->GetNumberOfV0s(); Int_t run = esdEvent->GetRunNumber(); Int_t time= esdEvent->GetTimeStamp(); Int_t evNr=esdEvent->GetEventNumberInFile(); Double_t bz = esdEvent->GetMagneticField(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); for (Int_t iv0=0; iv0GetV0(iv0); if (!v0) continue; AliESDtrack * track0 = esdEvent->GetTrack(v0->GetIndex(0)); AliESDtrack * track1 = esdEvent->GetTrack(v0->GetIndex(1)); if (!track0) continue; if (!track1) continue; AliESDfriendTrack* friendTrack0=NULL; AliESDfriendTrack* friendTrack1=NULL; if (esdFriend) { if (!esdFriend->TestSkipBit()){ Int_t ntracksFriend = esdFriend->GetNumberOfTracks(); if (v0->GetIndex(0)GetTrack(v0->GetIndex(0)); //this guy can be NULL } if (v0->GetIndex(1)GetTrack(v0->GetIndex(1)); //this guy can be NULL } } } if (track0->GetSign()<0) { track1 = esdEvent->GetTrack(v0->GetIndex(0)); track0 = esdEvent->GetTrack(v0->GetIndex(1)); } // AliESDfriendTrack *friendTrackStore0=friendTrack0; // store friend track0 for later processing AliESDfriendTrack *friendTrackStore1=friendTrack1; // store friend track1 for later processing if (fFriendDownscaling>=1){ // downscaling number of friend tracks if (gRandom->Rndm()>1./fFriendDownscaling){ friendTrackStore0 = 0; friendTrackStore1 = 0; } } if (fFriendDownscaling<=0){ if (((*fTreeSRedirector)<<"V0s").GetTree()){ TTree * tree = ((*fTreeSRedirector)<<"V0s").GetTree(); if (tree){ Double_t sizeAll=tree->GetZipBytes(); TBranch * br= tree->GetBranch("friendTrack0.fPoints"); Double_t sizeFriend=(br!=NULL)?br->GetZipBytes():0; br= tree->GetBranch("friendTrack0.fCalibContainer"); if (br) sizeFriend+=br->GetZipBytes(); if (sizeFriend*TMath::Abs(fFriendDownscaling)>sizeAll) { friendTrackStore0=0; friendTrackStore1=0; } } } } // Bool_t isDownscaled = IsV0Downscaled(v0); if (isDownscaled) continue; AliKFParticle kfparticle; // Int_t type=GetKFParticle(v0,esdEvent,kfparticle); if (type==0) continue; TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"V0s"<< "gid="<GetInputEventHandler(); if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // get reconstructed vertex AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } if(!vtxESD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); // check event cuts if(isEventOK && isEventTriggered) { Double_t vert[3] = {0}; vert[0] = vtxESD->GetX(); vert[1] = vtxESD->GetY(); vert[2] = vtxESD->GetZ(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); Double_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); // Global event id calculation using orbitID, bunchCrossingID and periodID ULong64_t orbitID = (ULong64_t)esdEvent->GetOrbitNumber(); ULong64_t bunchCrossID = (ULong64_t)esdEvent->GetBunchCrossNumber(); ULong64_t periodID = (ULong64_t)esdEvent->GetPeriodNumber(); ULong64_t gid = ((periodID << 36) | (orbitID << 12) | bunchCrossID); // large dEdx for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) { AliESDtrack *track = esdEvent->GetTrack(iTrack); if(!track) continue; AliESDfriendTrack* friendTrack=NULL; if (esdFriend && !esdFriend->TestSkipBit()) { Int_t ntracksFriend = esdFriend->GetNumberOfTracks(); if (iTrackGetTrack(iTrack); } //this guy can be NULL } if(track->Charge()==0) continue; if(!esdTrackCuts->AcceptTrack(track)) continue; if(!accCuts->AcceptTrack(track)) continue; if(!IsHighDeDxParticle(track)) continue; TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"dEdx"<< // high dEdx tree "gid="<GetOnFlyStatus() ==kFALSE) return 0; // // 1.) track cut // AliESDtrack * track0 = event->GetTrack(v0->GetIndex(0)); AliESDtrack * track1 = event->GetTrack(v0->GetIndex(1)); /* TCut cutD="abs(track0.fD/sqrt(track0.fCdd))>2&&abs(track1.fD/sqrt(track1.fCdd))>2"; TCut cutTheta="abs(track0.fP[3])<1&&abs(track1.fP[3])<1"; TCut cutNcl="track0.GetTPCClusterInfo(2,1)>100&&track1.GetTPCClusterInfo(2,1)>100"; */ if (TMath::Abs(track0->GetTgl())>1) return 0; if (TMath::Abs(track1->GetTgl())>1) return 0; if ((track0->GetTPCClusterInfo(2,1))<100) return 0; if ((track1->GetTPCClusterInfo(2,1))<100) return 0; Float_t pos0[2]={0}, cov0[3]={0}; Float_t pos1[2]={0}, cov1[3]={0}; track0->GetImpactParameters(pos0,cov0); track0->GetImpactParameters(pos1,cov1); // if (TMath::Abs(pos0[0])GetKFInfo(2,2,2); if (chi2KF>25) return 0; // // 4.) Rough mass cut - 0.200 GeV // static Double_t masses[2]={-1}; if (masses[0]<0){ masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass(); masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass(); } Double_t mass00= v0->GetEffMass(0,0); Double_t mass22= v0->GetEffMass(2,2); Double_t mass42= v0->GetEffMass(4,2); Double_t mass24= v0->GetEffMass(2,4); Bool_t massOK=kFALSE; Int_t type=0; Int_t ptype=0; Double_t dmass=1; Int_t p1=0, p2=0; if (TMath::Abs(mass00-0)GetParamP(); const AliExternalTrackParam *paramN = v0->GetParamN(); if (paramP->GetSign()<0){ paramP=v0->GetParamP(); paramN=v0->GetParamN(); } //Double_t *pparam1 = (Double_t*)paramP->GetParameter(); //Double_t *pparam2 = (Double_t*)paramN->GetParameter(); // AliKFParticle kfp1( *paramP, spdg[p1] ); AliKFParticle kfp2( *paramN, -1 *spdg[p2] ); AliKFParticle V0KF; (V0KF)+=kfp1; (V0KF)+=kfp2; kfparticle=V0KF; // // Pointing angle // Double_t errPhi = V0KF.GetErrPhi(); Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle()); if (pointAngle/errPhi>10) return 0; // return ptype; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0) { // // Downscale randomly low pt V0 // //return kFALSE; Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt()); Double_t scalempt= TMath::Min(maxPt,10.); Double_t downscaleF = gRandom->Rndm(); downscaleF *= fLowPtV0DownscaligF; // // Special treatment of the gamma conversion pt spectra is softer - Double_t mass00= v0->GetEffMass(0,0); const Double_t cutMass=0.2; if (TMath::Abs(mass00-0)Exp(1); Bool_t isDownscaled =TMath::Exp(rnd)<100*gRandom->Rndm(); his1->Fill(rnd); if (!isDownscaled) his2->Fill(rnd); }} */ } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::ConstrainTPCInner(AliExternalTrackParam *const tpcInnerC, const AliESDVertex* vtx, Double_t b[3]) { // Constrain TPC inner params constrained // if(!tpcInnerC) return kFALSE; if(!vtx) return kFALSE; Double_t dz[2],cov[3]; //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex(); //if(!tpcInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; //if(!tpcInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; if(!tpcInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; Double_t covar[6]; vtx->GetCovMatrix(covar); Double_t p[2]={tpcInnerC->GetParameter()[0]-dz[0],tpcInnerC->GetParameter()[1]-dz[1]}; Double_t c[3]={covar[2],0.,covar[5]}; Double_t chi2C=tpcInnerC->GetPredictedChi2(p,c); if (chi2C>kVeryBig) return kFALSE; if(!tpcInnerC->Update(p,c)) return kFALSE; return kTRUE; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::ConstrainTrackInner(AliExternalTrackParam *const trackInnerC, const AliESDVertex* vtx, Double_t mass, Double_t b[3]) { // Constrain TPC inner params constrained // if(!trackInnerC) return kFALSE; if(!vtx) return kFALSE; const Double_t kRadius = 2.8; const Double_t kMaxStep = 1.0; Double_t dz[2],cov[3]; //AliESDVertex *vtx= (AliESDVertex *)esdEvent->GetPrimaryVertex(); //if(!trackInnerC->PropagateToDCA(vtx, esdEvent->GetMagneticField(), 3, dz, cov)) return kFALSE; //if(!trackInnerC->PropagateToDCA(vtx, Bz, 3, dz, cov)) return kFALSE; if(!AliTracker::PropagateTrackToBxByBz(trackInnerC,kRadius,mass,kMaxStep,kFALSE)) return kFALSE; if(!trackInnerC->PropagateToDCABxByBz(vtx, b, 3, dz, cov)) return kFALSE; // Double_t covar[6]; vtx->GetCovMatrix(covar); Double_t p[2]={trackInnerC->GetParameter()[0]-dz[0],trackInnerC->GetParameter()[1]-dz[1]}; Double_t c[3]={covar[2],0.,covar[5]}; Double_t chi2C=trackInnerC->GetPredictedChi2(p,c); if (chi2C>kVeryBig) return kFALSE; if(!trackInnerC->Update(p,c)) return kFALSE; return kTRUE; } //_____________________________________________________________________________ TParticle *AliAnalysisTaskFilteredTree::GetMother(TParticle *const particle, AliStack *const stack) { if(!particle) return NULL; if(!stack) return NULL; Int_t motherLabel = TMath::Abs(particle->GetMother(0)); TParticle* mother = NULL; Int_t mcStackSize=stack->GetNtrack(); if (motherLabel>=mcStackSize) return NULL; mother = stack->Particle(motherLabel); return mother; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::IsFromConversion(Int_t label, AliStack *const stack) { Bool_t isFromConversion = kFALSE; if(stack) { Int_t mcStackSize=stack->GetNtrack(); if (label>=mcStackSize) return kFALSE; TParticle* particle = stack->Particle(label); if (!particle) return kFALSE; if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) { Int_t mech = particle->GetUniqueID(); // production mechanism Bool_t isPrim = stack->IsPhysicalPrimary(label); Int_t motherLabel = TMath::Abs(particle->GetMother(0)); if (motherLabel>=mcStackSize) return kFALSE; TParticle* mother = stack->Particle(motherLabel); if(mother) { Int_t motherPdg = mother->GetPdgCode(); if(!isPrim && mech==5 && motherPdg==kGamma) { isFromConversion=kTRUE; } } } } return isFromConversion; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::IsFromMaterial(Int_t label, AliStack *const stack) { Bool_t isFromMaterial = kFALSE; if(stack) { Int_t mcStackSize=stack->GetNtrack(); if (label>=mcStackSize) return kFALSE; TParticle* particle = stack->Particle(label); if (!particle) return kFALSE; if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) { Int_t mech = particle->GetUniqueID(); // production mechanism Bool_t isPrim = stack->IsPhysicalPrimary(label); Int_t motherLabel = TMath::Abs(particle->GetMother(0)); if (motherLabel>=mcStackSize) return kFALSE; TParticle* mother = stack->Particle(motherLabel); if(mother) { if(!isPrim && mech==13) { isFromMaterial=kTRUE; } } } } return isFromMaterial; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::IsFromStrangeness(Int_t label, AliStack *const stack) { Bool_t isFromStrangeness = kFALSE; if(stack) { Int_t mcStackSize=stack->GetNtrack(); if (label>=mcStackSize) return kFALSE; TParticle* particle = stack->Particle(label); if (!particle) return kFALSE; if(particle && particle->GetPDG() && particle->GetPDG()->Charge()!=0) { Int_t mech = particle->GetUniqueID(); // production mechanism Bool_t isPrim = stack->IsPhysicalPrimary(label); Int_t motherLabel = TMath::Abs(particle->GetMother(0)); if (motherLabel>=mcStackSize) return kFALSE; TParticle* mother = stack->Particle(motherLabel); if(mother) { Int_t motherPdg = mother->GetPdgCode(); // K+-, lambda, antilambda, K0s decays if(!isPrim && mech==4 && (TMath::Abs(motherPdg)==kKPlus || TMath::Abs(motherPdg)==kLambda0 || motherPdg==kK0Short)) { isFromStrangeness = kTRUE; } } } } return isFromStrangeness; } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::FinishTaskOutput() { // // Called one at the end // locally on working node // Bool_t deleteTrees=kTRUE; if ((AliAnalysisManager::GetAnalysisManager())) { if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == AliAnalysisManager::kProofAnalysis) deleteTrees=kFALSE; } if (deleteTrees) delete fTreeSRedirector; fTreeSRedirector=NULL; } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::Terminate(Option_t *) { // // Called one at the end // } //_____________________________________________________________________________ Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts) { // // calculate mc event true track multiplicity // if(!mcEvent) return 0; AliStack* stack = 0; Int_t mult = 0; // MC particle stack stack = mcEvent->Stack(); if (!stack) return 0; // //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv()); // Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent); if(!isEventOK) return 0; Int_t nPart = stack->GetNtrack(); for (Int_t iMc = 0; iMc < nPart; ++iMc) { TParticle* particle = stack->Particle(iMc); if (!particle) continue; // only charged particles if(!particle->GetPDG()) continue; Double_t charge = particle->GetPDG()->Charge()/3.; if (TMath::Abs(charge) < 0.001) continue; // physical primary Bool_t prim = stack->IsPhysicalPrimary(iMc); if(!prim) continue; // checked accepted without pt cut //if(accCuts->AcceptTrack(particle)) if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) { mult++; } } return mult; } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::FillHistograms(AliESDtrack* const ptrack, AliExternalTrackParam* const ptpcInnerC, Double_t centralityF, Double_t chi2TPCInnerC) { // // Fill pT relative resolution histograms for // TPC only, TPC only constrained to vertex and TPC+ITS tracking // if(!ptrack) return; if(!ptpcInnerC) return; const AliExternalTrackParam * innerParam = (AliExternalTrackParam *) ptrack->GetInnerParam(); if(!innerParam) return; Float_t dxy, dz; ptrack->GetImpactParameters(dxy,dz); // TPC+ITS primary tracks if( abs(ptrack->Eta())<0.8 && ptrack->GetTPCClusterInfo(3,1)>120 && ptrack->IsOn(0x40) && ptrack->GetTPCclusters(0)>0.0 && ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && //abs(innerParam->GetX())>0.0 && //abs(innerParam->GetY()/innerParam->GetX())<0.14 && //abs(innerParam->GetTgl())<0.85 && ptrack->IsOn(0x0004) && ptrack->GetNcls(0)>0 && ptrack->GetITSchi2()>0 && sqrt(ptrack->GetITSchi2()/ptrack->GetNcls(0))<6 && sqrt(chi2TPCInnerC)<6 && (ptrack->HasPointOnITSLayer(0) || ptrack->HasPointOnITSLayer(1)) && abs(dz)<2.0 && abs(dxy)<(0.018+0.035*abs(ptrack->GetSigned1Pt())) ) { fPtResPhiPtTPCITS->Fill(ptrack->Pt(),ptrack->Phi(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2())); fPtResEtaPtTPCITS->Fill(ptrack->Pt(),ptrack->Eta(),1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2())); fPtResCentPtTPCITS->Fill(ptrack->Pt(),centralityF,1./abs(ptrack->GetSigned1Pt())*TMath::Sqrt(ptrack->GetSigma1Pt2())); } // TPC primary tracks // and TPC constrained primary tracks AliExternalTrackParam *ptpcInner = (AliExternalTrackParam *) ptrack->GetTPCInnerParam(); if(!ptpcInner) return; Float_t dxyTPC, dzTPC; ptrack->GetImpactParametersTPC(dxyTPC,dzTPC); if( abs(ptrack->Eta())<0.8 && ptrack->GetTPCClusterInfo(3,1)>120 && ptrack->IsOn(0x40)&& ptrack->GetTPCclusters(0)>0.0 && ptrack->GetTPCnclsS()/ptrack->GetTPCclusters(0)<0.4 && //abs(innerParam->GetX())>0.0 && //abs(innerParam->GetY()/innerParam->GetX())<0.14 && //abs(innerParam->GetTgl())<0.85 && abs(dzTPC)<3.2 && abs(dxyTPC)<2.4 ) { // TPC only fPtResPhiPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Phi(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2())); fPtResEtaPtTPC->Fill(ptpcInner->Pt(),ptpcInner->Eta(),1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2())); fPtResCentPtTPC->Fill(ptpcInner->Pt(),centralityF,1./abs(ptpcInner->GetSigned1Pt())*TMath::Sqrt(ptpcInner->GetSigma1Pt2())); // TPC constrained to vertex fPtResPhiPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Phi(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2())); fPtResEtaPtTPCc->Fill(ptpcInnerC->Pt(),ptpcInnerC->Eta(),1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2())); fPtResCentPtTPCc->Fill(ptpcInnerC->Pt(),centralityF,1./abs(ptpcInnerC->GetSigned1Pt())*TMath::Sqrt(ptpcInnerC->GetSigma1Pt2())); } } void AliAnalysisTaskFilteredTree::ProcessITSTPCmatchOut(AliESDEvent *const esdEvent, AliESDfriend *const esdFriend){ // // Process ITS standalone tracks find match with closest TPC(or combined tracks) tracks // marian.ivanov@cern.ch // 0.) Init variables // 1.) GetTrack parameters at TPC inner wall // 2.) Match closest TPC track (STANDALONE/global) - chi2 match criteria // // Logic to be used in reco: // 1.) Find matching ITSalone->TPCalone // 2.) if (!TPCalone.FindClose(TPCother)) TPCalone.Addopt(ITSalone) // 3.) ff ((ITSalone.FindClose(Global)==0) CreateGlobaltrack const Double_t radiusMatch=84.; // redius to propagate // const Double_t dFastPhiCut=0.2; // 6 sigma (200 MeV) fast angular cut const Double_t dFastThetaCut=0.12; // 6 sigma (200 MeV) fast angular cut const Double_t dFastPosPhiCut=0.06; // 6 sigma (200 MeV) fast angular cut const Double_t dFastZCut=6; // 6 sigma (200 MeV) fast z difference cut const Double_t dFastPtCut=2.; // 6 sigma (200 MeV) fast 1/pt cut const Double_t chi2Cut=100; // chi2 matching cut // if (!esdFriend) return; // not ITS standalone track if (esdFriend->TestSkipBit()) return; // friends tracks not stored Int_t ntracks=esdEvent->GetNumberOfTracks(); Float_t bz = esdEvent->GetMagneticField(); // // 0.) Get parameters in reference radius TPC Inner wall // // TMatrixD vecPosR0(ntracks,6); // possition and momentum estimate at reference radius TMatrixD vecMomR0(ntracks,6); // TMatrixD vecPosR1(ntracks,6); // possition and momentum estimate at reference radius TPC track TMatrixD vecMomR1(ntracks,6); // Double_t xyz[3], pxyz[3]; // for (Int_t iTrack=0; iTrackGetTrack(iTrack); if(!track) continue; if (track->GetInnerParam()){ const AliExternalTrackParam *trackTPC=track->GetInnerParam(); trackTPC->GetXYZAt(radiusMatch,bz,xyz); trackTPC->GetPxPyPzAt(radiusMatch,bz,pxyz); for (Int_t i=0; i<3; i++){ vecPosR1(iTrack,i)=xyz[i]; vecMomR1(iTrack,i)=pxyz[i]; } vecPosR1(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); // phi pos angle vecMomR1(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); // phi mom angle vecMomR1(iTrack,4)= trackTPC->GetSigned1Pt();; vecMomR1(iTrack,5)= trackTPC->GetTgl();; } AliESDfriendTrack* friendTrack=esdFriend->GetTrack(iTrack); if(!friendTrack) continue; if (friendTrack->GetITSOut()){ const AliExternalTrackParam *trackITS=friendTrack->GetITSOut(); trackITS->GetXYZAt(radiusMatch,bz,xyz); trackITS->GetPxPyPzAt(radiusMatch,bz,pxyz); for (Int_t i=0; i<3; i++){ vecPosR0(iTrack,i)=xyz[i]; vecMomR0(iTrack,i)=pxyz[i]; } vecPosR0(iTrack,3)= TMath::ATan2(xyz[1],xyz[0]); vecMomR0(iTrack,3)= TMath::ATan2(pxyz[1],pxyz[0]); vecMomR0(iTrack,4)= trackITS->GetSigned1Pt();; vecMomR0(iTrack,5)= trackITS->GetTgl();; } } // // 1.) Find closest matching tracks, between the ITS standalone track // and the all other tracks // a.) caltegory - All // b.) category - without ITS // // Int_t ntracksPropagated=0; AliExternalTrackParam extTrackDummy; AliESDtrack esdTrackDummy; AliExternalTrackParam itsAtTPC; AliExternalTrackParam itsAtITSTPC; for (Int_t iTrack0=0; iTrack0GetTrack(iTrack0); if(!track0) continue; if (track0->IsOn(AliVTrack::kTPCin)) continue; AliESDfriendTrack* friendTrack0=esdFriend->GetTrack(iTrack0); if (!friendTrack0) continue; //if (!track0->IsOn(AliVTrack::kITSpureSA)) continue; //if (!friendTrack0->GetITSOut()) continue; // is there flag for ITS standalone? ntracksPropagated++; // // 2.) find clostest TPCtrack // a.) all tracks Double_t minChi2All=10000000; Double_t minChi2TPC=10000000; Double_t minChi2TPCITS=10000000; Int_t indexAll=-1; Int_t indexTPC=-1; Int_t indexTPCITS=-1; Int_t ncandidates0=0; // n candidates - rough cut Int_t ncandidates1=0; // n candidates - rough + chi2 cut itsAtTPC=*(friendTrack0->GetITSOut()); itsAtITSTPC=*(friendTrack0->GetITSOut()); for (Int_t iTrack1=0; iTrack1GetTrack(iTrack1); if(!track1) continue; if (!track1->IsOn(AliVTrack::kTPCin)) continue; // fast checks // if (TMath::Abs(vecPosR1(iTrack1,2)-vecPosR0(iTrack0,2))>dFastZCut) continue; if (TMath::Abs(vecPosR1(iTrack1,3)-vecPosR0(iTrack0,3))>dFastPosPhiCut) continue; if (TMath::Abs(vecMomR1(iTrack1,3)-vecMomR0(iTrack0,3))>dFastPhiCut) continue; if (TMath::Abs(vecMomR1(iTrack1,5)-vecMomR0(iTrack0,5))>dFastThetaCut) continue; if (TMath::Abs(vecMomR1(iTrack1,4)-vecMomR0(iTrack0,4))>dFastPtCut) continue; ncandidates0++; // const AliExternalTrackParam * param1= track1->GetInnerParam(); if (!friendTrack0->GetITSOut()) continue; AliExternalTrackParam outerITS = *(friendTrack0->GetITSOut()); if (!outerITS.Rotate(param1->GetAlpha())) continue; if (!outerITS.PropagateTo(param1->GetX(),bz)) continue; // assume track close to the TPC inner wall Double_t chi2 = outerITS.GetPredictedChi2(param1); if (chi2>chi2Cut) continue; ncandidates1++; if (chi2IsOn(AliVTrack::kITSin)==0){ minChi2TPC=chi2; indexTPC=iTrack1; itsAtTPC=outerITS; } if (chi2IsOn(AliVTrack::kITSin)){ minChi2TPCITS=chi2; indexTPCITS=iTrack1; itsAtITSTPC=outerITS; } } // AliESDtrack * trackAll= (indexAll>=0)? esdEvent->GetTrack(indexAll):&esdTrackDummy; AliESDtrack * trackTPC= (indexTPC>=0)? esdEvent->GetTrack(indexTPC):&esdTrackDummy; AliESDtrack * trackTPCITS= (indexTPCITS>=0)? esdEvent->GetTrack(indexTPCITS):&esdTrackDummy; (*fTreeSRedirector)<<"itsTPC"<< "indexAll="<