/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ #include "iostream" #include "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) , fProcessAll(kFALSE) , fProcessCosmics(kFALSE) , 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) , fDummyFriendTrack(0) , 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() { //the output trees not to be deleted in case of proof analysis //Bool_t deleteTrees=kTRUE; //if ((AliAnalysisManager::GetAnalysisManager())) //{ // if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == // AliAnalysisManager::kProofAnalysis) // deleteTrees=kFALSE; //} //if (deleteTrees) delete fTreeSRedirector; 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()); } 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 (!fDummyFriendTrack) { fDummyFriendTrack=new AliESDfriendTrack(); fDummyFriendTrack->SetTrackPointArray(new AliTrackPointArray(164)); printf("just made a new dummy friend track!"); } 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; } //// MC event //if(fUseMCInfo) { // fMC = MCEvent(); // if (!fMC) { // Printf("ERROR: MC 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"); } } //if set, use the environment variables to set the downscaling factors //AliAnalysisTaskFilteredTree_fLowPtTrackDownscaligF //AliAnalysisTaskFilteredTree_fLowPtV0DownscaligF 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)); } // 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);} printf("processed event %i\n", Entry()); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::ProcessCosmics(AliESDEvent *const event, AliESDfriend* esdFriend) { // // Select real events with high-pT tracks // if(!event) { AliDebug(AliLog::kError, "event not available"); return; } // AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } // get file name TTree *chain = (TChain*)GetInputData(0); if(!chain) { Printf("ERROR: Could not receive input chain"); return; } TObjString fileName(chain->GetCurrentFile()->GetName()); // check for cosmic pairs // // find cosmic pairs trigger by random trigger // // AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD(); AliESDVertex *vertexTPC = (AliESDVertex *)event->GetPrimaryVertexTPC(); const Double_t kMinPt=0.8; const Double_t kMinPtMax=0.8; const Double_t kMinNcl=50; const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1}; Int_t ntracks=event->GetNumberOfTracks(); // Float_t dcaTPC[2]={0,0}; // Float_t covTPC[3]={0,0,0}; UInt_t specie = event->GetEventSpecie(); // skip laser events if (specie==AliRecoParam::kCalib) return; for (Int_t itrack0=0;itrack0GetTrack(itrack0); if (!track0) continue; if (!track0->IsOn(AliESDtrack::kTPCrefit)) continue; if (TMath::Abs(AliTracker::GetBz())>1 && track0->Pt() < kMinPt) continue; if (track0->Pt() < kMinPt) continue; if (track0->GetTPCncls() < kMinNcl) continue; if (TMath::Abs(track0->GetY())GetKinkIndex(0)>0) continue; const Double_t * par0=track0->GetParameter(); //track param at rhe DCA //rm primaries // //track0->GetImpactParametersTPC(dcaTPC,covTPC); //if (TMath::Abs(dcaTPC[0])GetInnerParam(); AliESDfriendTrack* friendTrack0=NULL; if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack0 = esdFriend->GetTrack(itrack0);} //this guy can be NULL //Bool_t newFriendTrack0=kFALSE; //if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;} if (!friendTrack0) {friendTrack0=fDummyFriendTrack;} for (Int_t itrack1=itrack0+1;itrack1GetTrack(itrack1); if (!track1) continue; if (!track1->IsOn(AliESDtrack::kTPCrefit)) continue; if (track1->GetKinkIndex(0)>0) continue; if ((TMath::Abs(AliTracker::GetBz())>1) && (track1->Pt() < kMinPt)) continue; if (track1->Pt() < kMinPt) continue; if (track1->GetTPCncls()1 && TMath::Max(track1->Pt(), track0->Pt())GetY())GetImpactParametersTPC(dcaTPC,covTPC); // if (TMath::Abs(dcaTPC[0])GetParameter(); //track param at rhe DCA // Bool_t isPair=kTRUE; for (Int_t ipar=0; ipar<5; ipar++){ if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE; } if (!isPair) continue; if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE; //delta with correct sign /* TCut cut0="abs(t1.fP[0]+t0.fP[0])<2" TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02" TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2" */ if (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y opposite sign if (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign if (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign if (!isPair) continue; TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); Int_t eventNumber = event->GetEventNumberInFile(); //Bool_t hasFriend = kFALSE; //Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4); //printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS); // const AliExternalTrackParam * trackIn1 = track1->GetInnerParam(); // // Int_t ntracksSPD = vertexSPD->GetNContributors(); Int_t ntracksTPC = vertexTPC->GetNContributors(); Int_t runNumber = event->GetRunNumber(); Int_t timeStamp = event->GetTimeStamp(); ULong64_t triggerMask = event->GetTriggerMask(); Float_t magField = event->GetMagneticField(); TObjString triggerClass = event->GetFiredTriggerClasses().Data(); // 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); // // Dump to the tree // vertex // TPC-ITS tracks // //fCosmicPairsTree->Branch("fileName",&fileName,32000,0); //fCosmicPairsTree->Branch("runNumber",&runNumber,"runNumber/I"); //fCosmicPairsTree->Branch("timeStamp",&timeStamp,"timeStamp/I"); //fCosmicPairsTree->Branch("eventNumber",&eventNumber,"eventNumber/I"); //fCosmicPairsTree->Branch("triggerMask",&triggerMask,32000,0); //fCosmicPairsTree->Branch("triggerClass",&triggerClass,32000,0); //fCosmicPairsTree->Branch("magField",&magField,"magField/F"); //fCosmicPairsTree->Branch("ntracksSPD",&ntracksSPD,"ntracksSPD/I"); //fCosmicPairsTree->Branch("ntracksTPC",&ntracksTPC,"ntracksTPC/I"); //fCosmicPairsTree->Branch("vertexSPD",vertexSPD,32000,0); //fCosmicPairsTree->Branch("vertexTPC",vertexTPC,32000,0); //fCosmicPairsTree->Branch("track0",track0,32000,0); //fCosmicPairsTree->Branch("track1",track1,32000,0); // //fCosmicPairsTree->Fill(); AliESDfriendTrack* friendTrack1=NULL; if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack1 = esdFriend->GetTrack(itrack1);} //this guy can be NULL //Bool_t newFriendTrack1=kFALSE; //if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;} if (!friendTrack1) {friendTrack1=fDummyFriendTrack;} if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"CosmicPairs"<< "gid="<GetInputEventHandler(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } // get file name TTree *chain = (TChain*)GetInputData(0); if(!chain) { Printf("ERROR: Could not receive input chain"); return; } TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; //SetPhysicsTriggerSelection(physicsSelection); if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // use MC information AliHeader* header = 0; AliGenEventHeader* genHeader = 0; AliStack* stack = 0; TArrayF vtxMC(3); Int_t multMCTrueTracks = 0; if(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->GetZv()); //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->GetXv(); //vert[1] = vtxESD->GetYv(); //vert[2] = vtxESD->GetZv(); 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(); //fHighPtTree->Branch("fileName",&fileName,32000,0); //fHighPtTree->Branch("runNumber",&runNumber,"runNumber/I"); //fHighPtTree->Branch("evtTimeStamp",&evtTimeStamp,"evtTimeStamp/I"); //fHighPtTree->Branch("evtNumberInFile",&evtNumberInFile,"evtNumberInFile/I"); //fHighPtTree->Branch("triggerClass",&triggerClass,32000,0); //fHighPtTree->Branch("Bz",&bz,"Bz/F"); //fHighPtTree->Branch("vtxESD",vtxESD,32000,0); //fHighPtTree->Branch("IRtot",&ir1,"IRtot/I"); //fHighPtTree->Branch("IRint2",&ir2,"IRint2/I"); //fHighPtTree->Branch("mult",&mult,"mult/I"); //fHighPtTree->Branch("esdTrack",track,32000,0); //fHighPtTree->Branch("centralityF",¢ralityF,"centralityF/F"); //fHighPtTree->Fill(); //Double_t vtxX=vtxESD->GetX(); //Double_t vtxY=vtxESD->GetY(); //Double_t vtxZ=vtxESD->GetZ(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"highPt"<< "gid="<GetCurrentFile()->GetName()); if(!fFillTree) return; if(!fTreeSRedirector) return; // laser events 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; if (esdFriend) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL //Bool_t newFriendTrack=kFALSE; //if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;} if (!friendTrack) {friendTrack=fDummyFriendTrack;} (*fTreeSRedirector)<<"Laser"<< "gid="<GetInputEventHandler(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } // get file name TTree *chain = (TChain*)GetInputData(0); if(!chain) { Printf("ERROR: Could not receive input chain"); return; } TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) {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->GetXv(); vert[1] = vtxESD->GetYv(); vert[2] = vtxESD->GetZv(); 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;} if (!friendTrack) {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;} } //the following are guaranteed to exist: //if (!tpcInnerC) {tpcInnerC=&dummyexternaltrackparam;} //if (!trackInnerC) {trackInnerC=&dummyexternaltrackparam;} //if (!trackInnerC2) {trackInnerC2=&dummyexternaltrackparam;} //if (!outerITSc) {outerITSc=&dummyexternaltrackparam;} //if (!trackInnerC3) {trackInnerC3=&dummyexternaltrackparam;} ///////////////////////// //Double_t vtxX=vtxESD->GetX(); //Double_t vtxY=vtxESD->GetY(); //Double_t vtxZ=vtxESD->GetZ(); //AliESDVertex* pvtxESD = (AliESDVertex*)vtxESD->Clone(); //AliESDtrack* ptrack=(AliESDtrack*)track->Clone(); //AliExternalTrackParam* ptpcInnerC = (AliExternalTrackParam*)tpcInnerC->Clone(); //AliExternalTrackParam* ptrackInnerC = (AliExternalTrackParam*)trackInnerC->Clone(); //AliExternalTrackParam* ptrackInnerC2 = (AliExternalTrackParam*)trackInnerC2->Clone(); //AliExternalTrackParam* pouterITSc = (AliExternalTrackParam*)outerITSc->Clone(); //AliExternalTrackParam* ptrackInnerC3 = (AliExternalTrackParam*)trackInnerC3->Clone(); //AliESDVertex* pvtxESD = new AliESDVertex(*vtxESD); //AliESDtrack* ptrack= new AliESDtrack(*track); //AliExternalTrackParam* ptpcInnerC = new AliExternalTrackParam(*tpcInnerC); //AliExternalTrackParam* ptrackInnerC = new AliExternalTrackParam(*trackInnerC); //AliExternalTrackParam* ptrackInnerC2 = new AliExternalTrackParam(*trackInnerC2); //AliExternalTrackParam* pouterITSc = new AliExternalTrackParam(*outerITSc); //AliExternalTrackParam* ptrackInnerC3 = new AliExternalTrackParam(*trackInnerC3); 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) { //if (friendTrack) //{ // const AliTrackPointArray* array = friendTrack->GetTrackPointArray(); // if (!array) {printf("we have a friend, but the ponits are empty\n"); continue;} // if (friendTrack==fDummyFriendTrack) printf("using the dummy friend track\n"); // cout<GetX()[0]<<" "<GetX()[2]<GetInputEventHandler(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } // get file name TTree *chain = (TChain*)GetInputData(0); if(!chain) { Printf("ERROR: Could not receive input chain"); return; } TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // use MC information AliHeader* header = 0; AliGenEventHeader* genHeader = 0; AliStack* stack = 0; 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->GetXv(); vert[1] = vtxESD->GetYv(); vert[2] = vtxESD->GetZv(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); Double_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); // loop over MC stack for (Int_t iMc = 0; iMc < 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.="<<&fileName<< "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(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } // get file name TTree *chain = (TChain*)GetInputData(0); if(!chain) { Printf("ERROR: Could not receive input chain"); return; } TObjString fileName(chain->GetCurrentFile()->GetName()); // trigger if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; //SetPhysicsTriggerSelection(physicsSelection); if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // centrality determination Float_t centralityF = -1; AliCentrality *esdCentrality = esdEvent->GetCentrality(); centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data()); // get reconstructed vertex //const AliESDVertex* vtxESD = 0; AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } if(!vtxESD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); // check event cuts if(isEventOK && isEventTriggered) { // // Dump the pt downscaled V0 into the tree // Int_t ntracks = esdEvent->GetNumberOfTracks(); Int_t nV0s = esdEvent->GetNumberOfV0s(); Int_t run = esdEvent->GetRunNumber(); Int_t time= esdEvent->GetTimeStamp(); Int_t evNr=esdEvent->GetEventNumberInFile(); Double_t bz = esdEvent->GetMagneticField(); // 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()) { friendTrack0 = esdFriend->GetTrack(v0->GetIndex(0)); //this guy can be NULL friendTrack1 = esdFriend->GetTrack(v0->GetIndex(1)); //this guy can be NULL } } //Bool_t newFriendTrack0=kFALSE; //Bool_t newFriendTrack1=kFALSE; //if (!friendTrack0) {friendTrack0=new AliESDfriendTrack(); newFriendTrack0=kTRUE;} //if (!friendTrack1) {friendTrack1=new AliESDfriendTrack(); newFriendTrack1=kTRUE;} if (!friendTrack0) {friendTrack0=fDummyFriendTrack;} if (!friendTrack1) {friendTrack1=fDummyFriendTrack;} if (track0->GetSign()<0) { track1 = esdEvent->GetTrack(v0->GetIndex(0)); track0 = esdEvent->GetTrack(v0->GetIndex(1)); } // Bool_t isDownscaled = IsV0Downscaled(v0); if (isDownscaled) continue; AliKFParticle kfparticle; // Int_t type=GetKFParticle(v0,esdEvent,kfparticle); if (type==0) continue; TObjString triggerClass = esdEvent->GetFiredTriggerClasses().Data(); if(!fFillTree) return; if(!fTreeSRedirector) return; (*fTreeSRedirector)<<"V0s"<< "gid="<GetCurrentFile()->GetName()); // trigger Bool_t isEventTriggered = kTRUE; AliPhysicsSelection *physicsSelection = NULL; AliTriggerAnalysis* triggerAnalysis = NULL; // AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); if (!inputHandler) { Printf("ERROR: Could not receive input handler"); return; } if(evtCuts->IsTriggerRequired()) { // always MB isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB; physicsSelection = static_cast (inputHandler->GetEventSelection()); if(!physicsSelection) return; if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) { // set trigger (V0AND) triggerAnalysis = physicsSelection->GetTriggerAnalysis(); if(!triggerAnalysis) return; isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger()); } } // get reconstructed vertex AliESDVertex* vtxESD = 0; if(GetAnalysisMode() == kTPCAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTPC(); } else if(GetAnalysisMode() == kTPCITSAnalysisMode) { vtxESD = (AliESDVertex*)esdEvent->GetPrimaryVertexTracks(); } else { return; } if(!vtxESD) return; Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered); //printf("GetAnalysisMode() %d \n",GetAnalysisMode()); // check event cuts if(isEventOK && isEventTriggered) { Double_t vert[3] = {0}; vert[0] = vtxESD->GetXv(); vert[1] = vtxESD->GetYv(); vert[2] = vtxESD->GetZv(); Int_t mult = vtxESD->GetNContributors(); Double_t bz = esdEvent->GetMagneticField(); Double_t runNumber = esdEvent->GetRunNumber(); Double_t evtTimeStamp = esdEvent->GetTimeStamp(); Int_t evtNumberInFile = esdEvent->GetEventNumberInFile(); // 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) {if (!esdFriend->TestSkipBit()) friendTrack = esdFriend->GetTrack(iTrack);} //this guy can be NULL //Bool_t newFriendTrack=kFALSE; //if (!friendTrack) {friendTrack=new AliESDfriendTrack(); newFriendTrack=kTRUE;} if (!friendTrack) {friendTrack=fDummyFriendTrack;} 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"<< "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; //if ((track0->GetITSclusters(0))<2) return 0; //if ((track1->GetITSclusters(0))<2) return 0; Float_t pos0[2]={0}, cov0[3]={0}; Float_t pos1[2]={0}, cov1[3]={0}; track0->GetImpactParameters(pos0,cov0); track0->GetImpactParameters(pos1,cov1); // if (TMath::Abs(pos0[0])GetKFInfo(2,2,2); if (chi2KF>25) return 0; // // 4.) Rough mass cut - 0.200 GeV // static Double_t masses[2]={-1}; if (masses[0]<0){ masses[0] = TDatabasePDG::Instance()->GetParticle("K_S0")->Mass(); masses[1] = TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass(); } Double_t mass00= v0->GetEffMass(0,0); Double_t mass22= v0->GetEffMass(2,2); Double_t mass42= v0->GetEffMass(4,2); Double_t mass24= v0->GetEffMass(2,4); Bool_t massOK=kFALSE; Int_t type=0; Int_t ptype=0; Double_t dmass=1; Int_t p1=0, p2=0; if (TMath::Abs(mass00-0)GetParamP(); const AliExternalTrackParam *paramN = v0->GetParamN(); if (paramP->GetSign()<0){ paramP=v0->GetParamP(); paramN=v0->GetParamN(); } //Double_t *pparam1 = (Double_t*)paramP->GetParameter(); //Double_t *pparam2 = (Double_t*)paramN->GetParameter(); // AliKFParticle kfp1( *paramP, spdg[p1] ); AliKFParticle kfp2( *paramN, -1 *spdg[p2] ); AliKFParticle V0KF; (V0KF)+=kfp1; (V0KF)+=kfp2; kfparticle=V0KF; // // Pointing angle // Double_t errPhi = V0KF.GetErrPhi(); Double_t pointAngle= TMath::ACos(v0->GetV0CosineOfPointingAngle()); if (pointAngle/errPhi>10) return 0; // return ptype; } //_____________________________________________________________________________ Bool_t AliAnalysisTaskFilteredTree::IsV0Downscaled(AliESDv0 *const v0) { // // Downscale randomly low pt V0 // //return kFALSE; Double_t maxPt= TMath::Max(v0->GetParamP()->Pt(), v0->GetParamN()->Pt()); Double_t scalempt= TMath::Min(maxPt,10.); Double_t downscaleF = gRandom->Rndm(); downscaleF *= fLowPtV0DownscaligF; // // 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 // //// must be deleted to store trees //if(fTreeSRedirector) delete fTreeSRedirector; fTreeSRedirector=0; //// open temporary file and copy trees to the ouptut container //TChain* chain = 0; //// //chain = new TChain("highPt"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fHighPtTree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fHighPtTree) fHighPtTree->Print(); //// //chain = new TChain("V0s"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fV0Tree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fV0Tree) fV0Tree->Print(); //// //chain = new TChain("dEdx"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fdEdxTree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fdEdxTree) fdEdxTree->Print(); //// //chain = new TChain("Laser"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fLaserTree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fLaserTree) fLaserTree->Print(); //// //chain = new TChain("MCEffTree"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fMCEffTree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fMCEffTree) fMCEffTree->Print(); //// //chain = new TChain("CosmicPairs"); //if(chain) { // chain->Add("jotwinow_Temp_Trees.root"); // fCosmicPairsTree = chain->CopyTree("1"); // delete chain; chain=0; //} //if(fCosmicPairsTree) fCosmicPairsTree->Print(); Bool_t deleteTrees=kTRUE; if ((AliAnalysisManager::GetAnalysisManager())) { if (AliAnalysisManager::GetAnalysisManager()->GetAnalysisType() == AliAnalysisManager::kProofAnalysis) deleteTrees=kFALSE; } if (deleteTrees) delete fTreeSRedirector; fTreeSRedirector=NULL; //OpenFile(1); // Post output data. //PostData(1, fHighPtTree); //PostData(2, fV0Tree); //PostData(3, fdEdxTree); //PostData(4, fLaserTree); //PostData(5, fMCEffTree); //PostData(6, fCosmicPairsTree); } //_____________________________________________________________________________ void AliAnalysisTaskFilteredTree::Terminate(Option_t *) { // Called one at the end /* fOutputSummary = dynamic_cast (GetOutputData(1)); if(fOutputSummary) delete fOutputSummary; fOutputSummary=0; TChain* chain = new TChain("highPt"); if(!chain) return; chain->Add("jotwinow_HighPt_TrackAndV0_Trees.root"); TTree *tree = chain->CopyTree("1"); if (chain) { delete chain; chain=0; } if(!tree) return; tree->Print(); fOutputSummary = tree; if (!fOutputSummary) { Printf("ERROR: AliAnalysisTaskFilteredTree::Terminate(): Output data not avaiable %p \n", GetOutputData(1)); return; } */ } //_____________________________________________________________________________ Int_t AliAnalysisTaskFilteredTree::GetMCTrueTrackMult(AliMCEvent *const mcEvent, AliFilteredTreeEventCuts *const evtCuts, AliFilteredTreeAcceptanceCuts *const accCuts) { // // calculate mc event true track multiplicity // if(!mcEvent) return 0; AliStack* stack = 0; Int_t mult = 0; // MC particle stack stack = mcEvent->Stack(); if (!stack) return 0; // //printf("minZv %f, maxZv %f \n", evtCuts->GetMinZv(), evtCuts->GetMaxZv()); // Bool_t isEventOK = evtCuts->AcceptMCEvent(mcEvent); if(!isEventOK) return 0; Int_t nPart = stack->GetNtrack(); for (Int_t iMc = 0; iMc < nPart; ++iMc) { TParticle* particle = stack->Particle(iMc); if (!particle) continue; // only charged particles if(!particle->GetPDG()) continue; Double_t charge = particle->GetPDG()->Charge()/3.; if (TMath::Abs(charge) < 0.001) continue; // physical primary Bool_t prim = stack->IsPhysicalPrimary(iMc); if(!prim) continue; // checked accepted without pt cut //if(accCuts->AcceptTrack(particle)) if( particle->Eta() > accCuts->GetMinEta() && particle->Eta() < accCuts->GetMaxEta() ) { mult++; } } return mult; } //_____________________________________________________________________________ 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())); } }