// ************************************** // Task used for the correction of determiantion of reconstructed jet spectra // Compares input (gen) and output (rec) jets // ******************************************* /************************************************************************** * 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 #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "TDatabasePDG.h" #include "AliAnalysisTaskJetServices.h" #include "AliAnalysisDataContainer.h" #include "AliAnalysisDataSlot.h" #include "AliAnalysisManager.h" #include "AliJetFinder.h" #include "AliJetHeader.h" #include "AliJetReader.h" #include "AliJetReaderHeader.h" #include "AliUA1JetHeaderV1.h" #include "AliESDEvent.h" #include "AliAODEvent.h" #include "AliAODHandler.h" #include "AliAODTrack.h" #include "AliAODJet.h" #include "AliAODMCParticle.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliStack.h" #include "AliGenPythiaEventHeader.h" #include "AliJetKineReaderHeader.h" #include "AliGenCocktailEventHeader.h" #include "AliInputEventHandler.h" #include "AliPhysicsSelection.h" #include "AliAnalysisHelperJetTasks.h" ClassImp(AliAnalysisTaskJetServices) AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(), fUseAODInput(kFALSE), fUsePhysicsSelection(kFALSE), fRealData(kFALSE), fSelectionInfoESD(0), fAvgTrials(1), fVtxXMean(0), fVtxYMean(0), fVtxZMean(0), fVtxRCut(1.), fVtxZCut(8.), fPtMinCosmic(5.), fRIsolMinCosmic(3.), fMaxCosmicAngle(0.01), fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fh1PtHardTrials(0x0), fh1SelectionInfoESD(0x0), fh2TriggerCount(0x0), fh2ESDTriggerCount(0x0), fh2TriggerVtx(0x0), fh2ESDTriggerVtx(0x0), fh2ESDTriggerRun(0x0), fh2VtxXY(0x0), fh1NCosmicsPerEvent(0x0), fHistList(0x0) { fRunRange[0] = fRunRange[1] = 0; } AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name): AliAnalysisTaskSE(name), fUseAODInput(kFALSE), fUsePhysicsSelection(kFALSE), fRealData(kFALSE), fSelectionInfoESD(0), fAvgTrials(1), fVtxXMean(0), fVtxYMean(0), fVtxZMean(0), fVtxRCut(1.), fVtxZCut(8.), fPtMinCosmic(5.), fRIsolMinCosmic(3.), fMaxCosmicAngle(0.01), fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fh1PtHardTrials(0x0), fh1SelectionInfoESD(0x0), fh2TriggerCount(0x0), fh2ESDTriggerCount(0x0), fh2TriggerVtx(0x0), fh2ESDTriggerVtx(0x0), fh2ESDTriggerRun(0x0), fh2VtxXY(0x0), fh1NCosmicsPerEvent(0x0), fHistList(0x0) { fRunRange[0] = fRunRange[1] = 0; DefineOutput(1,TList::Class()); } Bool_t AliAnalysisTaskJetServices::Notify() { // // Implemented Notify() to read the cross sections // and number of trials from pyxsec.root // TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); Float_t xsection = 0; Float_t ftrials = 1; fAvgTrials = 1; if(tree){ TFile *curfile = tree->GetCurrentFile(); if (!curfile) { Error("Notify","No current file"); return kFALSE; } if(!fh1Xsec||!fh1Trials){ Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__); return kFALSE; } AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials); fh1Xsec->Fill("<#sigma>",xsection); // construct a poor man average trials Float_t nEntries = (Float_t)tree->GetTree()->GetEntries(); if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries; } return kTRUE; } void AliAnalysisTaskJetServices::UserCreateOutputObjects() { // // Create the output container // if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n"); OpenFile(1); if(!fHistList)fHistList = new TList(); Bool_t oldStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); fHistList->Add(fh1Xsec); fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); fHistList->Add(fh1Trials); const Int_t nBinPt = 100; Double_t binLimitsPt[nBinPt+1]; for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ if(iPt == 0){ binLimitsPt[iPt] = 0.0; } else {// 1.0 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5; } } fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); fHistList->Add(fh2TriggerCount); fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); fHistList->Add(fh2ESDTriggerCount); fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); fHistList->Add(fh2TriggerVtx); fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); fHistList->Add(fh2ESDTriggerVtx); fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); fHistList->Add(fh1PtHard); fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt); fHistList->Add(fh1PtHardTrials); fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info", AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5); fHistList->Add(fh1SelectionInfoESD); // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range // 3 triggers BB BE/EB EE fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5); fHistList->Add(fh2ESDTriggerRun); fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10); fHistList->Add(fh2VtxXY); // =========== Switch on Sumw2 for all histos =========== for (Int_t i=0; iGetEntries(); ++i) { TH1 *h1 = dynamic_cast(fHistList->At(i)); if (h1){ h1->Sumw2(); continue; } THnSparse *hn = dynamic_cast(fHistList->At(i)); if(hn)hn->Sumw2(); } fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.); fHistList->Add(fh1NCosmicsPerEvent), TH1::AddDirectory(oldStatus); } void AliAnalysisTaskJetServices::Init() { // // Initialization // if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n"); } void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/) { // // Execute analysis for current event // AliAODEvent *aod = 0; AliESDEvent *esd = 0; AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false fSelectionInfoESD = 0; // reset AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false if(fUseAODInput){ aod = dynamic_cast(InputEvent()); if(!aod){ Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput); return; } // fethc the header } else{ // assume that the AOD is in the general output... aod = AODEvent(); if(!aod){ Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__); return; } esd = dynamic_cast(InputEvent()); } fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNone; if(esd){ Float_t run = (Float_t)esd->GetRunNumber(); const AliESDVertex *vtxESD = esd->GetPrimaryVertex(); Float_t zvtx = -999; Float_t xvtx = -999; Float_t yvtx = -999; if(vtxESD->GetNContributors()>2){ zvtx = vtxESD->GetZ(); yvtx = vtxESD->GetY(); xvtx = vtxESD->GetX(); } Int_t iTrig = -1; if(esd->GetFiredTriggerClasses().Contains("CINT1B") ||esd->GetFiredTriggerClasses().Contains("CSMBB") ||esd->GetFiredTriggerClasses().Contains("MB1") ||esd->GetFiredTriggerClasses().Contains("CINT6B")){ iTrig = 0; fSelectionInfoESD |= AliAnalysisHelperJetTasks::kBunchBunch; } else if(esd->GetFiredTriggerClasses().Contains("CINT1A") ||esd->GetFiredTriggerClasses().Contains("CSMBA") ||esd->GetFiredTriggerClasses().Contains("CINT6A") ||esd->GetFiredTriggerClasses().Contains("CINT1C") ||esd->GetFiredTriggerClasses().Contains("CSMBC") ||esd->GetFiredTriggerClasses().Contains("CINT6C")){ // empty bunch or bunch empty iTrig = 1; fSelectionInfoESD |= AliAnalysisHelperJetTasks::kBunchEmpty; } else if(esd->GetFiredTriggerClasses().Contains("CINT1-E") ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){ iTrig = 2; fSelectionInfoESD |= AliAnalysisHelperJetTasks::kEmptyEmpty; } if(iTrig>=0){ iTrig *= 3; fh2ESDTriggerRun->Fill(run,iTrig+1); if(vtxESD->GetNContributors()>2){ fh2ESDTriggerRun->Fill(run,iTrig+2); fh2VtxXY->Fill(xvtx,yvtx); } xvtx -= fVtxXMean; yvtx -= fVtxYMean; zvtx -= fVtxZMean; Float_t r2 = xvtx *xvtx + yvtx *yvtx; if(TMath::Abs(zvtx)Fill(run,iTrig+3); } else{ fh2ESDTriggerRun->Fill(run,0); } } // Apply additional constraints Bool_t esdEventSelected = IsEventSelectedESD(esd); Bool_t esdEventPileUp = IsEventPileUpESD(esd); Bool_t esdEventCosmic = IsEventCosmicESD(esd); Bool_t aodEventSelected = IsEventSelectedAOD(aod); Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&AliVEvent::kMB); if(esdEventSelected) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kVertexIn; if(esdEventPileUp) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsPileUp; if(esdEventCosmic) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kIsCosmic; if(physicsSelection) fSelectionInfoESD |= AliAnalysisHelperJetTasks::kPhysicsSelection; // here we have all selection information, fill histogram for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){ if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i); } AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false if(esd&&aod&&fDebug){ if(esdEventSelected&&!aodEventSelected){ Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__); const AliESDVertex *vtxESD = esd->GetPrimaryVertex(); const AliAODVertex *vtxAOD = aod->GetPrimaryVertex(); Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors()); vtxESD->Print(); Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors()); vtxAOD->Print(); } } // loop over all possible triggers for esd if(esd){ const AliESDVertex *vtxESD = esd->GetPrimaryVertex(); // Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print(); TString vtxName(vtxESD->GetName()); for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){ Bool_t esdTrig = kFALSE; esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it); if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered); Bool_t cand = physicsSelection; if(cand){ fh2ESDTriggerCount->Fill(it,kSelectedALICE); } if(!fUsePhysicsSelection)cand = AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1); if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){ if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex); Float_t zvtx = vtxESD->GetZ(); if(esdEventSelected&&esdTrig){ fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn); fh2ESDTriggerVtx->Fill(it,zvtx); } } if(cand&&esdEventSelected){ fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn); fh2ESDTriggerCount->Fill(it,kSelected); AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event } } } if(aod){ const AliAODVertex *vtxAOD = aod->GetPrimaryVertex(); // Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print(); TString vtxTitle(vtxAOD->GetTitle()); for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){ Bool_t aodTrig = kFALSE; aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it); Bool_t cand = (aod->GetHeader()->GetOfflineTrigger())&(AliVEvent::kMB); if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered); if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){ if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex); Float_t zvtx = vtxAOD->GetZ(); if(aodTrig&&aodEventSelected){ fh2TriggerVtx->Fill(it,zvtx); fh2TriggerCount->Fill(it,kTriggeredVertexIn); } } } } if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry); Double_t ptHard = 0; Double_t nTrials = 1; // Trials for MC trigger fh1Trials->Fill("#sum{ntrials}",fAvgTrials); AliMCEvent* mcEvent = MCEvent(); // AliStack *pStack = 0; if(mcEvent){ AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent); if(pythiaGenHeader){ nTrials = pythiaGenHeader->Trials(); ptHard = pythiaGenHeader->GetPtHard(); int iProcessType = pythiaGenHeader->ProcessType(); // 11 f+f -> f+f // 12 f+barf -> f+barf // 13 f+barf -> g+g // 28 f+g -> f+g // 53 g+g -> f+barf // 68 g+g -> g+g if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType); if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent); fh1PtHard->Fill(ptHard); fh1PtHardTrials->Fill(ptHard,nTrials); }// if pythia gen header } // trigger selection PostData(1, fHistList); } Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){ if(!esd)return kFALSE; const AliESDVertex *vtx = esd->GetPrimaryVertex(); // We can check the type of the vertex by its name and title // if vertexer failed title is empty (default c'tor) and ncontributors is 0 // vtx name title // Tracks PrimaryVertex VertexerTracksNoConstraint // TPC TPCVertex VertexerTracksNoConstraint // SPD SPDVertex vertexer: 3D // SPD SPDVertex vertexer: Z if(vtx->GetNContributors()<3)return kFALSE; // do not want tpc only primary vertex TString vtxName(vtx->GetName()); if(vtxName.Contains("TPCVertex"))return kFALSE; Float_t zvtx = vtx->GetZ(); Float_t yvtx = vtx->GetY(); Float_t xvtx = vtx->GetX(); // here we may fill the histos for run dependence.... xvtx -= fVtxXMean; yvtx -= fVtxYMean; zvtx -= fVtxZMean; Float_t r2 = yvtx*yvtx+xvtx*xvtx; Bool_t eventSel = TMath::Abs(zvtx)IsPileupFromSPD(); } Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){ if(!esd)return kFALSE; // add track cuts for which we look for cosmics... Bool_t isCosmic = kFALSE; Int_t nTracks = esd->GetNumberOfTracks(); Int_t nCosmicCandidates = 0; for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) { AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1); if (!track1) continue; UInt_t status1 = track1->GetStatus(); //If track is ITS stand alone track, skip the track if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue; if(track1->Pt()GetTrack(iTrack2); if(!track2) continue; UInt_t status2 = track2->GetStatus(); //If track is ITS stand alone track, skip the track if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue; if(track2->Pt()GetPxPyPz(mom1); track2->GetPxPyPz(mom2); TVector3 momv1(mom1[0],mom1[1],mom1[2]); TVector3 momv2(mom2[0],mom2[1],mom2[2]); Float_t theta = (float)(momv1.Phi()-momv2.Phi()); if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi(); Float_t deltaPhi = track1->Phi()-track2->Phi(); if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi(); Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) )); if(rIsolFill((float)nCosmicCandidates); return isCosmic; } Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){ if(!aod)return kFALSE; const AliAODVertex *vtx = aod->GetPrimaryVertex(); // We can check the type of the vertex by its name and title // if vertexer failed title is empty (default c'tor) and ncontributors is 0 // vtx name title // Tracks PrimaryVertex VertexerTracksNoConstraint // TPC TPCVertex VertexerTracksNoConstraint // SPD SPDVertex vertexer: 3D // SPD SPDVertex vertexer: Z if(vtx->GetNContributors()<3)return kFALSE; // do not want tpc only primary vertex TString vtxName(vtx->GetName()); if(vtxName.Contains("TPCVertex"))return kFALSE; Float_t zvtx = vtx->GetZ(); Float_t yvtx = vtx->GetY(); Float_t xvtx = vtx->GetX(); // here we may fill the histos for run dependence.... xvtx -= fVtxXMean; yvtx -= fVtxYMean; zvtx -= fVtxZMean; Float_t r2 = yvtx*yvtx+xvtx*xvtx; Bool_t eventSel = TMath::Abs(zvtx)