PWGGA/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetSpectra.cxx

   1 // $Id$
   2 //
   3 // Jet spectrum task.
   4 //
   5 // Author: R.Reed, M.Connors
   6
   7 #include "AliAnalysisTaskEmcalJetSpectra.h"
   8
   9 #include <TCanvas.h>
  10 #include <TChain.h>
  11 #include <TClonesArray.h>
  12 #include <TH1F.h>
  13 #include <TH2F.h>
  14 #include <TList.h>
  15 #include <TLorentzVector.h>
  16 #include <TParameter.h>
  17 #include <TParticle.h>
  18 #include <TTree.h>
  19 #include <TVector3.h>
  20
  21 #include "AliAODEvent.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliAnalysisTask.h"
  24 #include "AliCentrality.h"
  25 #include "AliESDEvent.h"
  26 #include "AliESDInputHandler.h"
  27 #include "AliEmcalJet.h"
  28 #include "AliVCluster.h"
  29
  30 ClassImp(AliAnalysisTaskEmcalJetSpectra)
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra() :
  34   AliAnalysisTaskSE(),
  35   fTracksName("tracks"),
  36   fJetsName("jets"),
  37   fClustersName("clusters"),
  38   fRhos1Name("fArrRhos"),
  39   fRhos2Name("fArrRhos"),
  40   fRhos3Name("fArrRhos"),
  41   fPhimin(-10),
  42   fPhimax(10),
  43   fEtamin(-0.9),
  44   fEtamax(0.9),
  45   fAreacut(0.0),
  46   fESD(0),
  47   fOutputList(0),
  48   fHistCentrality(0),
  49   fHistDeltaRho12vsCent(0),
  50   fHistDeltaRho13vsCent(0),
  51   fHistDeltaRho23vsCent(0),
  52   fHistDeltaJetPt12vsCent(0),
  53   fHistDeltaJetPt13vsCent(0),
  54   fHistDeltaJetPt23vsCent(0),
  55   fHistRho1vsCent(0),
  56   fHistRho2vsCent(0),
  57   fHistRho3vsCent(0)
  58 {
  59   // Default constructor.
  60
  61   for (Int_t i = 0;i<6;++i) {
  62     fHistRawJetPt[i]       = 0;
  63     fHistAreavsRawPt[i]    = 0;
  64     for (Int_t j = 0;j<4;++j) {
  65       fHistNEFvsPt[i][j]   = 0;
  66        fHistZvsPt[i][j]     = 0;
  67        fHistZchvsPt[i][j]   = 0;
  68        fHistZemvsPt[i][j]   = 0;
  69        fHistJetPt[i][j]     = 0;
  70        fHistNconsvsPt[i][j] = 0;
  71        fHistJetPt3[i][j]    = 0;
  72        fHistJetPt5[i][j]    = 0;
  73        fHistJetPt7[i][j]    = 0;
  74        fHistJetPt9[i][j]    = 0;
  75     }
  76   }
  77 }
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskEmcalJetSpectra::AliAnalysisTaskEmcalJetSpectra(const char *name) :
  81   AliAnalysisTaskSE(name),
  82   fTracksName("tracks"),
  83   fJetsName("jets"),
  84   fClustersName("clusters"),
  85   fRhos1Name("fArrRhos"),
  86   fRhos2Name("fArrRhos"),
  87   fRhos3Name("fArrRhos"),
  88   fPhimin(-10),
  89   fPhimax(10),
  90   fEtamin(-0.9),
  91   fEtamax(0.9),
  92   fAreacut(0.0),
  93   fESD(0),
  94   fOutputList(0),
  95   fHistCentrality(0),
  96   fHistDeltaRho12vsCent(0),
  97   fHistDeltaRho13vsCent(0),
  98   fHistDeltaRho23vsCent(0),
  99   fHistDeltaJetPt12vsCent(0),
 100   fHistDeltaJetPt13vsCent(0),
 101   fHistDeltaJetPt23vsCent(0),
 102   fHistRho1vsCent(0),
 103   fHistRho2vsCent(0),
 104   fHistRho3vsCent(0)
 105 {
 106   // Constructor
 107
 108   for (Int_t i = 0;i<6;++i) {
 109     fHistRawJetPt[i]       = 0;
 110     fHistAreavsRawPt[i]    = 0;
 111     for (Int_t j = 0;j<4;++j) {
 112       fHistNEFvsPt[i][j]   = 0;
 113       fHistZvsPt[i][j]     = 0;
 114       fHistZchvsPt[i][j]   = 0;
 115       fHistZemvsPt[i][j]   = 0;
 116       fHistJetPt[i][j]     = 0;
 117       fHistNconsvsPt[i][j] = 0;
 118       fHistJetPt3[i][j]    = 0;
 119       fHistJetPt5[i][j]    = 0;
 120       fHistJetPt7[i][j]    = 0;
 121       fHistJetPt9[i][j]    = 0;
 122     }
 123   }
 124
 125   DefineInput(0, TChain::Class());
 126   DefineOutput(1, TList::Class());
 127 }
 128
 129 //________________________________________________________________________
 130 void AliAnalysisTaskEmcalJetSpectra::UserCreateOutputObjects()
 131 {
 132   // Called once at the beginning of the analysis.
 133
 134   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 135   if (!handler) {
 136     AliError("Input handler not available!");
 137     return;
 138   }
 139
 140   OpenFile(1);
 141   fOutputList = new TList();
 142   fOutputList->SetOwner();
 143
 144   fHistCentrality            = new TH1F("Centrality",             "Centrality",            101, -1,  100);
 145   fHistDeltaRho12vsCent      = new TH2F("DeltaRho12vsCent",       "DeltaRho12vsCent",      100, 0.0, 100.0, 500,-250,250);
 146   fHistDeltaRho13vsCent      = new TH2F("DeltaRho13vsCent",       "DeltaRho13vsCent",      100, 0.0, 100.0, 500,-250,250);
 147   fHistDeltaRho23vsCent      = new TH2F("DeltaRho23vsCent",       "DeltaRho23vsCent",      100, 0.0, 100.0, 500,-250,250);
 148   fHistDeltaJetPt12vsCent    = new TH2F("DeltaJetPt12vsCent",     "DeltaJetPt12vsCent",    100, 0.0, 100.0, 500,-250,250);
 149   fHistDeltaJetPt13vsCent    = new TH2F("DeltaJetPt13vsCent",     "DeltaJetPt13vsCent",    100, 0.0, 100.0, 500,-250,250);
 150   fHistDeltaJetPt23vsCent    = new TH2F("DeltaJetPt23vsCent",     "DeltaJetPt23vsCent",    100, 0.0, 100.0, 500,-250,250);
 151   fHistRho1vsCent            = new TH2F("Rho1vsCent",             "Rho1vsCent",            100, 0.0, 100.0, 500, 0, 500);
 152   fHistRho2vsCent            = new TH2F("Rho2vsCent",             "Rho2vsCent",            100, 0.0, 100.0, 500, 0, 500);
 153   fHistRho3vsCent            = new TH2F("Rho3vsCent",             "Rho3vsCent",            100, 0.0, 100.0, 500, 0, 500);
 154
 155   for (Int_t i = 0;i<6;++i) {
 156     TString name00(Form("fHistRawJetPt_%i",i));
 157     fHistRawJetPt[i] = new TH1F(name00,name00,250,0,500);
 158     fOutputList->Add(fHistRawJetPt[i]);
 159     TString name01(Form("fHistAreavsRawPt_%i",i));
 160     fHistAreavsRawPt[i] = new TH2F(name01,name01,250,0,500,100,0,1);
 161     fOutputList->Add(fHistAreavsRawPt[i]);
 162     for (Int_t j = 0;j<4;j++) {
 163       TString name0(Form("fHistNEFvsPt_%i_%i",i,j));
 164       fHistNEFvsPt[i][j] = new TH2F(name0,name0,250,-250,250,200,0,2);
 165       fOutputList->Add(fHistNEFvsPt[i][j]);
 166       TString name1(Form("fHistZvsPt_%i_%i",i,j));
 167       fHistZvsPt[i][j] = new TH2F(name1,name1,250,-250,250,200,0,2);
 168       fOutputList->Add(fHistZvsPt[i][j]);
 169       TString name2(Form("fHistZchvsPt_%i_%i",i,j));
 170       fHistZchvsPt[i][j] = new TH2F(name2,name2,250,-250,250,200,0,2);
 171       fOutputList->Add(fHistZchvsPt[i][j]);
 172       TString name3(Form("fHistZemvsPt_%i_%i",i,j));
 173       fHistZemvsPt[i][j] = new TH2F(name3,name3,250,-250,250,200,0,2);
 174       fOutputList->Add(fHistZemvsPt[i][j]);
 175       TString name5(Form("fHistJetPt_%i_%i",i,j));
 176       fHistJetPt[i][j] = new TH1F(name5,name5,250,-250,250);
 177       fOutputList->Add(fHistJetPt[i][j]);
 178       TString name6(Form("fHistNconsvsPt_%i_%i",i,j));
 179       fHistNconsvsPt[i][j] = new TH2F(name6,name6,250,-250,250,500,0,500);
 180       fOutputList->Add(fHistNconsvsPt[i][j]);
 181       TString name7(Form("fHistJetPt3_%i_%i",i,j));
 182       fHistJetPt3[i][j] = new TH1F(name7,name7,250,-250,250);
 183       fOutputList->Add(fHistJetPt3[i][j]);
 184       TString name8(Form("fHistJetPt5_%i_%i",i,j));
 185       fHistJetPt5[i][j] = new TH1F(name8,name8,250,-250,250);
 186       fOutputList->Add(fHistJetPt5[i][j]);
 187       TString name9(Form("fHistJetPt7_%i_%i",i,j));
 188       fHistJetPt7[i][j] = new TH1F(name9,name9,250,-250,250);
 189       fOutputList->Add(fHistJetPt7[i][j]);
 190       TString name10(Form("fHistJetPt9_%i_%i",i,j));
 191       fHistJetPt9[i][j] = new TH1F(name10,name10,250,-250,250);
 192       fOutputList->Add(fHistJetPt9[i][j]);
 193     }
 194   }
 195   fOutputList->Add(fHistCentrality);
 196   fOutputList->Add(fHistDeltaRho12vsCent);
 197   fOutputList->Add(fHistDeltaRho13vsCent);
 198   fOutputList->Add(fHistDeltaRho23vsCent);
 199   fOutputList->Add(fHistDeltaJetPt12vsCent);
 200   fOutputList->Add(fHistDeltaJetPt13vsCent);
 201   fOutputList->Add(fHistDeltaJetPt23vsCent);
 202   fOutputList->Add(fHistRho1vsCent);
 203   fOutputList->Add(fHistRho2vsCent);
 204   fOutputList->Add(fHistRho3vsCent);
 205
 206   PostData(1, fOutputList);
 207 }
 208
 209 //________________________________________________________________________
 210 Int_t AliAnalysisTaskEmcalJetSpectra::GetCentBin(Double_t cent) const
 211 {
 212   // Get centrality bin.
 213
 214   Int_t centbin = -1;
 215   if (cent>=0 && cent<10)
 216     centbin = 0;
 217   else if (cent>=10 && cent<20)
 218     centbin = 1;
 219   else if (cent>=20 && cent<30)
 220     centbin = 2;
 221   else if (cent>=30 && cent<40)
 222     centbin = 3;
 223   else if (cent>=40 && cent<50)
 224     centbin = 4;
 225   else if (cent>=50 && cent<90)
 226     centbin = 5;
 227   return centbin;
 228 }
 229
 230 //________________________________________________________________________
 231 void AliAnalysisTaskEmcalJetSpectra::UserExec(Option_t *)
 232 {
 233   // Main loop, called for each event.
 234
 235   // esd or aod mode
 236   Bool_t esdMode = kTRUE;
 237   if (dynamic_cast<AliAODEvent*>(InputEvent()))
 238     esdMode = kFALSE;
 239
 240   if (esdMode) {
 241     // optimization in case autobranch loading is off
 242     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 243     if (fTracksName == "Tracks")
 244       am->LoadBranch("Tracks");
 245   }
 246
 247   // get centrality
 248   Double_t fCent = -1;
 249   TList *l = InputEvent()->GetList();
 250   AliCentrality *centrality = InputEvent()->GetCentrality() ;
 251   if (centrality)
 252     fCent = centrality->GetCentralityPercentile("V0M");
 253   else
 254     fCent=99; // probably pp data
 255   if (fCent<0) {
 256     AliError(Form("Centrality negative: %f", fCent));
 257     return;
 258   }
 259
 260   Int_t centbin = GetCentBin(fCent);
 261
 262   //don't want to analyze events 90-100
 263   if (centbin < 0)
 264     return;
 265
 266   TClonesArray *jets = 0;
 267   TClonesArray *tracks = 0;
 268
 269   tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
 270   if (!tracks) {
 271     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 272     return;
 273   }
 274
 275   jets = dynamic_cast<TClonesArray*>(l->FindObject(fJetsName));
 276   if (!jets) {
 277     AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
 278     return;
 279   }
 280
 281   Double_t rho1 = 0;
 282   TParameter<Double_t> *Rho1Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos1Name));
 283   if (Rho1Param)
 284     rho1 = Rho1Param->GetVal();
 285
 286   Double_t rho2 = 0;
 287   TParameter<Double_t> *Rho2Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos2Name));
 288   if (Rho2Param)
 289     rho2 = Rho2Param->GetVal();
 290
 291   Double_t rho3 = 0;
 292   TParameter<Double_t> *Rho3Param = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhos3Name));
 293   if (Rho3Param)
 294     rho3 = Rho3Param->GetVal();
 295
 296   fHistRho1vsCent->Fill(fCent,rho1);
 297   fHistRho2vsCent->Fill(fCent,rho2);
 298   fHistRho3vsCent->Fill(fCent,rho3);
 299
 300   if (( rho1>0 ) && ( rho2>0 ))
 301     fHistDeltaRho12vsCent->Fill(fCent,rho1-rho2);
 302   if (( rho1>0 ) && ( rho3>0 ))
 303     fHistDeltaRho13vsCent->Fill(fCent,rho1-rho3);
 304   if (( rho2>0 ) && ( rho3>0 ))
 305     fHistDeltaRho23vsCent->Fill(fCent,rho2-rho3);
 306   fHistCentrality->Fill(fCent);
 307
 308   Double_t lJetPt   = -500;
 309   Double_t lJetPtun = -500;
 310   const Int_t Njets = jets->GetEntries();
 311   Int_t NjetAcc = 0;
 312
 313   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 314     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(iJets));
 315     if (!jet)
 316       continue;
 317     if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
 318       continue;
 319     if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
 320       continue;
 321     fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
 322     if (jet->Area()<fAreacut)
 323       continue;
 324     //prevents 0 area jets from sneaking by when area cut == 0
 325     if (jet->Area()==0)
 326       continue;
 327     if (jet->Pt()==0)
 328       continue;
 329     if (jet->MaxTrackPt()>100)
 330       continue;
 331     fHistRawJetPt[centbin]->Fill(jet->Pt());
 332
 333     NjetAcc++;
 334     Double_t Z;
 335     Double_t jetPt1 = -500;
 336     Double_t jetPt2 = -500;
 337     Double_t jetPt3 = -500;
 338     jetPt1 = jet->Pt()-jet->Area()*rho1;
 339     jetPt2 = jet->Pt()-jet->Area()*rho2;
 340     jetPt3 = jet->Pt()-jet->Area()*rho3;
 341     if (( rho1>0 ) && ( rho2>0 ))
 342       fHistDeltaJetPt12vsCent->Fill(fCent,jetPt1-jetPt2);
 343     if (( rho1>0 ) && ( rho3>0 ))
 344       fHistDeltaJetPt13vsCent->Fill(fCent,jetPt1-jetPt3);
 345     if (( rho2>0 ) && ( rho3>0 ))
 346       fHistDeltaJetPt23vsCent->Fill(fCent,jetPt2-jetPt3);
 347
 348     if (lJetPt < jetPt2) {
 349       lJetPt= jetPt2;
 350       lJetPtun = jet->Pt();
 351     }
 352
 353     Z = jet->MaxTrackPt();
 354
 355     fHistNEFvsPt[centbin][0]->Fill(jet->Pt(),jet->NEF());
 356     fHistZvsPt[centbin][0]->Fill(jet->Pt(),Z/jet->Pt());
 357     fHistZchvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxTrackPt()/jet->Pt());
 358     fHistZemvsPt[centbin][0]->Fill(jet->Pt(),jet->MaxClusterPt()/jet->Pt());
 359     fHistJetPt[centbin][0]->Fill(jet->Pt());
 360     fHistNconsvsPt[centbin][0]->Fill(jet->Pt(),jet->N());
 361     if (jet->MaxTrackPt()>3)
 362       fHistJetPt3[centbin][0]->Fill(jet->Pt());
 363     if (jet->MaxTrackPt()>5)
 364       fHistJetPt5[centbin][0]->Fill(jet->Pt());
 365     if (jet->MaxTrackPt()>7)
 366       fHistJetPt7[centbin][0]->Fill(jet->Pt());
 367     if (jet->MaxTrackPt()>9)
 368       fHistJetPt9[centbin][0]->Fill(jet->Pt());
 369
 370     if (jetPt1!=0){
 371       fHistNEFvsPt[centbin][1]->Fill(jetPt1,jet->NEF());
 372       fHistZvsPt[centbin][1]->Fill(jetPt1,Z/jetPt1);
 373       fHistZchvsPt[centbin][1]->Fill(jetPt1,jet->MaxTrackPt()/jetPt1);
 374       fHistZemvsPt[centbin][1]->Fill(jetPt1,jet->MaxClusterPt()/jetPt1);
 375       fHistJetPt[centbin][1]->Fill(jetPt1);
 376       fHistNconsvsPt[centbin][1]->Fill(jetPt1,jet->N());
 377       if (jet->MaxTrackPt()>3)
 378         fHistJetPt3[centbin][1]->Fill(jetPt1);
 379       if (jet->MaxTrackPt()>5)
 380         fHistJetPt5[centbin][1]->Fill(jetPt1);
 381       if (jet->MaxTrackPt()>7)
 382         fHistJetPt7[centbin][1]->Fill(jetPt1);
 383       if (jet->MaxTrackPt()>9)
 384         fHistJetPt9[centbin][1]->Fill(jetPt1);}
 385
 386     if (jetPt2!=0){
 387       fHistNEFvsPt[centbin][2]->Fill(jetPt2,jet->NEF());
 388       fHistZvsPt[centbin][2]->Fill(jetPt2,Z/jetPt2);
 389       fHistZchvsPt[centbin][2]->Fill(jetPt2,jet->MaxTrackPt()/jetPt2);
 390       fHistZemvsPt[centbin][2]->Fill(jetPt2,jet->MaxClusterPt()/jetPt2);
 391       fHistJetPt[centbin][2]->Fill(jetPt2);
 392       fHistNconsvsPt[centbin][2]->Fill(jetPt2,jet->N());
 393       if (jet->MaxTrackPt()>3)
 394         fHistJetPt3[centbin][2]->Fill(jetPt2);
 395       if (jet->MaxTrackPt()>5)
 396         fHistJetPt5[centbin][2]->Fill(jetPt2);
 397       if (jet->MaxTrackPt()>7)
 398         fHistJetPt7[centbin][2]->Fill(jetPt2);
 399       if (jet->MaxTrackPt()>9)
 400         fHistJetPt9[centbin][2]->Fill(jetPt2);}
 401
 402     if (jetPt3!=0){
 403       fHistNEFvsPt[centbin][3]->Fill(jetPt3,jet->NEF());
 404       fHistZvsPt[centbin][3]->Fill(jetPt3,Z/jetPt3);
 405       fHistZchvsPt[centbin][3]->Fill(jetPt3,jet->MaxTrackPt()/jetPt3);
 406       fHistZemvsPt[centbin][3]->Fill(jetPt3,jet->MaxClusterPt()/jetPt3);
 407       fHistJetPt[centbin][3]->Fill(jetPt3);
 408       fHistNconsvsPt[centbin][3]->Fill(jetPt3,jet->N());
 409       if (jet->MaxTrackPt()>3)
 410         fHistJetPt3[centbin][3]->Fill(jetPt3);
 411       if (jet->MaxTrackPt()>5)
 412         fHistJetPt5[centbin][3]->Fill(jetPt3);
 413       if (jet->MaxTrackPt()>7)
 414         fHistJetPt7[centbin][3]->Fill(jetPt3);
 415       if (jet->MaxTrackPt()>9)
 416         fHistJetPt9[centbin][3]->Fill(jetPt3);}
 417   }
 418
 419   PostData(1, fOutputList);
 420 }
 421
 422 //________________________________________________________________________
 423 void AliAnalysisTaskEmcalJetSpectra::Terminate(Option_t *)
 424 {
 425   // Called once at the end of the analysis.
 426 }
 427
 428