X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PWG0%2FdNdEta%2FAlidNdEtaTask.cxx;h=ae33c25ff276fa7362060f46f340c0a4d03a2244;hp=08c56501a69871dbc188607c14ddb04b67ad42fb;hb=81be4ee8f03ca1155e282f51da4ffd8aeb3cfdcb;hpb=1a278748965b53a65a121b8b45805b66c23d6fdd diff --git a/PWG0/dNdEta/AlidNdEtaTask.cxx b/PWG0/dNdEta/AlidNdEtaTask.cxx index 08c56501a69..ae33c25ff27 100644 --- a/PWG0/dNdEta/AlidNdEtaTask.cxx +++ b/PWG0/dNdEta/AlidNdEtaTask.cxx @@ -69,8 +69,9 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) : fUseMCKine(kFALSE), fCheckEventType(kFALSE), fSymmetrize(kFALSE), + fMultAxisEta1(kFALSE), + fDiffTreatment(AliPWG0Helper::kMCFlags), fEsdTrackCuts(0), - fPhysicsSelection(0), fTriggerAnalysis(0), fdNdEtaAnalysisESD(0), fMult(0), @@ -80,6 +81,7 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) : fdNdEtaAnalysis(0), fdNdEtaAnalysisND(0), fdNdEtaAnalysisNSD(0), + fdNdEtaAnalysisOnePart(0), fdNdEtaAnalysisTr(0), fdNdEtaAnalysisTrVtx(0), fdNdEtaAnalysisTracks(0), @@ -95,7 +97,6 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) : fFiredChips(0), fTrackletsVsClusters(0), fTrackletsVsUnassigned(0), - fTriggerVsTime(0), fStats(0), fStats2(0) { @@ -211,7 +212,7 @@ void AlidNdEtaTask::CreateOutputObjects() for (Int_t i=0; i<3; ++i) { - fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6); + fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3); fPartEta[i]->Sumw2(); fOutput->Add(fPartEta[i]); } @@ -219,22 +220,15 @@ void AlidNdEtaTask::CreateOutputObjects() fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40); fOutput->Add(fEvents); - Float_t resMax = (fAnalysisMode & AliPWG0Helper::kSPD) ? 0.2 : 2; - fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, resMax); + fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2); fOutput->Add(fVertexResolution); fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi()); fOutput->Add(fPhi); - fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi()); + fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi()); fOutput->Add(fEtaPhi); - fTriggerVsTime = new TGraph; //TH1F("fTriggerVsTime", "fTriggerVsTime;trigger time;count", 100, 0, 100); - fTriggerVsTime->SetName("fTriggerVsTime"); - fTriggerVsTime->GetXaxis()->SetTitle("trigger time"); - fTriggerVsTime->GetYaxis()->SetTitle("count"); - fOutput->Add(fTriggerVsTime); - fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5); fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d"); fStats->GetXaxis()->SetBinLabel(2, "vertexer z"); @@ -246,11 +240,11 @@ void AlidNdEtaTask::CreateOutputObjects() fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 10, -0.5, 9.5); fStats2->GetXaxis()->SetBinLabel(1, "No trigger"); - fStats2->GetXaxis()->SetBinLabel(2, "Splash identification"); - fStats2->GetXaxis()->SetBinLabel(3, "No Vertex"); - fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 10"); - fStats2->GetXaxis()->SetBinLabel(5, "0 tracklets"); - fStats2->GetXaxis()->SetBinLabel(6, "V0 Veto"); + fStats2->GetXaxis()->SetBinLabel(2, "No Vertex"); + fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad"); + fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15"); + fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10"); + fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets"); fStats2->GetXaxis()->SetBinLabel(7, "Selected"); fStats2->GetYaxis()->SetBinLabel(1, "n/a"); @@ -270,9 +264,9 @@ void AlidNdEtaTask::CreateOutputObjects() if (fAnalysisMode & AliPWG0Helper::kSPD) { - fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.2, 0.2); + fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16); fOutput->Add(fDeltaPhi); - fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.2, 0.2); + fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05); fOutput->Add(fDeltaTheta); fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5); fOutput->Add(fFiredChips); @@ -311,6 +305,9 @@ void AlidNdEtaTask::CreateOutputObjects() fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode); fOutput->Add(fdNdEtaAnalysisNSD); + fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode); + fOutput->Add(fdNdEtaAnalysisOnePart); + fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode); fOutput->Add(fdNdEtaAnalysisTr); @@ -331,17 +328,7 @@ void AlidNdEtaTask::CreateOutputObjects() fOutput->Add(fEsdTrackCuts); } - if (!fPhysicsSelection) - fPhysicsSelection = new AliPhysicsSelection; - - fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug); - //fOutput->Add(fPhysicsSelection); - - fTriggerAnalysis = new AliTriggerAnalysis; - fTriggerAnalysis->EnableHistograms(); - fTriggerAnalysis->SetSPDGFOThreshhold(2); - fOutput->Add(fTriggerAnalysis); } void AlidNdEtaTask::Exec(Option_t*) @@ -365,37 +352,24 @@ void AlidNdEtaTask::Exec(Option_t*) return; } -// if (fCheckEventType) -// eventTriggered = fPhysicsSelection->IsCollisionCandidate(fESD); - - // check event type (should be PHYSICS = 7) - if (fCheckEventType) + AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); + if (!inputHandler) { - UInt_t eventType = esdHeader->GetEventType(); - if (eventType != 7) - { - Printf("Skipping event because it is of type %d", eventType); - return; - } - - //Printf("Trigger classes: %s:", fESD->GetFiredTriggerClasses().Data()); - - Bool_t accept = kTRUE; - if (fRequireTriggerClass.Length() > 0 && !fESD->IsTriggerClassFired(fRequireTriggerClass)) - accept = kFALSE; - if (fRejectTriggerClass.Length() > 0 && fESD->IsTriggerClassFired(fRejectTriggerClass)) - accept = kFALSE; - - if (!accept) - { - Printf("Skipping event because it does not have the correct trigger class(es): %s", fESD->GetFiredTriggerClasses().Data()); - return; - } - - fStats->Fill(4); + Printf("ERROR: Could not receive input handler"); + return; } - fTriggerAnalysis->FillHistograms(fESD); + eventTriggered = inputHandler->IsEventSelected(); + + if (!fTriggerAnalysis) + { + AliPhysicsSelection* physicsSelection = dynamic_cast (inputHandler->GetEventSelection()); + if (physicsSelection) + fTriggerAnalysis = physicsSelection->GetTriggerAnalysis(); + } + + if (eventTriggered) + eventTriggered = fTriggerAnalysis->IsTriggerFired(fESD, fTrigger); AliTriggerAnalysis::V0Decision v0A = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE); AliTriggerAnalysis::V0Decision v0C = fTriggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE); @@ -413,134 +387,102 @@ void AlidNdEtaTask::Exec(Option_t*) if (v0A == AliTriggerAnalysis::kV0BB && v0C == AliTriggerAnalysis::kV0BG) vZero = 8; if (v0A == AliTriggerAnalysis::kV0BG && v0C == AliTriggerAnalysis::kV0BB) vZero = 9; } + + if (vZero == 0) + Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C); Bool_t filled = kFALSE; - // trigger - eventTriggered = fTriggerAnalysis->IsTriggerFired(fESD, fTrigger); - if (!eventTriggered) { fStats2->Fill(0.0, vZero); filled = kTRUE; } - if (v0A == AliTriggerAnalysis::kV0BG || v0C == AliTriggerAnalysis::kV0BG) - eventTriggered = kFALSE; - if (eventTriggered) - { fStats->Fill(3); - } - if (fCheckEventType) + /* + // ITS cluster tree + AliESDInputHandlerRP* handlerRP = dynamic_cast (AliAnalysisManager::GetAnalys +isManager()->GetInputEventHandler()); + if (handlerRP) { - /*const Int_t kMaxEvents = 1; - UInt_t maskedEvents[kMaxEvents][2] = { - {-1, -1} - }; - - for (Int_t i=0; iGetPeriodNumber() == maskedEvents[i][0] && fESD->GetOrbitNumber() == maskedEvents[i][1]) - { - Printf("Skipping event because it is masked: period: %d orbit: %x", fESD->GetPeriodNumber(), fESD->GetOrbitNumber()); - if (!filled) + TTree* itsClusterTree = handlerRP->GetTreeR("ITS"); + if (!itsClusterTree) + return; + + TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint"); + TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints"); + + itsClusterBranch->SetAddress(&itsClusters); + + Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries(); + + Int_t totalClusters = 0; + + // loop over the its subdetectors + for (Int_t iIts=0; iIts < nItsSubs; iIts++) { + + if (!itsClusterTree->GetEvent(iIts)) + continue; + + Int_t nClusters = itsClusters->GetEntriesFast(); + totalClusters += nClusters; + + #ifdef FULLALIROOT + if (fAnalysisMode & AliPWG0Helper::kSPD) { - fStats2->Fill(1, vZero); - filled = kTRUE; - } - return; - } - }*/ - - // ITS cluster tree - AliESDInputHandlerRP* handlerRP = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); - if (handlerRP) - { - TTree* itsClusterTree = handlerRP->GetTreeR("ITS"); - if (!itsClusterTree) - return; - - TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint"); - TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints"); - - itsClusterBranch->SetAddress(&itsClusters); - - Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries(); - - Int_t totalClusters = 0; - - // loop over the its subdetectors - for (Int_t iIts=0; iIts < nItsSubs; iIts++) { - - if (!itsClusterTree->GetEvent(iIts)) - continue; - - Int_t nClusters = itsClusters->GetEntriesFast(); - totalClusters += nClusters; - - #ifdef FULLALIROOT - if (fAnalysisMode & AliPWG0Helper::kSPD) - { - // loop over clusters - while (nClusters--) { - AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters); - - Int_t layer = cluster->GetLayer(); - - if (layer > 1) - continue; - - Float_t xyz[3] = {0., 0., 0.}; - cluster->GetGlobalXYZ(xyz); - - Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]); - Float_t z = xyz[2]; - - fZPhi[layer]->Fill(z, phi); - fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex()); - } + // loop over clusters + while (nClusters--) { + AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters); + + Int_t layer = cluster->GetLayer(); + + if (layer > 1) + continue; + + Float_t xyz[3] = {0., 0., 0.}; + cluster->GetGlobalXYZ(xyz); + + Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]); + Float_t z = xyz[2]; + + fZPhi[layer]->Fill(z, phi); + fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex()); } - #endif - } - - const AliMultiplicity* mult = fESD->GetMultiplicity(); - if (!mult) - return; - - fTrackletsVsClusters->Fill(mult->GetNumberOfTracklets(), totalClusters); - - Int_t limit = 80 + mult->GetNumberOfTracklets() * 220/20; - - if (totalClusters > limit) - { - Printf("Skipping event because %d clusters is above limit of %d from %d tracklets: Period number: %d Orbit number: %x", totalClusters, limit, mult->GetNumberOfTracklets(), fESD->GetPeriodNumber(), fESD->GetOrbitNumber()); - if (!filled) - { - fStats2->Fill(1, vZero); - filled = kTRUE; } - return; // TODO we skip this also for the MC. not good... - } + #endif } } - + */ + vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode); if (!vtxESD) { if (!filled) { - fStats2->Fill(2, vZero); + fStats2->Fill(1, vZero); filled = kTRUE; } } else { + if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode)) + { + if (!filled) + { + fStats2->Fill(2, vZero); + filled = kTRUE; + } + } + Double_t vtx[3]; vtxESD->GetXYZ(vtx); - if (TMath::Abs(vtx[2]) > 10) + // try to compare spd vertex and vertexer from tracks + // remove vertices outside +- 15 cm + if (TMath::Abs(vtx[2]) > 15) { if (!filled) { @@ -548,6 +490,15 @@ void AlidNdEtaTask::Exec(Option_t*) filled = kTRUE; } } + + if (TMath::Abs(vtx[2]) > 10) + { + if (!filled) + { + fStats2->Fill(4, vZero); + filled = kTRUE; + } + } const AliMultiplicity* mult = fESD->GetMultiplicity(); if (!mult) @@ -557,26 +508,17 @@ void AlidNdEtaTask::Exec(Option_t*) { if (!filled) { - fStats2->Fill(4, vZero); + fStats2->Fill(5, vZero); filled = kTRUE; } } } - if (fCheckEventType) - { - if (vZero >= 5) - { - if (!filled) - fStats2->Fill(5, vZero); - return; - } - } - if (!filled) + { fStats2->Fill(6, vZero); - - //Printf("Skipping event %d: Period number: %d Orbit number: %x", decision, fESD->GetPeriodNumber(), fESD->GetOrbitNumber()); + //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber()); + } if (eventTriggered) fStats->Fill(3); @@ -586,14 +528,16 @@ void AlidNdEtaTask::Exec(Option_t*) // get the ESD vertex vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode); - //Printf("vertex is: %s", fESD->GetPrimaryVertex()->GetTitle()); - Double_t vtx[3]; // fill z vertex resolution if (vtxESD) { - fVertexResolution->Fill(vtxESD->GetZRes()); + if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0) + fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion()); + else + fVertexResolution->Fill(vtxESD->GetZRes(), 0); + //if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0) { fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv()); @@ -607,7 +551,7 @@ void AlidNdEtaTask::Exec(Option_t*) if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0) { fStats->Fill(1); - if (fCheckEventType && TMath::Abs(vtx[0] > 0.3)) + if (fCheckEventType && TMath::Abs(vtx[0]) > 0.3) { Printf("Suspicious x-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber()); } @@ -616,10 +560,13 @@ void AlidNdEtaTask::Exec(Option_t*) Printf("Suspicious y-vertex x=%f y=%f z=%f (period: %d orbit %x)", vtx[0], vtx[1], vtx[2], fESD->GetPeriodNumber(), fESD->GetOrbitNumber()); } } - else if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0) - { + + if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0) fStats->Fill(2); - } + + // remove vertices outside +- 15 cm + if (TMath::Abs(vtx[2]) > 15) + vtxESD = 0; } else vtxESD = 0; @@ -673,7 +620,7 @@ void AlidNdEtaTask::Exec(Option_t*) return; } } - + // create list of (label, eta, pt) tuples Int_t inputCount = 0; Int_t* labelArr = 0; @@ -681,99 +628,114 @@ void AlidNdEtaTask::Exec(Option_t*) Float_t* thirdDimArr = 0; if (fAnalysisMode & AliPWG0Helper::kSPD) { - // get tracklets - const AliMultiplicity* mult = fESD->GetMultiplicity(); - if (!mult) - { - AliDebug(AliLog::kError, "AliMultiplicity not available"); - return; - } - - labelArr = new Int_t[mult->GetNumberOfTracklets()]; - etaArr = new Float_t[mult->GetNumberOfTracklets()]; - thirdDimArr = new Float_t[mult->GetNumberOfTracklets()]; - - // get multiplicity from SPD tracklets - for (Int_t i=0; iGetNumberOfTracklets(); ++i) + if (vtxESD) { - //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); - - if (fOnlyPrimaries) - if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0))) - continue; - - Float_t deltaPhi = mult->GetDeltaPhi(i); - // prevent values to be shifted by 2 Pi() - if (deltaPhi < -TMath::Pi()) - deltaPhi += TMath::Pi() * 2; - if (deltaPhi > TMath::Pi()) - deltaPhi -= TMath::Pi() * 2; - - if (TMath::Abs(deltaPhi) > 1) - printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi); - - Int_t label = mult->GetLabel(i, 0); - Float_t eta = mult->GetEta(i); - - // control histograms - Float_t phi = mult->GetPhi(i); - if (phi < 0) - phi += TMath::Pi() * 2; - fPhi->Fill(phi); - fEtaPhi->Fill(eta, phi); - -// if (deltaPhi < 0.01) -// { -// // layer 0 -// Float_t z = vtx[2] + 3.9 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta))); -// fZPhi[0]->Fill(z, phi); -// // layer 1 -// z = vtx[2] + 7.6 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta))); -// fZPhi[1]->Fill(z, phi); -// } - - if (vtxESD && TMath::Abs(vtx[2]) < 10) + // get tracklets + const AliMultiplicity* mult = fESD->GetMultiplicity(); + if (!mult) { - fDeltaPhi->Fill(deltaPhi); - fDeltaTheta->Fill(mult->GetDeltaTheta(i)); + AliDebug(AliLog::kError, "AliMultiplicity not available"); + return; + } + + Int_t arrayLength = mult->GetNumberOfTracklets(); + if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0) + arrayLength += mult->GetNumberOfSingleClusters(); + + labelArr = new Int_t[arrayLength]; + etaArr = new Float_t[arrayLength]; + thirdDimArr = new Float_t[arrayLength]; + + // get multiplicity from SPD tracklets + for (Int_t i=0; iGetNumberOfTracklets(); ++i) + { + //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i)); + + if (fOnlyPrimaries) + if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0))) + continue; + + Float_t deltaPhi = mult->GetDeltaPhi(i); + // prevent values to be shifted by 2 Pi() + if (deltaPhi < -TMath::Pi()) + deltaPhi += TMath::Pi() * 2; + if (deltaPhi > TMath::Pi()) + deltaPhi -= TMath::Pi() * 2; + + if (TMath::Abs(deltaPhi) > 1) + printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi); + + Int_t label = mult->GetLabel(i, 0); + Float_t eta = mult->GetEta(i); + + // control histograms + Float_t phi = mult->GetPhi(i); + if (phi < 0) + phi += TMath::Pi() * 2; + + // TEST exclude probably inefficient phi region + //if (phi > 5.70 || phi < 0.06) + // continue; + + fPhi->Fill(phi); + + if (vtxESD && TMath::Abs(vtx[2]) < 10) + { + fEtaPhi->Fill(eta, phi); + fDeltaPhi->Fill(deltaPhi); + fDeltaTheta->Fill(mult->GetDeltaTheta(i)); + } + + if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025)) + continue; + + if (fUseMCKine) + { + if (label > 0) + { + TParticle* particle = stack->Particle(label); + eta = particle->Eta(); + phi = particle->Phi(); + } + else + Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found"); + } + + if (fSymmetrize) + eta = TMath::Abs(eta); + + etaArr[inputCount] = eta; + labelArr[inputCount] = label; + thirdDimArr[inputCount] = phi; + ++inputCount; } - if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut) - continue; - - if (fUseMCKine) + if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0) { - if (label > 0) + // get additional clusters from L0 + for (Int_t i=0; iGetNumberOfSingleClusters(); ++i) { - TParticle* particle = stack->Particle(label); - eta = particle->Eta(); - phi = particle->Phi(); + etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.)); + labelArr[inputCount] = -1; + thirdDimArr[inputCount] = mult->GetPhiSingle(i); + + ++inputCount; } - else - Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found"); } - if (fSymmetrize) - eta = TMath::Abs(eta); - - etaArr[inputCount] = eta; - labelArr[inputCount] = label; - thirdDimArr[inputCount] = phi; - ++inputCount; - } - - if (!fFillPhi) - { - // fill multiplicity in third axis - for (Int_t i=0; iGetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1); + fFiredChips->Fill(firedChips, inputCount); + Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips); + + fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters()); } - - Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1); - fFiredChips->Fill(firedChips, inputCount); - Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips); - - fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters()); } else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) { @@ -783,6 +745,8 @@ void AlidNdEtaTask::Exec(Option_t*) return; } + Bool_t foundInEta10 = kFALSE; + if (vtxESD) { // get multiplicity from ESD tracks @@ -821,6 +785,13 @@ void AlidNdEtaTask::Exec(Option_t*) if (eventTriggered && vtxESD) fRawPt->Fill(pT); + if (esdTrack->Pt() < 0.15) + continue; + + // INEL>0 trigger + if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15) + foundInEta10 = kTRUE; + if (fOnlyPrimaries) { if (label == 0) @@ -850,13 +821,16 @@ void AlidNdEtaTask::Exec(Option_t*) ++inputCount; } - if (inputCount > 30) - Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber()); + //if (inputCount > 30) + // Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber()); // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1 delete list; } + + if (!foundInEta10) + eventTriggered = kFALSE; } else return; @@ -868,8 +842,6 @@ void AlidNdEtaTask::Exec(Option_t*) fMult->Fill(inputCount); fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount); - fTriggerVsTime->SetPoint(fTriggerVsTime->GetN(), fESD->GetTimeStamp(), 1); - if (vtxESD) { // control hist @@ -877,8 +849,8 @@ void AlidNdEtaTask::Exec(Option_t*) if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0) fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount); - fMultVtx->Fill(inputCount); - + Int_t part05 = 0; + Int_t part10 = 0; for (Int_t i=0; iFill(eta); } + + if (TMath::Abs(eta) < 0.5) + part05++; + if (TMath::Abs(eta) < 1.0) + part10++; } + + Int_t multAxis = inputCount; + if (fMultAxisEta1) + multAxis = part10; + + fMultVtx->Fill(multAxis); + //if (TMath::Abs(vtx[2]) < 10) + // fMultVtx->Fill(part05); // for event count per vertex - fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount); + fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis); // control hist - if (inputCount > 0) + if (multAxis > 0) fEvents->Fill(vtx[2]); if (fReadMC) @@ -931,7 +916,7 @@ void AlidNdEtaTask::Exec(Option_t*) thirdDim = particle->Phi(); } else - thirdDim = inputCount; + thirdDim = multAxis; } else thirdDim = particle->Pt(); @@ -943,11 +928,10 @@ void AlidNdEtaTask::Exec(Option_t*) } // end of track loop // for event count per vertex - fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount); + fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis); } } } - if (etaArr) delete[] etaArr; if (labelArr) @@ -994,18 +978,19 @@ void AlidNdEtaTask::Exec(Option_t*) TArrayF vtxMC(3); genHeader->PrimaryVertex(vtxMC); - + // get process type - Int_t processType = AliPWG0Helper::GetEventProcessType(header); + Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment); AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType)); - if (processType==AliPWG0Helper::kInvalidProcess) + if (processType == AliPWG0Helper::kInvalidProcess) AliDebug(AliLog::kError, Form("Unknown process type %d.", processType)); // loop over mc particles Int_t nPrim = stack->GetNprimary(); Int_t nAcceptedParticles = 0; + Bool_t oneParticleEvent = kFALSE; // count particles first, then fill for (Int_t iMc = 0; iMc < nPrim; ++iMc) @@ -1025,7 +1010,10 @@ void AlidNdEtaTask::Exec(Option_t*) //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID())); Float_t eta = particle->Eta(); - // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison)) + if (TMath::Abs(eta) < 1.0) + oneParticleEvent = kTRUE; + + // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison)) if (TMath::Abs(eta) < 1.5) // && pt > 0.3) nAcceptedParticles++; } @@ -1069,6 +1057,9 @@ void AlidNdEtaTask::Exec(Option_t*) if (processType == AliPWG0Helper::kND) fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim); + + if (oneParticleEvent) + fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim); if (eventTriggered) { @@ -1090,6 +1081,8 @@ void AlidNdEtaTask::Exec(Option_t*) fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles); if (processType == AliPWG0Helper::kND) fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles); + if (oneParticleEvent) + fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles); if (eventTriggered) { @@ -1118,7 +1111,7 @@ void AlidNdEtaTask::Terminate(Option_t *) for (Int_t i=0; i<3; ++i) fPartEta[i] = dynamic_cast (fOutput->FindObject(Form("dndeta_check_%d", i))); fEvents = dynamic_cast (fOutput->FindObject("dndeta_check_vertex")); - fVertexResolution = dynamic_cast (fOutput->FindObject("dndeta_vertex_resolution_z")); + fVertexResolution = dynamic_cast (fOutput->FindObject("dndeta_vertex_resolution_z")); fVertex = dynamic_cast (fOutput->FindObject("vertex_check")); fVertexVsMult = dynamic_cast (fOutput->FindObject("fVertexVsMult")); @@ -1133,13 +1126,10 @@ void AlidNdEtaTask::Terminate(Option_t *) fFiredChips = dynamic_cast (fOutput->FindObject("fFiredChips")); fTrackletsVsClusters = dynamic_cast (fOutput->FindObject("fTrackletsVsClusters")); fTrackletsVsUnassigned = dynamic_cast (fOutput->FindObject("fTrackletsVsUnassigned")); - fTriggerVsTime = dynamic_cast (fOutput->FindObject("fTriggerVsTime")); fStats = dynamic_cast (fOutput->FindObject("fStats")); fStats2 = dynamic_cast (fOutput->FindObject("fStats2")); fEsdTrackCuts = dynamic_cast (fOutput->FindObject("fEsdTrackCuts")); - fPhysicsSelection = dynamic_cast (fOutput->FindObject("AliPhysicsSelection_outputlist")); - fTriggerAnalysis = dynamic_cast (fOutput->FindObject("AliTriggerAnalysis")); } if (!fdNdEtaAnalysisESD) @@ -1164,10 +1154,11 @@ void AlidNdEtaTask::Terminate(Option_t *) if (fPartEta[0]) { - Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001)); - Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9)); + Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001)); + Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9)); Printf("%d events with vertex used", events1 + events2); + Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1)); if (events1 > 0 && events2 > 0) { @@ -1207,15 +1198,6 @@ void AlidNdEtaTask::Terminate(Option_t *) if (fEsdTrackCuts) fEsdTrackCuts->SaveHistograms("esd_track_cuts"); - if (fPhysicsSelection) - { - fPhysicsSelection->SaveHistograms("physics_selection"); - fPhysicsSelection->Print(); - } - - if (fTriggerAnalysis) - fTriggerAnalysis->SaveHistograms(); - if (fMult) fMult->Write(); @@ -1230,7 +1212,11 @@ void AlidNdEtaTask::Terminate(Option_t *) fEvents->Write(); if (fVertexResolution) + { fVertexResolution->Write(); + fVertexResolution->ProjectionX(); + fVertexResolution->ProjectionY(); + } if (fDeltaPhi) fDeltaPhi->Write(); @@ -1263,9 +1249,6 @@ void AlidNdEtaTask::Terminate(Option_t *) if (fTrackletsVsUnassigned) fTrackletsVsUnassigned->Write(); - if (fTriggerVsTime) - fTriggerVsTime->Write(); - if (fStats) fStats->Write(); @@ -1291,6 +1274,7 @@ void AlidNdEtaTask::Terminate(Option_t *) fdNdEtaAnalysis = dynamic_cast (fOutput->FindObject("dndeta")); fdNdEtaAnalysisND = dynamic_cast (fOutput->FindObject("dndetaND")); fdNdEtaAnalysisNSD = dynamic_cast (fOutput->FindObject("dndetaNSD")); + fdNdEtaAnalysisOnePart = dynamic_cast (fOutput->FindObject("dndetaOnePart")); fdNdEtaAnalysisTr = dynamic_cast (fOutput->FindObject("dndetaTr")); fdNdEtaAnalysisTrVtx = dynamic_cast (fOutput->FindObject("dndetaTrVtx")); fdNdEtaAnalysisTracks = dynamic_cast (fOutput->FindObject("dndetaTracks")); @@ -1306,6 +1290,7 @@ void AlidNdEtaTask::Terminate(Option_t *) fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone); fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone); fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone); + fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone); fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone); fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone); fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone); @@ -1319,6 +1304,7 @@ void AlidNdEtaTask::Terminate(Option_t *) fdNdEtaAnalysis->SaveHistograms(); fdNdEtaAnalysisND->SaveHistograms(); fdNdEtaAnalysisNSD->SaveHistograms(); + fdNdEtaAnalysisOnePart->SaveHistograms(); fdNdEtaAnalysisTr->SaveHistograms(); fdNdEtaAnalysisTrVtx->SaveHistograms(); fdNdEtaAnalysisTracks->SaveHistograms();