/************************************************************************** * Copyright(c) 1998-2009, 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 AliAnalysisTaskVertexESD // AliAnalysisTask to extract from ESD the information for the analysis // of the primary vertex reconstruction efficiency and resolution // (for MC events) and distributions (for real data). Three vertices: // - SPD tracklets // - ITS+TPC tracks // - TPC-only tracks // // Author: A.Dainese, andrea.dainese@pd.infn.it //************************************************************************* #include #include #include #include #include #include #include #include #include #include #include #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliMultiplicity.h" #include "AliVertexerTracks.h" #include "AliESDtrack.h" #include "AliExternalTrackParam.h" #include "AliESDVertex.h" #include "AliVEvent.h" #include "AliESDInputHandler.h" #include "AliTrackReference.h" //#include "AliTriggerAnalysis.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliStack.h" #include "AliLog.h" #include "AliGenEventHeader.h" #include "AliAnalysisTaskVertexESD.h" ClassImp(AliAnalysisTaskVertexESD) //________________________________________________________________________ AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : AliAnalysisTaskSE(name), fCheckEventType(kTRUE), fReadMC(kFALSE), fRecoVtxTPC(kTRUE), fRecoVtxITSTPC(kTRUE), fRecoVtxITSTPCHalfEvent(kFALSE), fOnlyITSTPCTracks(kFALSE), fOnlyITSSATracks(kFALSE), fFillNtuple(kFALSE), fFillTreeBeamSpot(kFALSE), fESD(0), fOutput(0), fNtupleVertexESD(0), fhSPDVertexX(0), fhSPDVertexY(0), fhSPDVertexZ(0), fhSPDVertexZonly(0), fhTRKVertexX(0), fhTRKVertexY(0), fhTRKVertexZ(0), fhTPCVertexX(0), fhTPCVertexY(0), fhTPCVertexZ(0), fhTrackRefs(0), fTreeBeamSpot(0), fhTriggeredTrklets(0), fhSPD3DTrklets(0), fhSPDZTrklets(0), fhTRKTrklets(0), fhTRKcTrklets(0), fhTRKncTrklets(0), fhTPCTrklets(0), fhTPCcTrklets(0), fhTPCncTrklets(0), fhSPD3DZreco(0), fhSPDZZreco(0), fhSPDVertexXPile(0), fhSPDVertexYPile(0), fhSPDVertexZPile(0), fhSPDVertexDiffZPileContr2(0), fhSPDVertexDiffZPileContr3(0), fhSPDVertexDiffZPileContr4(0), fhSPDVertexDiffZPileContr5(0), fhSPDVertexDiffZPileDefault(0), fhSPDContributorsPile(0), fhSPDDispContributors(0), fTriggerType(AliVEvent::kMB), fhntrksSPDvsSPDcls(0), fhntrksZvsSPDcls(0) { // Constructor // Define input and output slots here // Output slot #0 writes into a TList container DefineOutput(1, TList::Class()); //My private output } //________________________________________________________________________ AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD() { // Destructor // histograms are in the output list and deleted when the output // list is deleted by the TSelector dtor if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fOutput; fOutput = 0; } } //________________________________________________________________________ void AliAnalysisTaskVertexESD::UserCreateOutputObjects() { // Create histograms // Called once // Several histograms are more conveniently managed in a TList fOutput = new TList; fOutput->SetOwner(); fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","run:tstamp:bunchcross:triggered:dndygen:xdiam:ydiam:zdiam:xerrdiam:yerrdiam:zerrdiam:xtrue:ytrue:ztrue:xSPD:xerrSPD:ySPD:yerrSPD:zSPD:zerrSPD:ntrksSPD:SPD3D:dphiSPD:xTPC:xerrTPC:yTPC:yerrTPC:zTPC:zerrTPC:ntrksTPC:constrTPC:xTRK:xerrTRK:yTRK:yerrTRK:zTRK:zerrTRK:ntrksTRK:constrTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:SPD0cls:ispileupfromSPDdefault:xSPDp:xerrSPDp:ySPDp:yerrSPDp:zSPDp:zerrSPDp:ntrksSPDp:xTPCnc:xerrTPCnc:yTPCnc:yerrTPCnc:zTPCnc:zerrTPCnc:ntrksTPCnc:xTPCc:xerrTPCc:yTPCc:yerrTPCc:zTPCc:zerrTPCc:ntrksTPCc:xTRKnc:xerrTRKnc:yTRKnc:yerrTRKnc:zTRKnc:zerrTRKnc:ntrksTRKnc:xTRKc:xerrTRKc:yTRKc:yerrTRKc:zTRKc:zerrTRKc:ntrksTRKc:deltaxTRKnc:deltayTRKnc:deltazTRKnc:deltaxerrTRKnc:deltayerrTRKnc:deltazerrTRKnc:ntrksEvenTRKnc:ntrksOddTRKnc"); fOutput->Add(fNtupleVertexESD); fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1); fOutput->Add(fhSPDVertexX); fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1); fOutput->Add(fhSPDVertexY); fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",250,-25,25); fOutput->Add(fhSPDVertexZ); fhSPDVertexZonly = new TH1F("fhSPDVertexZonly","SPDVertexer z; z vertex [cm]; events",250,-25,25); fOutput->Add(fhSPDVertexZonly); fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",2000,-1,1); fOutput->Add(fhTRKVertexX); fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",2000,-1,1); fOutput->Add(fhTRKVertexY); fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",250,-25,25); fOutput->Add(fhTRKVertexZ); fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3); fOutput->Add(fhTPCVertexX); fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3); fOutput->Add(fhTPCVertexY); fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",250,-25,25); fOutput->Add(fhTPCVertexZ); fhSPDVertexXPile = new TH1F("fhSPDVertexXPile","SPDVertexPile x; x vertex [cm]; events",200,-20,20); fOutput->Add(fhSPDVertexXPile); fhSPDVertexYPile = new TH1F("fhSPDVertexYPile","SPDVertexPile y; y vertex [cm]; events",200,-20,20); fOutput->Add(fhSPDVertexYPile); fhSPDVertexZPile = new TH1F("fhSPDVertexZPile","SPDVertexPile z; z vertex [cm]; events",200,-40,40); fOutput->Add(fhSPDVertexZPile); fhSPDVertexDiffZPileContr2 = new TH1F("fhSPDVertexDiffZPileContr2","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100); fOutput->Add(fhSPDVertexDiffZPileContr2); fhSPDVertexDiffZPileContr3 = new TH1F("fhSPDVertexDiffZPileContr3","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100); fOutput->Add(fhSPDVertexDiffZPileContr3); fhSPDVertexDiffZPileContr4 = new TH1F("fhSPDVertexDiffZPileContr4","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100); fOutput->Add(fhSPDVertexDiffZPileContr4); fhSPDVertexDiffZPileContr5 = new TH1F("fhSPDVertexDiffZPileContr5","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100); fOutput->Add(fhSPDVertexDiffZPileContr5); fhSPDVertexDiffZPileDefault = new TH1F("fhSPDVertexDiffZPileDefault","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100); fOutput->Add(fhSPDVertexDiffZPileDefault); fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4); fOutput->Add(fhTrackRefs); fTreeBeamSpot = new TTree("fTreeBeamSpot", "BeamSpotTree"); UShort_t triggered, ntrkletsS, ntrksTRKnc; UInt_t run, bx; Float_t cetTimeLHC,xTRKnc, yTRKnc, zTRKnc; fTreeBeamSpot->Branch("run", &run, "run/i"); fTreeBeamSpot->Branch("cetTimeLHC", &cetTimeLHC, "cetTimeLHC/F"); fTreeBeamSpot->Branch("bx", &bx, "bx/i"); fTreeBeamSpot->Branch("triggered", &triggered, "triggered/s"); fTreeBeamSpot->Branch("ntrkletsS", &ntrkletsS, "ntrkletsS/s"); fTreeBeamSpot->Branch("xTRKnc", &xTRKnc, "xTRKnc/F"); fTreeBeamSpot->Branch("yTRKnc", &yTRKnc, "yTRKnc/F"); fTreeBeamSpot->Branch("zTRKnc", &zTRKnc, "zTRKnc/F"); fTreeBeamSpot->Branch("ntrksTRKnc", &ntrksTRKnc, "ntrksTRKnc/s"); fOutput->Add(fTreeBeamSpot); Int_t nbinTrklets=29; Float_t lowTrklets[30]={-0.5,0.5,1.5,2.5,3.5,4.5,5.5,7.5,10.5,25.5,50.5,75.5,100.5,150.5,200.5,300.5,400.5,500.5,600.5,800.5,1000.5,1500.5,2000.5,2500.5,3000.5,4000.5,5000.5,6000.5,8000.5,10000.5}; fhTriggeredTrklets = new TH1F("fhTriggeredTrklets","trklets dist for triggered ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTriggeredTrklets); fhSPD3DTrklets = new TH1F("fhSPD3DTrklets","trklets dist for SPD3D ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhSPD3DTrklets); fhSPDZTrklets = new TH1F("fhSPDZTrklets","trklets dist for SPDZ ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhSPDZTrklets); fhTRKTrklets = new TH1F("fhTRKTrklets","trklets dist for TRK ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTRKTrklets); fhTRKcTrklets = new TH1F("fhTRKcTrklets","trklets dist for TRKc ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTRKcTrklets); fhTRKncTrklets = new TH1F("fhTRKncTrklets","trklets dist for TRKnc ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTRKncTrklets); fhTPCTrklets = new TH1F("fhTPCTrklets","trklets dist for TPC ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTPCTrklets); fhTPCcTrklets = new TH1F("fhTPCcTrklets","trklets dist for TPCc ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTPCcTrklets); fhTPCncTrklets = new TH1F("fhTPCncTrklets","trklets dist for TPCnc ev.; ntrklets; entries",nbinTrklets,lowTrklets); fOutput->Add(fhTPCncTrklets); Int_t nbinZreco = 16; Double_t lowZreco[17]={-15.0,-10.0,-7.0,-5,-4,-3,-2,-1,0,1,2,3,4,5,7,10,15}; fhSPD3DZreco = new TH1F("fhSPD3DZreco","Zreco dist for SPD3D ev.; Zreco [cm]; entries",nbinZreco,lowZreco); fOutput->Add(fhSPD3DZreco); fhSPDZZreco = new TH1F("fhSPDZZreco","Zreco dist for SPDZ ev.; Zreco [cm]; entries",nbinZreco,lowZreco); fOutput->Add(fhSPDZZreco); fhSPDContributorsPile = new TH1F("fhSPDContributorsPile","ncontributors pile up vertex; ncontributors; entries",200,-0.5,199.5); fOutput->Add(fhSPDContributorsPile); fhSPDDispContributors = new TH2F("fhSPDDispContributors","ncontributors main-pile; ncontributors main; ncontributors pile",200,-0.5,199.5,200,-0.5,199.5); fOutput->Add(fhSPDDispContributors); fhntrksSPDvsSPDcls = new TH2F("fhntrksSPDvsSPDcls", "ncontributors SPD3D vs number of cluster SPD; nContributors 3D; nCluster SPD1",300,0.,3000.,400,0.,8000.); fOutput->Add(fhntrksSPDvsSPDcls); fhntrksZvsSPDcls = new TH2F("fhntrksZvsSPDcls", "ncontributors SPDZ vs number of cluster SPD; nContributors Z; nCluster SPD1",300,0.,3000.,400,0.,8000.); fOutput->Add(fhntrksZvsSPDcls); PostData(1, fOutput); return; } //________________________________________________________________________ void AliAnalysisTaskVertexESD::UserExec(Option_t *) { // Main loop // Called for each event if (!InputEvent()) { Printf("ERROR: fESD not available"); return; } AliESDEvent* esdE = (AliESDEvent*) InputEvent(); // Select PHYSICS events (type=7, for data) and MC events (type=0) // fCheckEventType is kFALSE if fReadMC is kTRUE, hence check is skipped if(fCheckEventType) { if(esdE->GetEventType()!=7 && esdE->GetEventType()!=0) return; } TArrayF mcVertex(3); mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.; Float_t dNchdy=-999.; // *********** MC info *************** if (fReadMC) { AliMCEventHandler* eventHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); if (!eventHandler) { Printf("ERROR: Could not retrieve MC event handler"); return; } AliMCEvent* mcEvent = eventHandler->MCEvent(); if (!mcEvent) { Printf("ERROR: Could not retrieve MC event"); return; } AliStack* stack = mcEvent->Stack(); if (!stack) { AliDebug(AliLog::kError, "Stack not available"); return; } AliHeader* header = mcEvent->Header(); if (!header) { AliDebug(AliLog::kError, "Header not available"); return; } AliGenEventHeader* genHeader = header->GenEventHeader(); genHeader->PrimaryVertex(mcVertex); Int_t ngenpart = (Int_t)stack->GetNtrack(); //printf("# generated particles = %d\n",ngenpart); dNchdy=0; for(Int_t ip=0; ipParticle(ip); // keep only electorns, muons, pions, kaons and protons Int_t apdg = TMath::Abs(part->GetPdgCode()); if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue; // reject secondaries if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue; // reject incoming protons Double_t energy = part->Energy(); if(energy>900.) continue; Double_t pz = part->Pz(); Double_t y = 100.; if(energy-pz+1.e-13>0.) y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13)); if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1 // tracks refs TClonesArray *trefs=0; Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs); if(ntrefs<0) continue; for(Int_t iref=0; irefAt(iref); if(tref->R()>10.) continue; fhTrackRefs->Fill(tref->X(),tref->Y()); } } //printf("# primary particles = %7.1f\n",dNchdy); } // *********** MC info *************** // Trigger //ULong64_t triggerMask; //ULong64_t spdFO = (1 << 14); //ULong64_t v0left = (1 << 10); //ULong64_t v0right = (1 << 11); //triggerMask=esdE->GetTriggerMask(); // MB1: SPDFO || V0L || V0R //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); //MB2: GFO && V0R //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right)) //Bool_t eventTriggered = (triggerMask & spdFO); //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis(); Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/); // use response of AliPhysicsSelection eventTriggered = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerType); Int_t nInFile = esdE->GetEventNumberInFile(); const AliMultiplicity *alimult = esdE->GetMultiplicity(); Int_t ntrklets=0,spd0cls=0, spd1cls=0; if(alimult) { ntrklets = alimult->GetNumberOfTracklets(); for(Int_t l=0;lGetNumberOfTracklets();l++){ if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--; } spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets; spd1cls = alimult->GetNumberOfITSClusters(1); } UShort_t triggered, ntrkletsS, ntrksTRKnc; UInt_t run, bx; Float_t cetTimeLHC,xTRKnc, yTRKnc, zTRKnc; fTreeBeamSpot->SetBranchAddress("run", &run); fTreeBeamSpot->SetBranchAddress("cetTimeLHC", &cetTimeLHC); fTreeBeamSpot->SetBranchAddress("bx", &bx); fTreeBeamSpot->SetBranchAddress("triggered", &triggered); fTreeBeamSpot->SetBranchAddress("ntrkletsS", &ntrkletsS); fTreeBeamSpot->SetBranchAddress("xTRKnc", &xTRKnc); fTreeBeamSpot->SetBranchAddress("yTRKnc", &yTRKnc); fTreeBeamSpot->SetBranchAddress("zTRKnc", &zTRKnc); fTreeBeamSpot->SetBranchAddress("ntrksTRKnc", &ntrksTRKnc); Double_t tstamp = esdE->GetTimeStamp(); Float_t cetTime =(tstamp-1262304000.)+7200.; Float_t cetTime1h =(tstamp-1262304000.)+3600.; cetTimeLHC = (Float_t)cetTime1h; Int_t ntracks = esdE->GetNumberOfTracks(); Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0; //printf("Tracks # = %d\n",esdE->GetNumberOfTracks()); for(Int_t itr=0; itrGetTrack(itr); if(t->GetNcls(0)>=5) nITS5or6++; Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0); if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) { nTPCin++; if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++; } } const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD(); const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0); const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC(); const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks(); // fill histos if(spdv) { if(spdv->GetNContributors()>0) { TString title=spdv->GetTitle(); if(title.Contains("3D")) { fhSPDVertexX->Fill(spdv->GetX()); fhSPDVertexY->Fill(spdv->GetY()); fhSPDVertexZ->Fill(spdv->GetZ()); fhntrksSPDvsSPDcls->Fill(spdv->GetNContributors(), spd0cls); } if(title.Contains("Z")){ fhSPDVertexZonly->Fill(spdv->GetZ()); fhntrksZvsSPDcls->Fill(spdv->GetNContributors(), spd1cls); } } } if(trkv) { if(trkv->GetNContributors()>0) { fhTRKVertexX->Fill(trkv->GetX()); fhTRKVertexY->Fill(trkv->GetY()); fhTRKVertexZ->Fill(trkv->GetZ()); } } if(tpcv) { if(tpcv->GetNContributors()>0) { fhTPCVertexX->Fill(tpcv->GetX()); fhTPCVertexY->Fill(tpcv->GetY()); fhTPCVertexZ->Fill(tpcv->GetZ()); } } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // filling Vertex reco efficiency plots if(eventTriggered ? 1. : 0.){ fhTriggeredTrklets->Fill(ntrklets); if(spdv){ if(spdv->GetNContributors()>0.5){ TString spdtitle = spdv->GetTitle(); if(spdtitle.Contains("vertexer: 3D") ? 1. : 0.){ fhSPD3DTrklets->Fill(ntrklets); fhSPD3DZreco->Fill(spdv->GetZ()); } else{ fhSPDZTrklets->Fill(ntrklets); fhSPDZZreco->Fill(spdv->GetZ()); } } } if(trkv){ if(trkv->GetNContributors()>0.5)fhTRKTrklets->Fill(ntrklets); } if(fRecoVtxITSTPC) { AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE); if(trkvc->GetNContributors()>0.5)fhTRKcTrklets->Fill(ntrklets); delete trkvc; trkvc=0; AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE); if(trkvnc->GetNContributors()>0.5)fhTRKncTrklets->Fill(ntrklets); delete trkvnc; trkvnc=0; } if(tpcv){ if(tpcv->GetNContributors()>0.5)fhTPCTrklets->Fill(ntrklets); } if(fRecoVtxTPC) { AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE); if(tpcvc->GetNContributors()>0.5)fhTPCcTrklets->Fill(ntrklets); delete tpcvc; tpcvc=0; AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE); if(tpcvnc->GetNContributors()>0.5)fhTPCncTrklets->Fill(ntrklets); delete tpcvnc; tpcvnc=0; } } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// Float_t xpile=-999.; Float_t ypile=-999.; Float_t zpile=-999.; Float_t expile=-999.; Float_t eypile=-999.; Float_t ezpile=-999.; Int_t ntrkspile=-1; Bool_t isPileUpfromSPD=esdE->IsPileupFromSPD(); if(esdE->GetNumberOfPileupVerticesSPD()>0 && spdvp && spdv){ if(spdvp->GetNContributors()>=2) { xpile=spdvp->GetX(); expile=spdvp->GetXRes(); ypile=spdvp->GetY(); eypile=spdvp->GetYRes(); zpile=spdvp->GetZ(); ezpile=spdvp->GetZRes(); ntrkspile=spdvp->GetNContributors(); if (esdE->IsPileupFromSPD(2,0.,3.,2.,5.)){ fhSPDVertexDiffZPileContr2->Fill(spdv->GetZ()-spdvp->GetZ()); if(spdvp->GetNContributors()>=3) {fhSPDVertexDiffZPileContr3->Fill(spdv->GetZ()-spdvp->GetZ());} if(spdvp->GetNContributors()>=4) {fhSPDVertexDiffZPileContr4->Fill(spdv->GetZ()-spdvp->GetZ());} if(spdvp->GetNContributors()>=5) {fhSPDVertexDiffZPileContr5->Fill(spdv->GetZ()-spdvp->GetZ());} }//end IsPileUpFromSPD if (isPileUpfromSPD){ fhSPDVertexXPile->Fill(spdvp->GetX()); fhSPDVertexYPile->Fill(spdvp->GetY()); fhSPDVertexZPile->Fill(spdvp->GetZ()); fhSPDContributorsPile->Fill(spdvp->GetNContributors()); fhSPDDispContributors->Fill(spdv->GetNContributors(),spdvp->GetNContributors()); fhSPDVertexDiffZPileDefault->Fill(spdv->GetZ()-spdvp->GetZ()); } } } // fill ntuple Int_t isize=89; Float_t xnt[89]; for(Int_t iii=0;iiiGetRunNumber(); secnt[indexSecNt++]=(Float_t)esdE->GetRunNumber(); run = (Int_t)esdE->GetRunNumber(); xnt[index++]=cetTime; //(Float_t)esdE->GetTimeStamp(); //secnt[indexSecNt++]=cetTime; secnt[indexSecNt++]=cetTime1h; xnt[index++]=(Float_t)esdE->GetBunchCrossNumber(); secnt[indexSecNt++]=(Float_t)esdE->GetBunchCrossNumber(); bx = (Int_t)esdE->GetBunchCrossNumber(); xnt[index++]=(eventTriggered ? 1. : 0.); secnt[indexSecNt++]=(eventTriggered ? 1. : 0.); triggered = (UShort_t)(eventTriggered ? 1 : 0); xnt[index++]=(Float_t)dNchdy; xnt[index++]=(Float_t)esdE->GetDiamondX(); xnt[index++]=(Float_t)esdE->GetDiamondY(); xnt[index++]=(Float_t)esdE->GetDiamondZ(); xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondX())); xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondY())); xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondZ())); xnt[index++]=mcVertex[0]; xnt[index++]=mcVertex[1]; xnt[index++]=mcVertex[2]; xnt[index++]=spdv->GetX(); xnt[index++]=spdv->GetXRes(); xnt[index++]=spdv->GetY(); xnt[index++]=spdv->GetYRes(); xnt[index++]=spdv->GetZ(); xnt[index++]=spdv->GetZRes(); xnt[index++]=spdv->GetNContributors(); TString spdtitle = spdv->GetTitle(); xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.); xnt[index++]=spdv->GetDispersion(); xnt[index++]=tpcv->GetX(); xnt[index++]=tpcv->GetXRes(); xnt[index++]=tpcv->GetY(); xnt[index++]=tpcv->GetYRes(); xnt[index++]=tpcv->GetZ(); xnt[index++]=tpcv->GetZRes(); xnt[index++]=tpcv->GetNContributors(); TString tpctitle = tpcv->GetTitle(); xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.); xnt[index++]=trkv->GetX(); xnt[index++]=trkv->GetXRes(); xnt[index++]=trkv->GetY(); xnt[index++]=trkv->GetYRes(); xnt[index++]=trkv->GetZ(); xnt[index++]=trkv->GetZRes(); xnt[index++]=trkv->GetNContributors();// tpccontrorig; TString trktitle = trkv->GetTitle(); xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.); xnt[index++]=float(ntrklets); secnt[indexSecNt++]=float(ntrklets); ntrkletsS = (UShort_t)ntrklets; xnt[index++]=float(ntracks); xnt[index++]=float(nITS5or6); xnt[index++]=float(nTPCin); xnt[index++]=float(nTPCinEta09); xnt[index++]=spd0cls; xnt[index++]=Float_t(isPileUpfromSPD); xnt[index++]=xpile; xnt[index++]=expile; xnt[index++]=ypile; xnt[index++]=eypile; xnt[index++]=zpile; xnt[index++]=ezpile; xnt[index++]=(Float_t)ntrkspile; // add recalculated vertices TRK and TPC if(fRecoVtxTPC) { AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE); xnt[index++]=tpcvnc->GetX(); xnt[index++]=tpcvnc->GetXRes(); xnt[index++]=tpcvnc->GetY(); xnt[index++]=tpcvnc->GetYRes(); xnt[index++]=tpcvnc->GetZ(); xnt[index++]=tpcvnc->GetZRes(); xnt[index++]=tpcvnc->GetNContributors(); delete tpcvnc; tpcvnc=0; AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE); xnt[index++]=tpcvc->GetX(); xnt[index++]=tpcvc->GetXRes(); xnt[index++]=tpcvc->GetY(); xnt[index++]=tpcvc->GetYRes(); xnt[index++]=tpcvc->GetZ(); xnt[index++]=tpcvc->GetZRes(); xnt[index++]=tpcvc->GetNContributors(); delete tpcvc; tpcvc=0; } else index+=14; if(fRecoVtxITSTPC) { AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE); xnt[index++]=trkvnc->GetX(); xnt[index++]=trkvnc->GetXRes(); xnt[index++]=trkvnc->GetY(); xnt[index++]=trkvnc->GetYRes(); xnt[index++]=trkvnc->GetZ(); xnt[index++]=trkvnc->GetZRes(); xnt[index++]=trkvnc->GetNContributors(); secnt[indexSecNt++]=trkvnc->GetX(); secnt[indexSecNt++]=trkvnc->GetY(); secnt[indexSecNt++]=trkvnc->GetZ(); secnt[indexSecNt++]=trkvnc->GetNContributors(); xTRKnc = (Float_t)trkvnc->GetX(); yTRKnc = (Float_t)trkvnc->GetY(); zTRKnc = (Float_t)trkvnc->GetZ(); ntrksTRKnc = (trkvnc->GetNContributors()<0 ? 0 : (UShort_t)trkvnc->GetNContributors()); delete trkvnc; trkvnc=0; AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE); xnt[index++]=trkvc->GetX(); xnt[index++]=trkvc->GetXRes(); xnt[index++]=trkvc->GetY(); xnt[index++]=trkvc->GetYRes(); xnt[index++]=trkvc->GetZ(); xnt[index++]=trkvc->GetZRes(); xnt[index++]=trkvc->GetNContributors(); delete trkvc; trkvc=0; } else index+=14; if(fRecoVtxITSTPCHalfEvent) { AliESDVertex *trkvncodd = ReconstructPrimaryVertexITSTPC(kFALSE,1); AliESDVertex *trkvnceven = ReconstructPrimaryVertexITSTPC(kFALSE,2); if(trkvncodd->GetNContributors()>0 && trkvnceven->GetNContributors()>0) { xnt[index++]=trkvncodd->GetX()-trkvnceven->GetX(); xnt[index++]=trkvncodd->GetY()-trkvnceven->GetY(); xnt[index++]=trkvncodd->GetZ()-trkvnceven->GetZ(); xnt[index++]=TMath::Sqrt(trkvncodd->GetXRes()*trkvncodd->GetXRes()+trkvnceven->GetXRes()*trkvnceven->GetXRes()); xnt[index++]=TMath::Sqrt(trkvncodd->GetYRes()*trkvncodd->GetYRes()+trkvnceven->GetYRes()*trkvnceven->GetYRes()); xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()*trkvncodd->GetZRes()+trkvnceven->GetZRes()*trkvnceven->GetZRes()); xnt[index++]=trkvnceven->GetNContributors(); xnt[index++]=trkvncodd->GetNContributors(); } else { xnt[index++]=0.; xnt[index++]=0.; xnt[index++]=0.; xnt[index++]=0.; xnt[index++]=0.; xnt[index++]=0.; xnt[index++]=-1; xnt[index++]=-1; } delete trkvncodd; trkvncodd=0; delete trkvnceven; trkvnceven=0; } else index+=8; if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n"); if(fFillNtuple) fNtupleVertexESD->Fill(xnt); if(indexSecNt>isizeSecNt) printf("AliAnalysisTaskVertexESD: ERROR, indexSecNt!=isizeSecNt\n"); // if(indexTree>isizeTree) printf("AliAnalysisTaskVertexESD: ERROR, indexTree!=isizeTree\n"); // only every second event (to reduce output size) if(fFillTreeBeamSpot && (nInFile%2 == 0)) fTreeBeamSpot->Fill(); // Post the data already here PostData(1, fOutput); return; } //________________________________________________________________________ void AliAnalysisTaskVertexESD::Terminate(Option_t *) { // Draw result to the screen // Called once at the end of the query fOutput = dynamic_cast (GetOutputData(1)); if (!fOutput) { Printf("ERROR: fOutput not available"); return; } ////////////////////////////////////////////////////// /* TH1F *fhTriggeredTrklets=(TH1F*)fOutput->FindObject("fhTriggeredTrklets"); TH1F *fhSPDZTrklets=(TH1F*)fOutput->FindObject("fhSPDZTrklets"); TH1F *fhSPD3DTrklets=(TH1F*)fOutput->FindObject("fhSPD3DTrklets"); TH1F *fhTRKTrklets=(TH1F*)fOutput->FindObject("fhTRKTrklets"); TH1F *fhTRKcTrklets=(TH1F*)fOutput->FindObject("fhTRKcTrklets"); TH1F *fhTRKncTrklets=(TH1F*)fOutput->FindObject("fhTRKncTrklets"); TH1F *fhSPDZZreco=(TH1F*)fOutput->FindObject("fhSPDZZreco"); TH1F *fhSPD3DZreco=(TH1F*)fOutput->FindObject("fhSPD3DZreco"); TGraphAsymmErrors *fhSPDZEffTrklets=new TGraphAsymmErrors(fhSPDZTrklets,fhTriggeredTrklets,"w"); fhSPDZEffTrklets->SetName("fhSPDZEffTrklets"); fhSPDZEffTrklets->SetDrawOption("AP"); TGraphAsymmErrors *fhSPD3DEffTrklets=new TGraphAsymmErrors(fhSPD3DTrklets,fhTriggeredTrklets,"w"); fhSPD3DEffTrklets->SetName("fhSPD3DEffTrklets"); TH1F * fhSPDOverallTrklets=(TH1F*)fhSPDZTrklets->Clone("fhSPDOverallTrklets"); fhSPDOverallTrklets->Add(fhSPD3DTrklets); TGraphAsymmErrors *fhSPDOverallEffTrklets=new TGraphAsymmErrors(fhSPDOverallTrklets,fhTriggeredTrklets,"w"); fhSPDOverallEffTrklets->SetName("fhSPDOverallEffTrklets"); TGraphAsymmErrors *fhTRKEffTrklets=new TGraphAsymmErrors(fhTRKTrklets,fhTriggeredTrklets,"w"); fhTRKEffTrklets->SetName("fhTRKEffTrklets"); TGraphAsymmErrors *fhTRKcEffTrklets=new TGraphAsymmErrors(fhTRKcTrklets,fhTriggeredTrklets,"w"); fhTRKcEffTrklets->SetName("fhTRKcEffTrklets"); TGraphAsymmErrors *fhTRKncEffTrklets=new TGraphAsymmErrors(fhTRKncTrklets,fhTriggeredTrklets,"w"); fhTRKncEffTrklets->SetName("fhTRKncEffTrklets"); TH1F * fhSPDOverallZreco=(TH1F*)fhSPDZZreco->Clone("fhSPDOverallZreco"); fhSPDOverallZreco->Add(fhSPD3DZreco); TGraphAsymmErrors *fhSPDEffZreco=new TGraphAsymmErrors(fhSPD3DZreco,fhSPDOverallZreco,"w"); fhSPDEffZreco->SetName("fhSPDEffZreco"); TH1F *fhEff = new TH1F("hEff","hEff",6,0.5,6.5); Int_t count=1; if(fhSPDZTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhSPDZTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPDZTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"SPDZ"); count++; if(fhSPD3DTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhSPD3DTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPD3DTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"SPD3D"); count++; if(fhSPDOverallTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhSPDOverallTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPDOverallTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"SPD Overall"); count++; if(fhTRKTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhTRKTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"TRK"); count++; if(fhTRKcTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhTRKcTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKcTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"TRKc"); count++; if(fhTRKncTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){ fhEff->Fill(count,fhTRKncTrklets->GetEntries()/fhTriggeredTrklets->GetEntries()); fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKncTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries())); } fhEff->GetXaxis()->SetBinLabel(count,"TRKnc"); count++; fhEff->Print("all"); TFile* fileEff = new TFile("VtxEff.root","recreate"); fhSPDZEffTrklets->Write(); fhSPD3DEffTrklets->Write(); fhSPDOverallEffTrklets->Write(); fhTRKEffTrklets->Write(); fhTRKcEffTrklets->Write(); fhTRKncEffTrklets->Write(); fhSPDEffZreco->Write(); fhEff->Write(); fileEff->Close(); delete fileEff; ///////////////////////////////////////// */ if (!fNtupleVertexESD){ Printf("ERROR: fNtuple not available"); return; } fNtupleVertexESD = dynamic_cast(fOutput->FindObject("fNtupleVertexESD")); return; } //_________________________________________________________________________ AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const { // On the fly reco of TPC vertex from ESD AliESDEvent* evt = (AliESDEvent*) fInputEvent; AliVertexerTracks vertexer(evt->GetMagneticField()); if(evt->GetNumberOfTracks()<500) { vertexer.SetTPCMode(); // defaults } else { vertexer.SetTPCMode(0.1,1.0,5.0,10,1,3.,0.1,1.5,3.,30.,1,1);// PbPb } Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy); Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.}; if (cov[0]<1e-10) cov[0] = 1e-2; // dummy error of 1mm if (cov[2]<1e-10) cov[2] = 1e-2; // dummy error of 1mm AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1); vertexer.SetVtxStart(initVertex); delete initVertex; if(!constr) vertexer.SetConstraintOff(); return vertexer.FindPrimaryVertex(evt); } //_________________________________________________________________________ AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const { // On the fly reco of ITS+TPC vertex from ESD // mode=0 use all tracks // mode=1 use odd-number tracks // mode=2 use even-number tracks AliESDEvent* evt = (AliESDEvent*) fInputEvent; AliVertexerTracks vertexer(evt->GetMagneticField()); if(evt->GetNumberOfTracks()<500) { vertexer.SetITSMode(); // defaults vertexer.SetMinClusters(3); // default is 5 } else { vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb } //vertexer.SetITSpureSA(); Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy); Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.}; if (cov[0]<1e-10) cov[0] = 1e-2; // dummy error of 1mm if (cov[2]<1e-10) cov[2] = 1e-2; // dummy error of 1mm AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1); vertexer.SetVtxStart(initVertex); delete initVertex; if(!constr) vertexer.SetConstraintOff(); // use only ITS-TPC or only ITS-SA tracks if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) { Int_t iskip=0; Int_t ntracks = evt->GetNumberOfTracks(); Int_t *skip = new Int_t[ntracks]; for(Int_t i=0;iGetTrack(itr); if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA skip[iskip++]=itr; continue; } if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC skip[iskip++]=itr; continue; } if(mode==1 && itr%2==0) skip[iskip++]=itr; if(mode==2 && itr%2==1) skip[iskip++]=itr; } vertexer.SetSkipTracks(iskip,skip); delete [] skip; skip=NULL; } return vertexer.FindPrimaryVertex(evt); }