PWGJE/FlavourJetTasks/AliAnalysisTaskEmcalJetFlavourTagExample.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /*
  17 *   This sample task demostrates the basics of tagging a jet. The PID response
  18 *   is retrived for both the TPC.  A Hadronic tag is implemented for
  19 *   clusters matched to tracks in a jet with a E/P < 0.2
  20 *
  21 *   Author: Andrew Castro (UTK) and Joel Mazer (UTK)
  22 */
  23
  24 #include "AliAnalysisTaskEmcalJetFlavourTagExample.h"
  25
  26 // general ROOT includes
  27 #include <TCanvas.h>
  28 #include <TChain.h>
  29 #include <TClonesArray.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <THnSparse.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TParameter.h>
  37 #include <TParticle.h>
  38 #include <TTree.h>
  39 #include <TVector3.h>
  40 #include <TObjArray.h>
  41
  42 // AliROOT includes
  43 #include "AliAODEvent.h"
  44 #include "AliESDEvent.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAnalysisTask.h"
  47 #include "AliCentrality.h"
  48 #include "AliEmcalJet.h"
  49 #include "AliAODJet.h"
  50 #include "AliVCluster.h"
  51 #include "AliVTrack.h"
  52 #include <AliVEvent.h>
  53 #include <AliVParticle.h>
  54 #include "AliRhoParameter.h"
  55 #include "AliLog.h"
  56 #include "AliJetContainer.h"
  57 #include "AliParticleContainer.h"
  58 #include "AliClusterContainer.h"
  59 #include "AliEmcalParticle.h"
  60 #include "AliESDCaloCluster.h"
  61 #include <AliESDtrackCuts.h>
  62 #include "AliPID.h"
  63
  64 // event handler (and pico's) includes
  65 #include <AliInputEventHandler.h>
  66 #include <AliVEventHandler.h>
  67 #include "AliESDInputHandler.h"
  68 #include "AliPicoTrack.h"
  69 #include "AliEventPoolManager.h"
  70 #include "AliAODTrack.h"
  71 #include "AliESDtrack.h"
  72
  73 // PID includes
  74 #include "AliPIDResponse.h"
  75 #include "AliTPCPIDResponse.h"
  76 #include "AliESDpid.h"
  77
  78 // magnetic field includes
  79 #include "TGeoGlobalMagField.h"
  80 #include "AliMagF.h"
  81
  82 using std::cout;
  83 using std::endl;
  84
  85 ClassImp(AliAnalysisTaskEmcalJetFlavourTagExample)
  86
  87 //________________________________________________________________________
  88 AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample() :
  89   AliAnalysisTaskEmcalJet("heavyF",kFALSE),
  90   event(0),
  91   fCuts(0),
  92   fPhimin(-10), fPhimax(10),
  93   fEtamin(-0.9), fEtamax(0.9),
  94   fAreacut(0.0),
  95   fJetHIpt(50.0),
  96   fTrackPtCut(2.0),
  97   fTrackEta(0.9),
  98   fesdTrackCuts(0),
  99   fPIDResponse(0x0), fTPCResponse(),
 100   fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
 101   fESD(0), fAOD(0),
 102   fHistJetPhi(0),
 103   fHistCorJetPt(0), fHistJetPt(0),
 104   fHistHighJetPt(0),
 105   fHistnSigElecPt(0),
 106   fHistdEdXvsPt(0),
 107   fHistnJetTrackvnJetClusters(0)
 108 {
 109   // Default constructor.
 110
 111
 112   SetMakeGeneralHistograms(kTRUE);
 113
 114 }
 115
 116 //________________________________________________________________________
 117 AliAnalysisTaskEmcalJetFlavourTagExample::AliAnalysisTaskEmcalJetFlavourTagExample(const char *name) :
 118   AliAnalysisTaskEmcalJet(name,kTRUE),
 119   event(0),
 120   fCuts(0),
 121   fPhimin(-10), fPhimax(10),
 122   fEtamin(-0.9), fEtamax(0.9),
 123   fAreacut(0.0),
 124   fJetHIpt(50.0),
 125   fTrackPtCut(2.0),
 126   fTrackEta(0.9),
 127   fesdTrackCuts(0),
 128   fPIDResponse(0x0), fTPCResponse(),
 129   fJetsCont(0), fTracksCont(0), fCaloClustersCont(0), fTracksJetCont(0), fCaloClustersJetCont(0),
 130   fESD(0), fAOD(0),
 131   fHistJetPhi(0),
 132   fHistCorJetPt(0), fHistJetPt(0),
 133   fHistHighJetPt(0),
 134   fHistnSigElecPt(0),
 135   fHistdEdXvsPt(0),
 136   fHistnJetTrackvnJetClusters(0)
 137 {
 138
 139    SetMakeGeneralHistograms(kTRUE);
 140
 141    DefineInput(0,TChain::Class());
 142    DefineOutput(1, TList::Class());
 143 }
 144
 145 //_______________________________________________________________________
 146 AliAnalysisTaskEmcalJetFlavourTagExample::~AliAnalysisTaskEmcalJetFlavourTagExample()
 147 {
 148   // destructor
 149   //
 150   if (fOutput) {
 151     delete fOutput;
 152     fOutput = 0;
 153   }
 154 }
 155
 156 //________________________________________________________________________
 157 void AliAnalysisTaskEmcalJetFlavourTagExample::UserCreateOutputObjects()
 158 {
 159   if (! fCreateHisto)
 160     return;
 161   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 162
 163   //fJetsCont           = GetJetContainer(0);
 164   if(fJetsCont) { //get particles and clusters connected to jets
 165     fTracksJetCont       = fJetsCont->GetParticleContainer();
 166     fCaloClustersJetCont = fJetsCont->GetClusterContainer();
 167   }
 168  else {        //no jets, just analysis tracks and clusters
 169   fTracksCont       = GetParticleContainer(0);
 170   fCaloClustersCont = GetClusterContainer(0);
 171 }
 172 fTracksCont->SetClassName("AliVTrack");
 173 fCaloClustersCont->SetClassName("AliVCluster");
 174
 175   fHistJetPhi                = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
 176   fHistJetPt                 = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
 177   fOutput->Add(fHistJetPhi);
 178   fOutput->Add(fHistJetPt);
 179
 180
 181   TString histname;
 182
 183   fHistCorJetPt                     = new TH1F("CorrJetPt", "CorrJetPt", 300, -100, 200);
 184   fHistnSigElecPt             = new TH2F("nsigvsPt(TPC)","nSigmaElectronTPC_v_Pt",60,0,60,40,-10,10);
 185   fHistdEdXvsPt               = new TH2F("dEdXvsTrackPt","dEdXvsTrackPt",60,0,60,80,0,80);
 186   fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
 187   fHistHighJetPt              = new TH1F("HighestPtJetPerEvent","HighJetPt",80,0,80);
 188
 189   TString name;
 190   TString title;
 191
 192   fOutput->Add(fHistCorJetPt);
 193   fOutput->Add(fHistnSigElecPt);
 194   fOutput->Add(fHistdEdXvsPt);
 195   fOutput->Add(fHistnJetTrackvnJetClusters);
 196   fOutput->Add(fHistHighJetPt);
 197
 198
 199   // ****************************** PID *****************************************************
 200   // set up PID handler
 201   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 202   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 203   if(!inputHandler) {
 204     AliFatal("Input handler needed");
 205   }
 206
 207   // PID response object
 208   //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
 209   //  inputHandler->CreatePIDResponse(fIsMC);         // needed to create object, why though?
 210   fPIDResponse = inputHandler->GetPIDResponse();
 211   if (!fPIDResponse) {
 212     AliError("PIDResponse object was not created");
 213   }
 214   // ***************************************************************************************
 215
 216   PostData(1, fOutput);
 217
 218 }
 219
 220 //________________________________________________________
 221 void AliAnalysisTaskEmcalJetFlavourTagExample::ExecOnce()
 222 {
 223   //  Initialize the analysis
 224   AliAnalysisTaskEmcalJet::ExecOnce();
 225
 226   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 227   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 228   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 229
 230
 231 } // end of ExecOnce
 232
 233 //________________________________________________________________________
 234 Bool_t AliAnalysisTaskEmcalJetFlavourTagExample::Run()
 235 {
 236   // check to see if we have any tracks
 237   if (!fTracks)  return kTRUE;
 238   if (!fJets)  return kTRUE;
 239
 240   // what kind of event do we have: AOD or ESD?
 241   Bool_t useAOD;
 242   if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
 243   else useAOD = kFALSE;
 244
 245   // if we have ESD event, set up ESD object
 246   if(!useAOD){
 247     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 248     if (!fESD) {
 249       AliError(Form("ERROR: fESD not available\n"));
 250       return kTRUE;
 251     }
 252   }
 253
 254   // if we have AOD event, set up AOD object
 255   if(useAOD){
 256     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 257     if(!fAOD) {
 258       AliError(Form("ERROR: fAOD not available\n"));
 259       return kTRUE;
 260     }
 261   }
 262
 263   // get magnetic field info for DCA
 264   Double_t  MagF = fESD->GetMagneticField();
 265   Double_t MagSign = 1.0;
 266   if(MagF<0)MagSign = -1.0;
 267   // set magnetic field
 268   if (!TGeoGlobalMagField::Instance()->GetField()) {
 269     AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
 270     TGeoGlobalMagField::Instance()->SetField(field);
 271   }
 272
 273   // get vertex information
 274   Double_t fvertex[3]={0,0,0};
 275   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 276   //Double_t zVtx=fvertex[2];
 277
 278   // create pointer to list of input event
 279   TList *list = InputEvent()->GetList();
 280   if(!list) {
 281     AliError(Form("ERROR: list not attached\n"));
 282     return kTRUE;
 283   }
 284
 285   // background density
 286   fRhoVal = fRho->GetVal();
 287
 288   // initialize TClonesArray pointers to jets and tracks
 289   TClonesArray *jets = 0;
 290   //TClonesArray *tracks = 0;
 291   //TClonesArray *clusters = 0;
 292
 293   // get Jets object
 294   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 295   if(!jets){
 296     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
 297     return kTRUE;
 298   } // verify existence of jets
 299
 300   // get number of jets and tracks
 301   const Int_t Njets = jets->GetEntries();
 302   if(Njets<1)     return kTRUE;
 303
 304   if (fTracksCont) {
 305     AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 306     while(track) {
 307       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 308     }
 309   }
 310   if (fCaloClustersCont) {
 311     AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 312     while(cluster) {
 313       TLorentzVector nPart;
 314       cluster->GetMomentum(nPart, fVertex);
 315       cluster = fCaloClustersCont->GetNextAcceptCluster();
 316     }
 317   }
 318
 319     //  Start Jet Analysis
 320     // initialize jet parameters
 321     Int_t ijethi=-1;
 322     Double_t highestjetpt=0.0;
 323
 324   // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
 325   for (Int_t ijet = 0; ijet < Njets; ijet++){
 326     // get our jets
 327     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 328     if (!jet) continue;
 329
 330     // apply jet cuts
 331     if(!AcceptMyJet(jet)) continue;
 332
 333     if(highestjetpt<jet->Pt()){
 334       ijethi=ijet;
 335       highestjetpt=jet->Pt();
 336     }
 337   } // end of looping over jets
 338
 339     fHistHighJetPt->Fill(ijethi);
 340
 341  // **********************************************************************
 342  //                JET LOOP
 343  // **********************************************************************
 344
 345     // loop over jets in the event and make appropriate cuts
 346     for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 347       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 348       if (!jet)  // see if we have a jet
 349         continue;
 350
 351       // apply jet cuts
 352       if(!AcceptMyJet(jet)) continue;
 353
 354       //AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
 355       //jet->AddFlavourTag(tag);
 356
 357
 358       Int_t JetClusters = jet->GetNumberOfClusters();
 359       Int_t JetTracks = jet -> GetNumberOfTracks();
 360       fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
 361       // Initializations and Calculations
 362       Double_t jetPt = -500;                                     // initialize corr jet pt LOCAL
 363       Double_t jetarea = -500;                                   // initialize jet area
 364       jetPt = jet->Pt() - jetarea*fRhoVal;                       // semi-corrected pT of jet from GLOBAL rho value
 365       fHistCorJetPt->Fill(jetPt);
 366
 367       Bool_t bkgrnd1 = kFALSE;
 368
 369       if(jet->Pt() > fJetHIpt) {
 370         if(!fTracksCont || !fCaloClustersCont) continue;
 371
 372         //******************************Cluster Matched To Closest Track
 373         //**************************************************************
 374         Int_t NumbTrackContainer = -999;
 375         NumbTrackContainer = fTracksCont->GetNParticles();
 376         for(int iTracks = 0; iTracks <= NumbTrackContainer; iTracks++){
 377           AliVTrack *AcceptedTrack =static_cast<AliVTrack*>(fTracksCont->GetParticle(iTracks));
 378           if(!AcceptedTrack){
 379             AliError(Form("Couldn't get AliVTrack Container %d\n", iTracks));
 380             continue;
 381           }
 382           if(!IsJetTrack(jet,iTracks,kFALSE))continue;
 383           //Get matched cluster
 384           Int_t emc1 = AcceptedTrack->GetEMCALcluster();
 385
 386           Double_t acceptTrackP = AcceptedTrack->P();
 387           Double_t acceptTrackPt = AcceptedTrack->Pt();
 388           Double_t nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(AcceptedTrack,AliPID::kElectron);
 389           fHistnSigElecPt->Fill(acceptTrackPt,nSigmaElectron_TPC);
 390
 391           AliESDtrack *ESDacceptedTrack = static_cast<AliESDtrack*>(AcceptedTrack);
 392           if(!ESDacceptedTrack){
 393             AliError(Form("Couldn't get AliESDTrack %d\n", iTracks));
 394             continue;
 395           }
 396           Double_t dEdx = AcceptedTrack->GetTPCsignal();
 397                               fHistdEdXvsPt->Fill(acceptTrackPt,dEdx);
 398           if(fCaloClustersCont && emc1>=0) {
 399             AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
 400             if(!clusMatch){
 401               AliError(Form("Couldn't get matched AliVCluster %d\n", emc1));
 402               continue;
 403             }
 404
 405             Double_t mClusterE = clusMatch->E();
 406             Double_t EovP_mc = -999;
 407             EovP_mc = mClusterE/acceptTrackP;
 408
 409             if(EovP_mc < 0.2) bkgrnd1 = kTRUE;
 410           }
 411
 412       } //loop over tracks for matching to closest cluster
 413
 414
 415
 416
 417     } // highest pt jet cut
 418
 419
 420       if(bkgrnd1 == kTRUE) {
 421         AliEmcalJet::EFlavourTag tag=AliEmcalJet::kBckgrd1;
 422         jet->AddFlavourTag(tag);
 423     }
 424
 425   } // LOOP over JETS in event
 426
 427
 428   return kTRUE;
 429
 430 }
 431 //________________________________________________________________________
 432 void AliAnalysisTaskEmcalJetFlavourTagExample::Terminate(Option_t *)
 433 {
 434   cout<<"###########################"<<endl;
 435   cout<<"####   Task Finished   ####"<<endl;
 436   cout<<"###########################"<<endl;
 437   cout<<"###########################"<<endl;
 438 } // end of terminate
 439
 440 //________________________________________________________________________
 441 Int_t AliAnalysisTaskEmcalJetFlavourTagExample::AcceptMyJet(AliEmcalJet *jet) {
 442   //applies all jet cuts except pt
 443   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
 444   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
 445   if (jet->Area()<fAreacut) return 0;
 446   // prevents 0 area jets from sneaking by when area cut == 0
 447   if (jet->Area()==0) return 0;
 448   //exclude jets with extremely high pt tracks which are likely misreconstructed
 449   if(jet->MaxTrackPt()>100) return 0;
 450
 451   //passed all above cuts
 452   return 1;
 453 }
 454
 455
 456
 457
 458