PWG4/JetTasks/AliAnalysisTaskJetResponse.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TMath.h"
   4 #include "TH1F.h"
   5 #include "TH2F.h"
   6 #include "TCanvas.h"
   7
   8 #include "AliLog.h"
   9
  10 #include "AliAnalysisTask.h"
  11 #include "AliAnalysisManager.h"
  12
  13 #include "AliVEvent.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliCentrality.h"
  17 #include "AliAnalysisHelperJetTasks.h"
  18 #include "AliInputEventHandler.h"
  19
  20 #include "AliAODEvent.h"
  21 #include "AliAODJet.h"
  22
  23 #include "AliAnalysisTaskJetResponse.h"
  24
  25 ClassImp(AliAnalysisTaskJetResponse)
  26
  27 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
  28   AliAnalysisTaskSE(),
  29   fESD(0x0),
  30   fAOD(0x0),
  31   fOfflineTrgMask(AliVEvent::kAny),
  32   fMinContribVtx(1),
  33   fVtxZMin(-8.),
  34   fVtxZMax(8.),
  35   fEvtClassMin(0),
  36   fEvtClassMax(10),
  37   fCentMin(0.),
  38   fCentMax(100.),
  39   fJetEtaMin(-.5),
  40   fJetEtaMax(.5),
  41   fJetDeltaEta(.4),
  42   fJetDeltaPhi(.4),
  43   fkNbranches(2),
  44   fkEvtClasses(10),
  45   fOutputList(0x0),
  46   fHistEvtSelection(0x0),
  47   fHistPtLeadingJet(new TH1F*[fkNbranches]),
  48   fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
  49   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
  50   fHistPtPtExtra(0x0),
  51   fHistPtResponse(new TH2F*[fkEvtClasses]),
  52   fHistPtSmearing(new TH2F*[fkEvtClasses])
  53 {
  54   // default Constructor
  55
  56   fJetBranchName[0] = "";
  57   fJetBranchName[1] = "";
  58
  59   fListJets[0] = new TList;
  60   fListJets[1] = new TList;
  61 }
  62
  63 AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
  64   AliAnalysisTaskSE(name),
  65   fESD(0x0),
  66   fAOD(0x0),
  67   fOfflineTrgMask(AliVEvent::kAny),
  68   fMinContribVtx(1),
  69   fVtxZMin(-8.),
  70   fVtxZMax(8.),
  71   fEvtClassMin(0),
  72   fEvtClassMax(10),
  73   fCentMin(0.),
  74   fCentMax(100.),
  75   fJetEtaMin(-.5),
  76   fJetEtaMax(.5),
  77   fJetDeltaEta(.4),
  78   fJetDeltaPhi(.4),
  79   fkNbranches(2),
  80   fkEvtClasses(10),
  81   fOutputList(0x0),
  82   fHistEvtSelection(0x0),
  83   fHistPtLeadingJet(new TH1F*[fkNbranches]),
  84   fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
  85   fHistDeltaEtaDeltaPhiLeadingJet(0x0),
  86   fHistPtPtExtra(0x0),
  87   fHistPtResponse(new TH2F*[fkEvtClasses]),
  88   fHistPtSmearing(new TH2F*[fkEvtClasses])
  89 {
  90   // Constructor
  91
  92   fJetBranchName[0] = "";
  93   fJetBranchName[1] = "";
  94
  95   fListJets[0] = new TList;
  96   fListJets[1] = new TList;
  97
  98   DefineOutput(1, TList::Class());
  99 }
 100
 101 AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
 102 {
 103   delete fListJets[0];
 104   delete fListJets[1];
 105 }
 106
 107 void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
 108 {
 109   fJetBranchName[0] = branch1;
 110   fJetBranchName[1] = branch2;
 111 }
 112
 113 void AliAnalysisTaskJetResponse::Init()
 114 {
 115
 116    // check for jet branches
 117    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 118       AliError("Jet branch name not set.");
 119    }
 120
 121 }
 122
 123 void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
 124 {
 125   // Create histograms
 126   // Called once
 127
 128   fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
 129   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 130   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 131   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 132   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 133   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 134
 135   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 136     fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
 137                                           Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
 138                                           300, 0., 300.);
 139     fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
 140
 141     fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
 142                                               Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
 143                                               100, -2., 2., 100, 0., 2*TMath::Pi());
 144
 145   }
 146   fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
 147                                              "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
 148                                              100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());
 149
 150   fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
 151                             300, 0., 300., 300, 0., 300.);
 152
 153   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 154     fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 155                                           300, 0., 300., 300, 0., 300.);
 156     fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
 157                                           200, -50., 150., 300, 0., 300.);
 158   }
 159
 160   OpenFile(1);
 161   fOutputList = new TList;
 162   fOutputList->Add(fHistEvtSelection);
 163   for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
 164     fOutputList->Add(fHistPtLeadingJet[iBranch]);
 165     fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
 166   }
 167   fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
 168   fOutputList->Add(fHistPtPtExtra);
 169   for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
 170     fOutputList->Add(fHistPtResponse[iEvtClass]);
 171     fOutputList->Add(fHistPtSmearing[iEvtClass]);
 172   }
 173 }
 174
 175 void AliAnalysisTaskJetResponse::UserExec(Option_t *)
 176 {
 177   // load events, apply event cuts, then compare leading jets
 178
 179   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 180       AliError("Jet branch name not set.");
 181       return;
 182    }
 183
 184   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 185   if (!fESD) {
 186     AliError("ESD not available");
 187     return;
 188   }
 189   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 190   if (!fAOD) {
 191     AliError("AOD not available");
 192     return;
 193   }
 194
 195   fHistEvtSelection->Fill(1); // number of events before event selection
 196
 197   // physics selection
 198   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 199     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 200   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 201     Printf(" Trigger Selection: event REJECTED ... ");
 202     fHistEvtSelection->Fill(2);
 203     PostData(1, fOutputList);
 204     return;
 205   }
 206
 207   // vertex selection
 208   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 209   Int_t nTracksPrim = primVtx->GetNContributors();
 210   if ((nTracksPrim < fMinContribVtx) ||
 211       (primVtx->GetZ() < fVtxZMin) ||
 212       (primVtx->GetZ() > fVtxZMax) ){
 213     Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 214     fHistEvtSelection->Fill(3);
 215     PostData(1, fOutputList);
 216     return;
 217   }
 218
 219   // event class selection (from jet helper task)
 220   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 221   Printf("Event class %d", eventClass);
 222   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 223     fHistEvtSelection->Fill(4);
 224     PostData(1, fOutputList);
 225     return;
 226   }
 227
 228   // centrality selection
 229   AliCentrality *cent = fESD->GetCentrality();
 230   Float_t centValue = cent->GetCentralityPercentile("TRK");
 231   printf("centrality: %f\n", centValue);
 232   if (centValue < fCentMin || centValue > fCentMax){
 233     fHistEvtSelection->Fill(5);
 234     PostData(1, fOutputList);
 235     return;
 236   }
 237
 238   fHistEvtSelection->Fill(0); // accepted events
 239
 240
 241   TClonesArray *aodJets[2];
 242   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data()));
 243   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));
 244   AliAODJet *leadingJet[2] = { 0x0, 0x0 };
 245
 246   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 247     fListJets[iJetType]->Clear();
 248     if (!aodJets[iJetType])
 249       continue;
 250     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 251       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 252       if (jet)
 253         fListJets[iJetType]->Add(jet);
 254     }
 255     fListJets[iJetType]->Sort();
 256     leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
 257
 258     if (leadingJet[iJetType])
 259       fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
 260   }
 261
 262   // check if leading jets are close in eta-phi and compare their Pt
 263   if (leadingJet[0] && leadingJet[1]) {
 264
 265     // leading jets in "jet" acceptance
 266     if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
 267        leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
 268        Printf("Jet not in eta acceptance.");
 269     }
 270     else{
 271        // check association of jets
 272        Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
 273        Float_t deltaPhi = leadingJet[0]->Phi() - leadingJet[1]->Phi();
 274        fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
 275
 276        if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
 277           printf("leading jets two far apart\n");
 278        else {
 279           fHistPtLeadingJet[0]->Fill(leadingJet[0]->Pt());
 280           fHistPtLeadingJet[1]->Fill(leadingJet[1]->Pt());
 281
 282           fHistPtPtExtra->Fill(leadingJet[0]->Pt(), leadingJet[1]->Pt());
 283
 284           if (eventClass > -1 && eventClass < fkEvtClasses){
 285               fHistPtResponse[eventClass]->Fill(leadingJet[1]->Pt(), leadingJet[0]->Pt());
 286               fHistPtSmearing[eventClass]->Fill(leadingJet[1]->Pt()-leadingJet[0]->Pt(), leadingJet[0]->Pt());
 287           }
 288        }
 289     }
 290   }
 291
 292   PostData(1, fOutputList);
 293 }
 294
 295 void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
 296 {
 297   // Draw result to the screen
 298   // Called once at the end of the query
 299
 300   if (!GetOutputData(1))
 301     return;
 302 }