PWGJE/FlavourJetTasks/AliAnalysisTaskEmcalJetHF.cxx

   1 //
   2 // Author: A Castro
   3
   4 #include "AliAnalysisTaskEmcalJetHF.h"
   5
   6 // general ROOT includes
   7 #include <TCanvas.h>
   8 #include <TChain.h>
   9 #include <TClonesArray.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TH3F.h>
  13 #include <THnSparse.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 #include <TObjArray.h>
  21
  22 // AliROOT includes
  23 #include "AliAODEvent.h"
  24 #include "AliESDEvent.h"
  25 #include "AliAnalysisManager.h"
  26 #include "AliAnalysisTask.h"
  27 #include "AliCentrality.h"
  28 #include "AliEmcalJet.h"
  29 #include "AliAODJet.h"
  30 #include "AliVCluster.h"
  31 #include "AliVTrack.h"
  32 #include <AliVEvent.h>
  33 #include <AliVParticle.h>
  34 #include "AliRhoParameter.h"
  35 #include "AliEmcalParticle.h"
  36
  37 // event handler (and pico's) includes
  38 #include <AliInputEventHandler.h>
  39 #include <AliVEventHandler.h>
  40 #include "AliESDInputHandler.h"
  41 #include "AliPicoTrack.h"
  42 #include "AliEventPoolManager.h"
  43
  44 // PID includes
  45 #include "AliPIDResponse.h"
  46 #include "AliTPCPIDResponse.h"
  47 #include "AliESDpid.h"
  48
  49 #include <AliInputEventHandler.h>
  50 #include <AliVEventHandler.h>
  51
  52 // magnetic field includes
  53 #include "TGeoGlobalMagField.h"
  54 #include "AliMagF.h"
  55
  56 using std::cout;
  57 using std::endl;
  58
  59 ClassImp(AliAnalysisTaskEmcalJetHF)
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF() :
  63   AliAnalysisTaskEmcalJet("heavyF",kFALSE),
  64   event(0),
  65   fillHist(0),
  66   //isESD(0),
  67   doGlobalPID(0),
  68   fCuts(0),
  69   fPhimin(-10), fPhimax(10),
  70   fEtamin(-0.9), fEtamax(0.9),
  71   fAreacut(0.0),
  72   fJetHIpt(50.0),
  73   fTrackPtCut(2.0),
  74   fTrackEta(0.9),
  75   fTrkQAcut(0),
  76   fPIDResponse(0x0), fTPCResponse(),
  77   fEsdtrackCutsITSTPC(),
  78   fEsdtrackCutsTPC(),
  79   fEsdtrackCutsITS(),
  80   fESD(0), fAOD(0),
  81   fHistRhovsCent(0),
  82   fHistJetPhi(0),
  83   fHistCorJetPt(0), fHistJetPt(0),
  84   fHistdEdx(0),
  85   fHistdEdxvPt(0),
  86   fHistClusE(0),
  87   fHistEovPTracks(0),
  88   fHistEovPvsPtTracks(0),
  89   fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
  90   fHistnsigelectron(0),
  91   fHistnSigElecPt(0),
  92   fHistTrackPhivEta(0),
  93   fHistClusterPhivEta(0),
  94   fHistnJetTrackvnJetClusters(0),
  95   fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
  96 {
  97   // Default constructor.
  98   for (Int_t i = 0;i<6;++i){
  99     fHistJetPtvsTrackPt[i]      = 0;
 100     fHistTrackPt[i]             = 0;
 101     fHistEP0[i]                 = 0;
 102     fHistEP0A[i]                = 0;
 103     fHistEP0C[i]                = 0;
 104     fHistEPAvsC[i]              = 0;
 105   }
 106
 107   SetMakeGeneralHistograms(kTRUE);
 108
 109 }
 110
 111 //________________________________________________________________________
 112 AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF(const char *name) :
 113   AliAnalysisTaskEmcalJet(name,kTRUE),
 114   event(0),
 115   fillHist(0),
 116   //isESD(0),
 117   doGlobalPID(0),
 118   fCuts(0),
 119   fPhimin(-10), fPhimax(10),
 120   fEtamin(-0.9), fEtamax(0.9),
 121   fAreacut(0.0),
 122   fJetHIpt(50.0),
 123   fTrackPtCut(2.0),
 124   fTrackEta(0.9),
 125   fTrkQAcut(0),
 126   fPIDResponse(0x0), fTPCResponse(),
 127   fEsdtrackCutsITSTPC(),
 128   fEsdtrackCutsTPC(),
 129   fEsdtrackCutsITS(),
 130   fESD(0), fAOD(0),
 131   fHistRhovsCent(0),
 132   fHistJetPhi(0),
 133   fHistCorJetPt(0), fHistJetPt(0),
 134   fHistdEdx(0),
 135   fHistdEdxvPt(0),
 136   fHistClusE(0),
 137   fHistEovPTracks(0),
 138   fHistEovPvsPtTracks(0),
 139   fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
 140   fHistnsigelectron(0),
 141   fHistnSigElecPt(0),
 142   fHistTrackPhivEta(0),
 143   fHistClusterPhivEta(0),
 144   fHistnJetTrackvnJetClusters(0),
 145   fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
 146 {
 147   for (Int_t i = 0;i<6;++i){
 148     fHistJetPtvsTrackPt[i]      = 0;
 149     fHistTrackPt[i]             = 0;
 150     fHistEP0[i]                 = 0;
 151     fHistEP0A[i]                = 0;
 152     fHistEP0C[i]                = 0;
 153     fHistEPAvsC[i]              = 0;
 154    }
 155    SetMakeGeneralHistograms(kTRUE);
 156
 157    DefineInput(0,TChain::Class());
 158    DefineOutput(1, TList::Class());
 159 }
 160
 161 //_______________________________________________________________________
 162 AliAnalysisTaskEmcalJetHF::~AliAnalysisTaskEmcalJetHF()
 163 {
 164   // destructor
 165   //
 166   if (fOutput) {
 167     delete fOutput;
 168     fOutput = 0;
 169   }
 170 }
 171
 172 //________________________________________________________________________
 173 void AliAnalysisTaskEmcalJetHF::UserCreateOutputObjects()
 174 {
 175   if (! fCreateHisto)
 176     return;
 177   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 178
 179   fOutput = new TList();
 180   fOutput->SetOwner(kTRUE);
 181
 182   fHistJetPhi                = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
 183   fHistJetPt                 = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
 184   fOutput->Add(fHistJetPhi);
 185   fOutput->Add(fHistJetPt);
 186
 187   fillHist = 1;
 188
 189   fHistClusE = new TH1F("NumberClustersvsEnergy","NumberClustersvsEnergy", 500, 0, 10);
 190
 191
 192   if(fillHist>0){
 193   fHistRhovsCent              = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 400, 0, 400);
 194   fHistCorJetPt                     = new TH1F("NjetvsCorrJetPt", "NjetvsCorrJetPt", 300, -100, 200);
 195   fHistdEdx                             = new TH1F("dEdxSignal", "dEdxSignal", 500, 0, 500);
 196   fHistdEdxvPt                = new TH2F("dEdxvPt", "dEdxvPt", 200, 0, 100,500, 0 ,500);
 197   fHistEovPTracks             = new TH1F("EovPTracks","EovPTracks",200,0.0,2.0);
 198   fHistEovPvsPtTracks         = new TH2F("E/p_vs_Pt","E/p_vs_Pt",200, 0 ,100 ,50, 0, 2);
 199   fHistPID                    = new TH1F("fHistPID", "Detector PID", 8, 0, 8);
 200   fHistPIDtpc                 = new TH1F("fHistPIDtpc", "TPC pid", 4, 0, 4);
 201   fHistPIDits                 = new TH1F("fHistPIDits", "ITS pid", 4, 0, 4);
 202   fHistPIDtof                 = new TH1F("fHistPIDtof", "TOF pid", 4, 0, 4);
 203   fHistnsigelectron           = new TH1F("nsig_elec_TPC","nsig_elec_TPC",500,-10,10);
 204   fHistnSigElecPt             = new TH2F("nsig_v_pt(TPC)","nsig_v_pt(TPC)",200,0,100,100,-10,10);
 205   fHistTrackPhivEta           = new TH2F("TrackPhi_v_Eta","TrackPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
 206   fHistClusterPhivEta         = new TH2F("ClusterPhi_v_Eta","ClusterPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
 207   fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
 208
 209   // PT bins used to be (2000, -100, 300)
 210   TString name;
 211   TString title;
 212
 213   // creating centrality dependent histos that don't involve Global Rho
 214   for (Int_t i = 0;i<6;++i){
 215     name = TString(Form("JetPtvsTrackPt_%i",i));
 216     title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
 217     fHistJetPtvsTrackPt[i] = new TH2F(name,title, 500, -100, 400, 100,0,100);
 218     fOutput->Add(fHistJetPtvsTrackPt[i]);
 219
 220     name = TString(Form("TrackPt_%i",i));
 221     title = TString(Form("Track pT cent bin %i",i));
 222     fHistTrackPt[i] = new TH1F(name,title,400,0,200);
 223     fOutput->Add(fHistTrackPt[i]);
 224
 225     name = TString(Form("EP0_%i",i));
 226     title = TString(Form("EP VZero cent bin %i",i));
 227     fHistEP0[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
 228     fOutput->Add(fHistEP0[i]);
 229
 230     name = TString(Form("EP0A_%i",i));
 231     title = TString(Form("EP VZero cent bin %i",i));
 232     fHistEP0A[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
 233     fOutput->Add(fHistEP0A[i]);
 234
 235     name = TString(Form("EP0C_%i",i));
 236     title = TString(Form("EP VZero cent bin %i",i));
 237     fHistEP0C[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
 238     fOutput->Add(fHistEP0C[i]);
 239
 240     name = TString(Form("EPAvsC_%i",i));
 241     title = TString(Form("EP VZero cent bin %i",i));
 242     fHistEPAvsC[i] = new TH2F(name,title,144,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
 243     fOutput->Add(fHistEPAvsC[i]);
 244
 245   }
 246
 247   fOutput->Add(fHistRhovsCent);
 248   fOutput->Add(fHistCorJetPt);
 249   fOutput->Add(fHistdEdx);
 250   fOutput->Add(fHistdEdxvPt);
 251   fOutput->Add(fHistEovPTracks);
 252   fOutput->Add(fHistEovPvsPtTracks);
 253   fOutput->Add(fHistPID);
 254   fOutput->Add(fHistPIDtpc);
 255   fOutput->Add(fHistPIDits);
 256   fOutput->Add(fHistPIDtof);
 257   fOutput->Add(fHistnsigelectron);
 258   fOutput->Add(fHistnSigElecPt);
 259   fOutput->Add(fHistTrackPhivEta);
 260   fOutput->Add(fHistClusterPhivEta);
 261   fOutput->Add(fHistnJetTrackvnJetClusters);
 262
 263   }//Fill Histograms
 264
 265     fOutput->Add(fHistClusE);
 266
 267
 268   // ****************************** PID *****************************************************
 269   // set up PID handler
 270   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 271   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 272   if(!inputHandler) {
 273     AliFatal("Input handler needed");
 274   }
 275
 276   // PID response object
 277   //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
 278   //  inputHandler->CreatePIDResponse(fIsMC);         // needed to create object, why though?
 279   fPIDResponse = inputHandler->GetPIDResponse();
 280   if (!fPIDResponse) {
 281     AliError("PIDResponse object was not created");
 282   }
 283   // ****************************************************************************************
 284   UInt_t bitcoded = 0;  // bit coded, see GetDimParamsPID() below
 285   bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
 286   fhnPIDHF = NewTHnSparseDHF("fhnPIDHF", bitcoded);
 287
 288   UInt_t bitcoded1 = 0;  // bit coded, see GetDimParamsPID() below
 289   bitcoded1 = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4;
 290   fhnJetQA = NewTHnSparseDJetQA("fhnJetQA", bitcoded1);
 291
 292   UInt_t bitcoded2 = 0;  // bit coded, see GetDimParamsPID() below
 293   bitcoded2 = 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15;
 294   fhnClusQA = NewTHnSparseDJetQA("fhnClusQA", bitcoded2);
 295
 296   UInt_t bitcoded3 = 0;  // bit coded, see GetDimParamsPID() below
 297   bitcoded3 = 1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14;
 298   fhnTrackQA = NewTHnSparseDJetQA("fhnTrackQA", bitcoded3);
 299
 300   cout << "_______________Created Sparse__________________" << endl;
 301
 302   fOutput->Add(fhnPIDHF);
 303   fOutput->Add(fhnJetQA);
 304   fOutput->Add(fhnClusQA);
 305   fOutput->Add(fhnTrackQA);
 306   fOutput->Add(fhnGlobalPID);
 307
 308
 309   PostData(1, fOutput);
 310 }
 311
 312 //________________________________________________________
 313 void AliAnalysisTaskEmcalJetHF::ExecOnce()
 314 {
 315   //  Initialize the analysis
 316   AliAnalysisTaskEmcalJet::ExecOnce();
 317
 318 } // end of ExecOnce
 319
 320 //________________________________________________________________________
 321 Bool_t AliAnalysisTaskEmcalJetHF::Run()
 322 {
 323   // check to see if we have any tracks
 324   if (!fTracks)  return kTRUE;
 325   if (!fJets)  return kTRUE;
 326
 327   // what kind of event do we have: AOD or ESD?
 328   Bool_t useAOD;
 329   if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
 330   else useAOD = kFALSE;
 331
 332   // if we have ESD event, set up ESD object
 333   if(!useAOD){
 334     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 335     if (!fESD) {
 336       AliError(Form("ERROR: fESD not available\n"));
 337       return kTRUE;
 338     }
 339   }
 340
 341   // if we have AOD event, set up AOD object
 342   if(useAOD){
 343     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 344     if(!fAOD) {
 345       AliError(Form("ERROR: fAOD not available\n"));
 346       return kTRUE;
 347     }
 348   }
 349
 350   // get magnetic field info for DCA
 351   Double_t  MagF = fESD->GetMagneticField();
 352   Double_t MagSign = 1.0;
 353   if(MagF<0)MagSign = -1.0;
 354   // set magnetic field
 355   if (!TGeoGlobalMagField::Instance()->GetField()) {
 356     AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
 357     TGeoGlobalMagField::Instance()->SetField(field);
 358   }
 359
 360   // get centrality bin
 361   Int_t centbin = GetCentBin(fCent);
 362   //for pp analyses we will just use the first centrality bin
 363   if (centbin == -1)  centbin = 0;
 364
 365   // get vertex information
 366   Double_t fvertex[3]={0,0,0};
 367   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 368   //Double_t zVtx=fvertex[2];
 369
 370   // create pointer to list of input event
 371   TList *list = InputEvent()->GetList();
 372   if(!list) {
 373     AliError(Form("ERROR: list not attached\n"));
 374     return kTRUE;
 375   }
 376
 377   // background density
 378   fRhoVal = fRho->GetVal();
 379
 380   // initialize TClonesArray pointers to jets and tracks
 381   TClonesArray *jets = 0;
 382   TClonesArray *tracks = 0;
 383   TClonesArray *clusters = 0;
 384
 385   // get Jets object
 386   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 387   if(!jets){
 388     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
 389     return kTRUE;
 390   } // verify existence of jets
 391
 392   // get Tracks object
 393   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 394   if (!tracks) {
 395     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
 396     return kTRUE;
 397   } // verify existence of tracks  (fTracksName.Data())
 398
 399   //get Clusters object
 400   clusters = dynamic_cast<TClonesArray*>(list->FindObject(fCaloClusters));
 401   if (!clusters){
 402     AliError(Form("Pointer to tracks %s == 0",fCaloClusters->GetName()));
 403     return kTRUE;
 404   }  //verify cluster existence
 405
 406
 407   //const Int_t nclus = clusters->GetEntries();
 408   /* Globalcluster
 409   for(int icluster = 0; icluster < nclus; icluster++) {
 410     AliVCluster* andyCluster = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
 411     if (!icluster) continue;
 412   }
 413   */
 414   // get all emcal clusters
 415   TRefArray* caloClusters = new TRefArray();
 416   fESD->GetEMCALClusters( caloClusters );
 417
 418   //TObjArray* listcuts = fEsdtrackCutsTPC->GetAcceptedTracks(fESD);
 419   //Int_t nGoodTracks = list->GetEntries();
 420   Int_t nCluster = caloClusters->GetEntries();
 421
 422   // event plane histograms filled
 423   if(fillHist>0) fHistEP0[centbin]->Fill(fEPV0);
 424   if(fillHist>0) fHistEP0A[centbin]->Fill(fEPV0A);
 425   if(fillHist>0) fHistEP0C[centbin]->Fill(fEPV0C);
 426   if(fillHist>0) fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
 427
 428   event++;
 429   cout<<"Event #: "<<event<<endl;
 430
 431   // get number of jets and tracks
 432   const Int_t Njets = jets->GetEntries();
 433   const Int_t Ntracks = tracks->GetEntries();
 434   //const Int_t Nclusters = clusters->GetEntries();
 435   if(Ntracks<1)   return kTRUE;
 436   if(Njets<1)     return kTRUE;
 437   //if(Nclusters<1)  return kTRUE;
 438
 439   if(doGlobalPID){
 440   // loop over Globap tracks in the event for PID: Runs ON ESD's
 441   for (int i = 0;i<Ntracks;i++){
 442     //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
 443     AliVTrack *trackGlobal = static_cast<AliVTrack*>(fTracks->At(i));
 444     //AliESDtrack* track = fESD->GetTrack(i);
 445     if (! trackGlobal) {
 446       AliError(Form("Couldn't get ESD track %d\n", i));
 447       continue;
 448     }
 449
 450
 451     // apply track cuts
 452     if(TMath::Abs(trackGlobal->Eta())>fTrackEta) continue;
 453     if (trackGlobal->Pt()<0.15) continue;
 454
 455     // initialize track info
 456     Double_t dEdxGlobal = -99;
 457     Double_t trackphiGlobal = -99;
 458     Double_t trackptGlobal = 0;
 459     Double_t p = trackGlobal->P();
 460     Double_t fClsEGlobal = -99;
 461     Double_t EovPGlobal = -99;
 462
 463     // distance of closest approach
 464     Float_t DCAxyGlobal = -999;
 465     Float_t DCAzGlobal = -999;
 466     //Double_t DCAXYGlobal = -999;
 467     //Double_t DCAZGlobal = -999;
 468
 469     // track info
 470     trackGlobal->GetTPCsignal();
 471     //trackGlobal->GetImpactParameters(DCAxyGlobal, DCAzGlobal);
 472     dEdxGlobal = trackGlobal->GetTPCsignal();
 473     trackptGlobal = trackGlobal->Pt();
 474     trackphiGlobal = trackGlobal->Phi();
 475
 476     // fill track histo's
 477     //if(fillHist>0) fHistdEdx->Fill(dEdx);
 478     //if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
 479     //if(fillHist>0) fHistTrackPt[centbin]->Fill(track->Pt());
 480
 481     // clusters
 482     Int_t nMatchClusGlobal = -1;
 483     AliESDCaloCluster *matchedClusGlobal = NULL;
 484
 485     //////////////////////////////////////////////////////////////////////////
 486
 487     // cut on 1.5 GeV for EMCal Cluster                                                                                                                        //if(track->GetEMCALcluster()<0 || pt<1.5) continue;
 488     //////////////////////////////////////////////////////////////////////////
 489     // particles in TOF
 490     //Double_t nSigmaPion_TOFGlobal, nSigmaProton_TOFGlobal, nSigmaKaon_TOFGlobal = -1.;
 491
 492     // misc quantities
 493     //Double_t TOFsigGlobal  = -1.;
 494     //Double_t nClustersTPCGlobal = -1;
 495     //Int_t chargeGlobal     = 0;
 496     //Int_t trackCutsGlobal  = 0;     // set to 0 to get to run
 497     //Int_t myPIDGlboal      = 0;     // set to 0 because myPID is to compare with MC unless I hard code in cuts
 498     Double_t nSigmaElectron_TPCGlobal = -999;
 499     Double_t nSigmaElectron_TOFGlobal = -999;
 500     //Double_t nSigmaElectron_EMCALGlobal = -999;
 501
 502     // get clusters
 503     //Int_t clusiGlobal = trackGlobal->GetEMCALcluster();
 504     //nClustersTPCGlobal = trackGlobal->GetTPCclusters(0);
 505
 506     AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
 507
 508     nSigmaElectron_TPCGlobal = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
 509     nSigmaElectron_TOFGlobal = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
 510
 511     //for EMCAL
 512     nMatchClusGlobal = trackGlobal->GetEMCALcluster();
 513     if(nMatchClusGlobal > 0){
 514       matchedClusGlobal = (AliESDCaloCluster*)fESD->GetCaloCluster(nMatchClusGlobal);
 515       //AliESDCaloCluster* matchedClus = fESD->GetCaloCluster(clusi);
 516       //double eop2Global = -1;
 517       //double ssGl[4]={0.,0.,0.,0.};
 518       //Double_t nSigmaEop = fPID->GetPIDResponse()-m >NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
 519
 520       fClsEGlobal       = matchedClusGlobal->E();
 521       EovPGlobal        = fClsEGlobal/p;
 522       //nSigmaElectron_EMCALGlobal = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
 523       //if(fillHist>0) fHistEovPTracks->Fill(EovP);
 524       //if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
 525     }
 526
 527     //if(fillHist>0) fHistnsigelectron->Fill(nSigmaElectron_TPC);
 528     //if(fillHist>0) fHistnSigElecPt->Fill(trackpt,nSigmaElectron_TPC);
 529
 530     // extra attempt
 531     AliVEvent *eventQA=InputEvent();
 532     if (!eventQA||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
 533     //////////////////////////////////////////////////////////////////////////
 534         // cut on 1.5 GeV for EMCal Cluster
 535     if(trackGlobal->Pt() < fTrackPtCut) continue;
 536
 537     Double_t HF_tracks[12] = {fCent, trackGlobal->Pt(), trackGlobal->P() ,trackGlobal->Eta(), trackGlobal->Phi(), EovPGlobal, DCAxyGlobal, DCAzGlobal, dEdxGlobal, nSigmaElectron_TPCGlobal, nSigmaElectron_TOFGlobal, 0};//nSigmaElectron_EMCAL};
 538     fhnGlobalPID->Fill(HF_tracks);    // fill Sparse Histo with trigger entries
 539
 540   } // Loop over tracks
 541   }//PID Switch
 542
 543     //  Start Jet Analysis
 544     // initialize jet parameters
 545     Int_t ijethi=-1;
 546     Double_t highestjetpt=0.0;
 547
 548     // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
 549     for (Int_t ijet = 0; ijet < Njets; ijet++){
 550         // get our jets
 551         AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 552         if (!jet) continue;
 553
 554         // apply jet cuts
 555         if(!AcceptMyJet(jet)) continue;
 556
 557         Double_t JetQA[5] = {0, static_cast<Double_t>(jet->GetNumberOfTracks()), static_cast<Double_t>(jet->GetNumberOfClusters()),jet->Eta(), jet->Phi()};
 558         fhnJetQA->Fill(JetQA);
 559
 560
 561         // Loop over clusters for JetQA
 562         for(int iCluster = 0; iCluster < nCluster; iCluster++) {
 563           AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
 564           //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(jet->GetNumberOfClusters());
 565           //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
 566           //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
 567           if (!caloCluster){
 568             AliError(Form("ERROR: Could not get cluster %d", iCluster));
 569             continue;
 570           }
 571           if(!IsJetCluster(jet,iCluster,kFALSE)) continue ;
 572           Float_t pos[3];
 573           caloCluster->GetPosition(pos);  // Get cluster position
 574           TVector3 cp(pos);
 575           Double_t NtrMatched = -999.0;
 576           NtrMatched = caloCluster->GetNTracksMatched();
 577
 578           //loop over tracks for Jet QA
 579           for(int itrack = 0; itrack < NtrMatched; itrack++){
 580             //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
 581             AliVTrack *trackcluster = static_cast<AliVTrack*>(fTracks->At(itrack));
 582             //AliESDtrack* track = fESD->GetTrack(i);
 583             if (! trackcluster) {
 584               AliError(Form("Couldn't get ESD track %d\n", itrack));
 585               continue;
 586             }
 587             if(!IsJetTrack(jet,itrack,kFALSE)) continue;
 588             Double_t ClusQA[6] = {caloCluster->E(),cp.PseudoRapidity() ,cp.Phi(), static_cast<Double_t>(caloCluster->IsEMCAL()), NtrMatched, trackcluster->Phi()};
 589             fhnClusQA->Fill(ClusQA); //,1./Njets);    // fill Sparse Histo with trigger entries
 590           }//loop over tracks for JetQA
 591
 592         }//loop over clusters for JetQA
 593
 594       // loop over tracks in the event
 595       for (int iTrack = 0; iTrack<jet->GetNumberOfTracks(); iTrack++){  //loop over tracks in jet for TrackQA
 596         //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
 597         AliVTrack *jetTrack = static_cast<AliVTrack*>(fTracks->At(iTrack));
 598         //AliESDtrack* track = fESD->GetTrack(i);
 599         if (! jetTrack) {
 600           AliError(Form("Couldn't get ESD track %d\n", iTrack));
 601           continue;
 602         }
 603         if(!IsJetTrack(jet,iTrack,kFALSE)) continue;
 604         Double_t trackQA[5] = {jetTrack->Pt(), jetTrack->Eta(), jetTrack->Phi(), static_cast<Double_t>(jetTrack->IsEMCAL()), static_cast<Double_t>(jetTrack->GetEMCALcluster())};
 605         fhnTrackQA->Fill(trackQA); //,1./Njets);
 606
 607       }//track loop for TrackQA
 608
 609         if(highestjetpt<jet->Pt()){
 610             ijethi=ijet;
 611             highestjetpt=jet->Pt();
 612         }
 613     } // end of looping over jets
 614
 615   // loop over jets in the event and make appropriate cuts
 616   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 617      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 618      if (!jet)  // see if we have a jet
 619        continue;
 620
 621      // phi of jet, constrained to 1.6 < Phi < 2.94
 622      float jetphi = jet->Phi();      // phi of jet
 623
 624     // apply jet cuts
 625     if(!AcceptMyJet(jet)) continue;
 626
 627     Int_t JetClusters = jet->GetNumberOfClusters();
 628     Int_t JetTracks = jet -> GetNumberOfTracks();
 629     fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
 630
 631      // Initializations and Calculations
 632      //Double_t jeteta = jet->Eta();                                    // ETA of jet
 633      Double_t jetptraw = jet->Pt();                             // raw pT of jet
 634      Double_t jetPt = -500;                                     // initialize corr jet pt LOCAL
 635      Double_t jetarea = -500;                                   // initialize jet area
 636      jetarea = jet->Area();                                     // jet area
 637      jetPt = jet->Pt() - jetarea*fRhoVal;                 // semi-corrected pT of jet from GLOBAL rho value
 638
 639      // make histo's
 640      if(fillHist>0) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
 641      fHistJetPhi->Fill(jetphi);
 642      if(fillHist>0) fHistCorJetPt->Fill(jetPt);
 643      fHistJetPt->Fill(jetptraw);
 644
 645       if(jet->Pt() > fJetHIpt) {
 646
 647       //loop over clusters
 648       //for (int i = 0; i < njetclusters; i++){
 649       for(int iCluster = 0; iCluster < nCluster; iCluster++) {
 650           AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
 651           //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(iCluster);
 652           //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
 653           //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
 654           if (!caloCluster){
 655               AliError(Form("ERROR: Could not get cluster %d", iCluster));
 656               continue;
 657           }
 658
 659           //if (!caloCluster -> IsEMCAL()) continue; //Check that Cluster is EMCal Cluster
 660           if(! IsJetCluster(jet, iCluster, kFALSE)) continue;
 661
 662           AliESDtrack *track = 0;
 663           if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
 664           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
 665           Double_t fclusE = -999;
 666           fclusE = caloCluster->E();
 667
 668           if (fclusE<0.) continue;  //Check that cluster has positive energy
 669
 670           fHistClusE->Fill(fclusE);
 671           //Int_t  cluslabel = caloCluster->GetID();
 672           //Int_t nmatched = caloCluster->GetNTracksMatched();
 673
 674           Float_t pos[3];
 675           caloCluster->GetPosition(pos);  // Get cluster position
 676           TVector3 cp(pos);
 677           if(fillHist>0) fHistClusterPhivEta->Fill(cp.PseudoRapidity(),cp.Phi());
 678
 679           Int_t trackMatchedindex=caloCluster->GetTrackMatchedIndex();
 680           if(trackMatchedindex<0)continue;  // Make sure we don't have a bad index
 681
 682           if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
 683           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
 684           Double_t NtrMatched = -999;
 685           NtrMatched = caloCluster->GetNTracksMatched();
 686           for(int itrack = 0; itrack < NtrMatched; itrack++){   // Loop over tracks matched to clusters from jets
 687             //AliVParticle *trackCluster = static_cast<AliVParticle*>(fTracks->At(i));
 688             AliVTrack *trackCluster = static_cast<AliVTrack*>(fTracks->At(itrack));
 689             //AliESDtrack* trackCluster = fESD->GetTrack(itrack);
 690             if (! trackCluster) {
 691               AliError(Form("Couldn't get ESD track %d\n", itrack));
 692               continue;
 693             }
 694             if(!IsJetTrack(jet,itrack,kFALSE)) continue;  // Check that track is still part of ith jet
 695
 696               //if (trackCluster->Phi()>3.2 || trackCluster->Phi()<1.4)cout <<"Out of Range Track!    Track Phi:  " << trackCluster->Phi() << endl;
 697
 698               //if(track->Pt() < 4.0) continue;
 699               // initialize track info
 700               Double_t dEdx = -99;
 701               //Double_t trackphi = -99;
 702               Double_t trackpt = 0;
 703               Double_t p = trackCluster->P();
 704               Double_t EovP = -99;
 705
 706               // distance of closest approach
 707               Float_t DCAxy = -999;
 708               Float_t DCAz = -999;
 709               //Double_t DCAXY = -999;
 710               //Double_t DCAZ = -999;
 711
 712               // track info
 713               //trackCluster->GetImpactParameters(DCAxy, DCAz);
 714               trackpt = trackCluster->Pt();
 715               //trackphi = track->Phi();
 716
 717               // fill track histo's
 718               if(fillHist>0) fHistTrackPhivEta->Fill(trackCluster->Eta(),trackCluster->Phi());
 719               if(fillHist>0) fHistdEdx->Fill(dEdx);
 720               if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
 721               if(fillHist>0) fHistTrackPt[centbin]->Fill(trackCluster->Pt());
 722
 723               Double_t nSigmaElectron_TPC = -999;
 724               Double_t nSigmaElectron_TOF = -999;
 725               Double_t nSigmaElectron_EMCAL = -999;
 726
 727             if(!useAOD){
 728
 729               AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
 730
 731               dEdx = trackESD->GetTPCsignal();
 732               trackESD->GetImpactParameters(DCAxy, DCAz);
 733               nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
 734               nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
 735             }
 736            if (useAOD) {
 737               AliAODTrack *trackAOD = fAOD->GetTrack(Ntracks);
 738
 739               // get detector signals
 740               dEdx = trackAOD->GetTPCsignal();
 741
 742              //trackAOD->GetImpactParameters(DCAxy,DCAz);
 743               nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kElectron);
 744               nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kElectron);
 745
 746             } // end of AOD pid
 747
 748               //double eop2 = -1;
 749               //double ss[4]={0.,0.,0.,0.};
 750               EovP        = fclusE/p;
 751               //nSigmaElectron_EMCAL = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
 752               if(fillHist>0) fHistEovPTracks->Fill(EovP);
 753               if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
 754               Double_t HF_tracks[12] = {fCent, trackCluster->Pt(), trackCluster->P() ,trackCluster->Eta(), trackCluster->Phi(), EovP, 0/*DCAxy*/, 0/*DCAz*/, dEdx, nSigmaElectron_TPC, nSigmaElectron_TOF, nSigmaElectron_EMCAL};
 755               fhnPIDHF->Fill(HF_tracks);    // fill Sparse Histo with trigger entries
 756
 757           //Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();//find the index of the matched track. This by default returns the best match
 758           //AliESDtrack *track = event->GetTrack(trackMatchedIndex);
 759           //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
 760           //Int_t matched = fESD->AcceptTrack(track);//If the track is bad, don't count it.  By default even bad tracks are accepted
 761
 762         }  //loop over tracks matched to clusters in jet
 763
 764       } // loop over jet clusters
 765
 766
 767       } // highest pt jet cut
 768   } // LOOP over JETS in event
 769
 770   return kTRUE;
 771
 772 }
 773 //________________________________________________________________________
 774 void AliAnalysisTaskEmcalJetHF::Terminate(Option_t *)
 775 {
 776   cout<<"###########################"<<endl;
 777   cout<<"####   Task Finished   ####"<<endl;
 778   cout<<"###########################"<<endl;
 779   cout<<"###########################"<<endl;
 780 } // end of terminate
 781
 782 //________________________________________________________________________
 783 Int_t AliAnalysisTaskEmcalJetHF::AcceptMyJet(AliEmcalJet *jet) {
 784   //applies all jet cuts except pt
 785   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
 786   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
 787   if (jet->Area()<fAreacut) return 0;
 788   // prevents 0 area jets from sneaking by when area cut == 0
 789   if (jet->Area()==0) return 0;
 790   //exclude jets with extremely high pt tracks which are likely misreconstructed
 791   if(jet->MaxTrackPt()>100) return 0;
 792
 793   //passed all above cuts
 794   return 1;
 795 }
 796
 797 //________________________________________________________________________
 798 Int_t AliAnalysisTaskEmcalJetHF::GetCentBin(Double_t cent) const
 799 {  // Get centrality bin.
 800
 801   Int_t centbin = -1;
 802   if (cent>=0 && cent<10)
 803     centbin = 0;
 804   else if (cent>=10 && cent<20)
 805     centbin = 1;
 806   else if (cent>=20 && cent<30)
 807     centbin = 2;
 808   else if (cent>=30 && cent<40)
 809     centbin = 3;
 810   else if (cent>=40 && cent<50)
 811     centbin = 4;
 812   else if (cent>=50 && cent<90)
 813     centbin = 5;
 814   return centbin;
 815 }
 816 //____________________________________________________________________________________________
 817 THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDHF(const char* name, UInt_t entries)
 818 {
 819   // generate new THnSparseD PID, axes are defined in GetDimParams()
 820   Int_t count = 0;
 821   UInt_t tmp = entries;
 822   while(tmp!=0){
 823     count++;
 824     tmp = tmp &~ -tmp;  // clear lowest bit
 825   }
 826
 827   TString hnTitle(name);
 828   const Int_t dim = count;
 829   Int_t nbins[dim];
 830   Double_t xmin[dim];
 831   Double_t xmax[dim];
 832
 833   Int_t i=0;
 834   Int_t c=0;
 835   while(c<dim && i<32){
 836     if(entries&(1<<i)){
 837
 838       TString label("");
 839       GetDimParamsHF(i, label, nbins[c], xmin[c], xmax[c]);
 840       hnTitle += Form(";%s",label.Data());
 841       c++;
 842     }
 843
 844     i++;
 845   }
 846   hnTitle += ";";
 847
 848   return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 849 } // end of NewTHnSparseF PID
 850
 851 THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDJetQA(const char* name, UInt_t entries)
 852 {
 853   // generate new THnSparseD JetQA, axes are defined in GetDimParamsJetQA()
 854   Int_t count = 0;
 855   UInt_t tmp = entries;
 856   while(tmp!=0){
 857     count++;
 858     tmp = tmp &~ -tmp;  // clear lowest bit
 859   }
 860
 861   TString hnTitle(name);
 862   const Int_t dim = count;
 863   Int_t nbins[dim];
 864   Double_t xmin[dim];
 865   Double_t xmax[dim];
 866
 867   Int_t i=0;
 868   Int_t c=0;
 869   while(c<dim && i<32){
 870     if(entries&(1<<i)){
 871
 872       TString label("");
 873       GetDimParamsJetQA(i, label, nbins[c], xmin[c], xmax[c]);
 874       hnTitle += Form(";%s",label.Data());
 875       c++;
 876     }
 877
 878     i++;
 879   }
 880   hnTitle += ";";
 881
 882   return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 883 } // end of NewTHnSparseF JetQA
 884
 885
 886 void AliAnalysisTaskEmcalJetHF::GetDimParamsHF(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 887 {
 888   // stores label and binning of axis for THnSparse
 889   const Double_t pi = TMath::Pi();
 890
 891   switch(iEntry){
 892
 893   case 0:
 894     label = "V0 centrality (%)";
 895     nbins = 10;
 896     xmin = 0.;
 897     xmax = 100.;
 898     break;
 899
 900   case 1:
 901     label = "Track p_{T}";
 902     nbins = 750;
 903     xmin = 0.;
 904     xmax = 75.;
 905     break;
 906
 907   case 2:
 908     label = "Track p";
 909     nbins = 750;
 910     xmin = 0.;
 911     xmax = 75.;
 912     break;
 913
 914   case 3:
 915     label = "Track Eta";
 916     nbins = 64;
 917     xmin = -1.6;
 918     xmax = 1.6;
 919     break;
 920
 921   case 4:
 922     label = "Track Phi";
 923     nbins = 72;
 924     xmin = 0;
 925     xmax = 2*pi;
 926     break;
 927
 928   case 5:
 929     label = "E/p of track";
 930     nbins = 400;
 931     xmin = 0;
 932     xmax = 4.0;
 933     break;
 934
 935  case 6:
 936     label = "DCA xy";
 937     nbins = 20;
 938     xmin = -10;
 939     xmax =  10;
 940     break;
 941
 942   case 7:
 943     label = "DCA z";
 944     nbins = 20;
 945     xmin = -10;
 946     xmax = 10;
 947     break;
 948
 949   case 8:
 950     label = "dEdX of track - TPC";
 951     nbins = 300;
 952     xmin = 0;
 953     xmax = 300;
 954     break;
 955
 956   case 9:
 957     label = "nSigma electron TPC";
 958     nbins = 50;
 959     xmin = -5;
 960     xmax = 5;
 961     break;
 962
 963    case 10:
 964     label = "nSigma electron TOF";
 965     nbins = 50;
 966     xmin = -5;
 967     xmax = 5;
 968     break;
 969
 970    case 11:
 971     label = "nSigma electron Emcal";
 972     nbins = 50;
 973     xmin = -5;
 974     xmax = 5;
 975     break;
 976
 977
 978   } // end of switch
 979 } // end of getting dim-params
 980
 981 void AliAnalysisTaskEmcalJetHF::GetDimParamsJetQA(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 982 {
 983   // stores label and binning of axis for THnSparse
 984   const Double_t pi = TMath::Pi();
 985
 986   switch(iEntry){
 987
 988     case 0:
 989       label = "number of Jets in Event";
 990       nbins = 2;
 991       xmin = 0.;
 992       xmax = 1.;
 993       break;
 994
 995     case 1:
 996       label = "number of Clusters in a Jet";
 997       nbins = 100;
 998       xmin = 0.;
 999       xmax = 100.;
1000       break;
1001
1002     case 2:
1003       label = "number of Tracks in a Jet";
1004       nbins = 100;
1005       xmin = 0.;
1006       xmax = 100.;
1007       break;
1008
1009     case 3:
1010       label = "Jet Eta";
1011       nbins = 24;
1012       xmin = -1.2;
1013       xmax = 1.2;
1014       break;
1015
1016     case 4:
1017       label = "Jet Phi";
1018       nbins = 72;
1019       xmin = 0;
1020       xmax = 2*pi;
1021       break;
1022
1023     case 5:
1024       label = "Cluster Energy";
1025       nbins = 1000;
1026       xmin = 0;
1027       xmax = 10;
1028       break;
1029
1030     case 6:
1031       label = "Cluster Eta";
1032       nbins = 24;
1033       xmin = -1.2;
1034       xmax =  1.2;
1035       break;
1036
1037     case 7:
1038       label = "Cluster Phi";
1039       nbins = 72;
1040       xmin = 0;
1041       xmax = 2*pi;
1042       break;
1043
1044     case 8:
1045       label = "Is EMCalCluster";
1046       nbins = 2;
1047       xmin = 0;
1048       xmax = 2;
1049       break;
1050
1051     case 9:
1052       label = "Number of Matched Tracks to Cluster";
1053       nbins = 60;
1054       xmin = 0;
1055       xmax = 60;
1056       break;
1057
1058     case 10:
1059       label = "Track Pt";
1060       nbins = 750;
1061       xmin = 0;
1062       xmax = 75;
1063       break;
1064
1065     case 11:
1066       label = "Track Eta";
1067       nbins = 24;
1068       xmin = -1.2;
1069       xmax = 1.2;
1070       break;
1071
1072     case 12:
1073       label= "Track Phi";
1074       nbins = 72;
1075       xmin = 0;
1076       xmax = 2*pi;
1077       break;
1078
1079     case 13:
1080       label="Is Track EMCal";
1081       nbins = 2;
1082       xmin = 0;
1083       xmax = 2;
1084       break;
1085
1086     case 14:
1087       label = "Get Track EMCal Cluster";
1088       nbins = 100;
1089       xmin = 0;
1090       xmax = 100;
1091       break;
1092
1093     case 15:
1094       label = "Track Matched Phi";
1095       nbins = 72;
1096       xmin = 0;
1097       xmax = 2*pi;
1098
1099
1100
1101
1102   } // end of switch
1103 } // end of getting dim-params
1104
1105
1106