HMPID/AliAnalysisTaskJetsHMPID.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // AliAnalysysTaskJetsHMPID - Class performing PID analysis in jets with the HMPID
  18 // A set of histograms is created.
  19 //=================================================================================
  20 //
  21 // By means of AnalysisTrainHMPID.C macro it is possible to use this class
  22 // to perform the analysis on local data, on data on alien using local machine
  23 // and on CAF.
  24
  25 #include <TClonesArray.h>
  26 #include <TRandom2.h>
  27 #include "AliAnalysisTaskJetsHMPID.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliAODHandler.h"
  30 #include "AliAODInputHandler.h"
  31 #include "AliAODEvent.h"
  32 #include "AliAODJet.h"
  33 #include "AliAODJetEventBackground.h"
  34
  35 ClassImp(AliAnalysisTaskJetsHMPID)
  36
  37 //__________________________________________________________________________
  38 AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID() :
  39   fJetBranch("jets"),
  40   fBkgBranch(""),
  41   fJetPtCut(10.),
  42   fAOD(0x0),
  43   fHistList(0x0),
  44   fThetaChJet(0x0),
  45   fThetaChBkg(0x0),
  46   fThetaChRndCone(0x0),
  47   fEvSelDPhi(0x0),
  48   fJetsPt(0x0),
  49   fRndConePt(0x0),
  50   fAwayJetPt(0x0),
  51   fJetsEtaPhi(0x0),
  52   fTrksEtaPhiJet(0x0),
  53   fTrksEtaPhiBkg(0x0),
  54   fTree(0x0),
  55   fTrackPt(0),
  56   fJetPt(0),
  57   fPionBkg(1.1),
  58   fKaonBkg(1.1),
  59   fProtBkg(1.1),
  60   fPionJet(1.1),
  61   fKaonJet(1.1),
  62   fProtJet(1.1)
  63 {
  64   //
  65   //Default ctor
  66   //
  67 }
  68
  69 //__________________________________________________________________________
  70 AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID(const Char_t* name) :
  71   AliAnalysisTaskSE(name),
  72   fJetBranch("jets"),
  73   fBkgBranch(""),
  74   fJetPtCut(10.),
  75   fAOD(0x0),
  76   fHistList(0x0),
  77   fThetaChJet(0x0),
  78   fThetaChBkg(0x0),
  79   fThetaChRndCone(0x0),
  80   fEvSelDPhi(0x0),
  81   fJetsPt(0x0),
  82   fRndConePt(0x0),
  83   fAwayJetPt(0x0),
  84   fJetsEtaPhi(0x0),
  85   fTrksEtaPhiJet(0x0),
  86   fTrksEtaPhiBkg(0x0),
  87   fTree(0x0),
  88   fTrackPt(0),
  89   fJetPt(0),
  90   fPionBkg(1.1),
  91   fKaonBkg(1.1),
  92   fProtBkg(1.1),
  93   fPionJet(1.1),
  94   fKaonJet(1.1),
  95   fProtJet(1.1)
  96 {
  97   //
  98   // Constructor. Initialization of Inputs and Outputs
  99   //
 100
 101   DefineOutput(1,TList::Class());
 102   DefineOutput(2,TTree::Class());
 103 }
 104
 105 //___________________________________________________________________________
 106 AliAnalysisTaskJetsHMPID& AliAnalysisTaskJetsHMPID::operator=(const AliAnalysisTaskJetsHMPID& c)
 107 {
 108   //
 109   // Assignment operator
 110   //
 111   if (this!=&c) {
 112     AliAnalysisTaskSE::operator=(c);
 113     fJetBranch       = c.fJetBranch;
 114     fBkgBranch       = c.fBkgBranch;
 115     fJetPtCut        = c.fJetPtCut;
 116     fAOD             = c.fAOD;
 117     fHistList        = c.fHistList;
 118     fThetaChJet      = c.fThetaChJet;
 119     fThetaChBkg      = c.fThetaChBkg;
 120     fThetaChRndCone  = c.fThetaChRndCone;
 121     fEvSelDPhi       = c.fEvSelDPhi;
 122     fJetsPt          = c.fJetsPt;
 123     fRndConePt       = c.fRndConePt;
 124     fAwayJetPt       = c.fAwayJetPt;
 125     fJetsEtaPhi      = c.fJetsEtaPhi;
 126     fTrksEtaPhiJet   = c.fTrksEtaPhiJet;
 127     fTrksEtaPhiBkg   = c.fTrksEtaPhiBkg;
 128     fTree            = c.fTree;
 129     fTrackPt         = c.fTrackPt;
 130     fJetPt           = c.fJetPt;
 131     fPionBkg         = c.fPionBkg;
 132     fKaonBkg         = c.fKaonBkg;
 133     fProtBkg         = c.fProtBkg;
 134     fPionJet         = c.fPionJet;
 135     fKaonJet         = c.fKaonJet;
 136     fProtJet         = c.fProtJet;
 137   }
 138   return *this;
 139 }
 140
 141 //___________________________________________________________________________
 142 AliAnalysisTaskJetsHMPID::AliAnalysisTaskJetsHMPID(const AliAnalysisTaskJetsHMPID& c) :
 143   AliAnalysisTaskSE(c),
 144   fJetBranch(c.fJetBranch),
 145   fBkgBranch(c.fBkgBranch),
 146   fJetPtCut(c.fJetPtCut),
 147   fAOD(c.fAOD),
 148   fHistList(c.fHistList),
 149   fThetaChJet(c.fThetaChJet),
 150   fThetaChBkg(c.fThetaChBkg),
 151   fThetaChRndCone(c.fThetaChRndCone),
 152   fEvSelDPhi(c.fEvSelDPhi),
 153   fJetsPt(c.fJetsPt),
 154   fRndConePt(c.fRndConePt),
 155   fAwayJetPt(c.fAwayJetPt),
 156   fJetsEtaPhi(c.fJetsEtaPhi),
 157   fTrksEtaPhiJet(c.fTrksEtaPhiJet),
 158   fTrksEtaPhiBkg(c.fTrksEtaPhiBkg),
 159   fTree(c.fTree),
 160   fTrackPt(c.fTrackPt),
 161   fJetPt(c.fJetPt),
 162   fPionBkg(c.fPionBkg),
 163   fKaonBkg(c.fKaonBkg),
 164   fProtBkg(c.fProtBkg),
 165   fPionJet(c.fPionJet),
 166   fKaonJet(c.fKaonJet),
 167   fProtJet(c.fProtJet)
 168 {
 169   //
 170   // Copy Constructor
 171   //
 172 }
 173
 174 //___________________________________________________________________________
 175 AliAnalysisTaskJetsHMPID::~AliAnalysisTaskJetsHMPID()
 176 {
 177   //
 178   //destructor
 179   //
 180   Info("~AliAnalysisTaskJetsHMPID","Calling Destructor");
 181   if (fHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHistList;
 182 }
 183
 184 //__________________________________________________________________________
 185 void AliAnalysisTaskJetsHMPID::ConnectInputData(Option_t *)
 186 {
 187   if (fDebug) AliInfo("ConnectInputData() ");
 188
 189   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 190
 191   if(handler && handler->InheritsFrom("AliAODInputHandler"))
 192   { // input AOD
 193     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
 194     if (fDebug) AliInfo("Tracks and Jets from AliAODInputHandler");
 195   }
 196   else
 197   {  //output AOD
 198     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 199     if (handler && handler->InheritsFrom("AliAODHandler"))
 200     {
 201       fAOD = ((AliAODHandler*)handler)->GetAOD();
 202       if (fDebug > 1) AliInfo("Tracks and Jets from AliAODHandler");
 203     }
 204     else
 205     {  // no AOD
 206       AliWarning("No AOD event in the input or in the output");
 207     }
 208   }
 209
 210 }
 211
 212 //___________________________________________________________________________
 213 void AliAnalysisTaskJetsHMPID::UserCreateOutputObjects()
 214 {
 215   if (!fHistList) fHistList = new TList();
 216   fHistList->SetOwner();
 217
 218   fThetaChJet = new TH2F("ThetaChJet","Theta Cherenkov distribution in the jets",1000,0.,10.,400,0.,0.8);
 219   fHistList->Add(fThetaChJet);
 220
 221   fThetaChBkg = new TH2F("ThetaChBkg","Theta Cherenkov distribution out of the jets",1000,0.,10.,400,0.,0.8);
 222   fHistList->Add(fThetaChBkg);
 223
 224   fThetaChRndCone = new TH2F("ThetaChRndCone","Theta Cherenkov distribution in random cones",1000,0.,10.,400,0.,0.8);
 225   fHistList->Add(fThetaChRndCone);
 226
 227   fEvSelDPhi = new TH1F("PhiJets","Jets phi in different selections",600,0.,6.);
 228   fHistList->Add(fEvSelDPhi);
 229
 230   fJetsPt = new TH1F("JetsPt","Jets spectrum in the HMPID",200,0.,200.);
 231   fHistList->Add(fJetsPt);
 232
 233   fRndConePt = new TH1F("RndConePt","Pt of random cones in the HMPID",200,0.,200.);
 234   fHistList->Add(fRndConePt);
 235
 236   fAwayJetPt = new TH1F("AwayJetPt","Pt of jets opposite to the HMPID",200,0.,200.);
 237   fHistList->Add(fAwayJetPt);
 238
 239   fJetsEtaPhi = new TH2F("JetsEtaPhi","Eta-Phi of jets in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
 240   fHistList->Add(fJetsEtaPhi);
 241
 242   fTrksEtaPhiJet = new TH2F("TrksEtaPhiJet","Eta-Phi of jet tracks in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
 243   fHistList->Add(fTrksEtaPhiJet);
 244
 245   fTrksEtaPhiBkg = new TH2F("TrksEtaPhiBkg","Eta-Phi of bkg tracks in the HMPID",180,-0.9,0.9,200,-0.25*TMath::Pi(),0.75*TMath::Pi());
 246   fHistList->Add(fTrksEtaPhiBkg);
 247
 248   fTree = new TTree("Tree","Tree with data");
 249   fTree->Branch("Pt",&fTrackPt);
 250   fTree->Branch("PtJet",&fJetPt);
 251   fTree->Branch("PionBkg",&fPionBkg);
 252   fTree->Branch("KaonBkg",&fKaonBkg);
 253   fTree->Branch("ProtBkg",&fProtBkg);
 254   fTree->Branch("PionJet",&fPionJet);
 255   fTree->Branch("KaonJet",&fKaonJet);
 256   fTree->Branch("ProtJet",&fProtJet);
 257
 258   PostData(1,fHistList);
 259   PostData(2,fTree);
 260 }
 261
 262 //___________________________________________________________________________
 263 void AliAnalysisTaskJetsHMPID::UserExec(Option_t *)
 264 {
 265   if (!fAOD){
 266     PostData(1,fHistList);
 267     PostData(2,fTree);
 268     return;
 269   }
 270
 271   const Float_t pi = TMath::Pi(), jetR = 0.4;
 272   const Float_t hmpPhiMid = 30*TMath::DegToRad();
 273   const Float_t hmpPhiWid = 25*TMath::DegToRad();
 274
 275   TClonesArray* jets = (TClonesArray*) fAOD->FindListObject(fJetBranch.Data());
 276   if (fDebug) printf("Jet Branch: %s\n",fJetBranch.Data());
 277   if (fDebug) printf("Bkg Branch: %s\n",fBkgBranch.Data());
 278
 279   if (!jets) AliError("Jet branch not found");
 280
 281   Int_t nj = jets->GetEntriesFast();
 282   Float_t *jetPt = new Float_t[nj];
 283   Int_t *iPerm = new Int_t[nj];
 284 //  Correct jets Pt
 285   static AliAODJetEventBackground* externalBackground = 0;
 286   if(fBkgBranch.Length()) externalBackground = (AliAODJetEventBackground*)(fAOD->FindListObject(fBkgBranch.Data()));
 287   Float_t pTbkg = 0;
 288   AliAODJet *jet = 0;
 289   for(Int_t iJet=0; iJet<nj; iJet++){
 290     jet = (AliAODJet*)jets->At(iJet);
 291     jetPt[iJet] = jet->Pt();
 292     pTbkg = (externalBackground) ? externalBackground->GetBackground(2)*jet->EffectiveAreaCharged() : 0;
 293     jetPt[iJet] -= pTbkg;
 294   }
 295   TMath::Sort(nj,jetPt,iPerm,kTRUE);
 296
 297 //  Event characterisation:
 298 //  iEvSel = 0: no jets in the HMPID, consider it as background
 299 //  iEvSel = 1: jet in the HMPID, use it for PID in jets
 300 //  iEvSel = 2: HMPID is in the away side of a leading jet, but no jet found around
 301 //              the detector; dismiss the event as it might alter the analysis
 302
 303   Int_t iEvSel = -1, iJetSel = -1;
 304   Float_t phiTrig = 0;
 305   if (nj && jetPt[iPerm[0]] > fJetPtCut){
 306     phiTrig = ((AliAODJet*)jets->At(iPerm[0]))->Phi();
 307     if (phiTrig < pi+hmpPhiMid) phiTrig-=hmpPhiMid;
 308     else if (phiTrig > pi+hmpPhiMid) phiTrig-=(2*pi+hmpPhiMid);
 309     if (TMath::Abs(phiTrig) < hmpPhiWid+jetR) iEvSel = 1, iJetSel = 0;
 310     else if (TMath::Abs(phiTrig) > pi/3 && TMath::Abs(phiTrig) < 2*pi/3) iEvSel = 0;
 311     else if (TMath::Abs(phiTrig) > 2*pi/3) iEvSel = 2, iJetSel = 0;
 312   }
 313   else iEvSel = 0;
 314
 315   for (Int_t iJet=1; iJet<jets->GetEntriesFast(); iJet++){
 316     if (jetPt[iPerm[iJet]] < fJetPtCut || iEvSel == 1) break;
 317     jet = (AliAODJet*)jets->At(iPerm[iJet]);
 318     phiTrig = jet->Phi();
 319     if (phiTrig < pi+hmpPhiMid) phiTrig-=hmpPhiMid;
 320     else if (phiTrig > pi+hmpPhiMid) phiTrig-=(2*pi+hmpPhiMid);
 321     if (TMath::Abs(phiTrig) < hmpPhiWid+jetR) iEvSel = 1, iJetSel = iJet;
 322   }
 323
 324   TRandom2 rndm(0);
 325   Float_t etaRnd = rndm.Uniform(-0.8,0.8);
 326   Float_t phiRnd = rndm.Uniform(hmpPhiMid-hmpPhiWid-jetR,hmpPhiMid+hmpPhiWid+jetR);
 327   AliAODJet *rndJet = new AliAODJet();
 328   rndJet->SetPtEtaPhiM(10.,etaRnd,phiRnd,0.);
 329   rndJet->SetEffArea(pi*jetR*jetR,pi*jetR*jetR);
 330   Float_t ptSum=0., phi=0;
 331   TRefArray *jetTracks = 0;
 332   Double_t probs[5]={0, 0, 0, 0, 0};
 333   switch(iEvSel){
 334   case 0:
 335     for (Int_t iTr = 0; iTr < fAOD->GetNumberOfTracks(); iTr++){
 336       AliAODTrack *tr=(AliAODTrack*)fAOD->GetTrack(iTr);
 337       if (rndJet->DeltaR(tr) < jetR) ptSum+=tr->Pt();
 338       AliAODPid *pid = tr->GetDetPid();
 339       if (!pid) continue;
 340       if (pid->GetHMPIDsignal() > 0.){
 341         phi = tr->Phi()>pi ? tr->Phi()-2*pi : tr->Phi();
 342         fThetaChBkg->Fill(tr->Pt(),pid->GetHMPIDsignal());
 343         if (rndJet->DeltaR(tr) < jetR) fThetaChRndCone->Fill(tr->Pt(),pid->GetHMPIDsignal());
 344         fTrksEtaPhiBkg->Fill(tr->Eta(),phi);
 345         pid->GetHMPIDprobs(probs);
 346         Short_t charge = tr->Charge();
 347         fTrackPt = tr->Pt();
 348         fJetPt   = 0;
 349         fPionBkg = charge*(probs[0] + probs[1] + probs[2]);
 350         fKaonBkg = charge*probs[3];
 351         fProtBkg = charge*probs[4];
 352         fPionJet = charge*1.2;
 353         fKaonJet = charge*1.2;
 354         fProtJet = charge*1.2;
 355         fTree->Fill();
 356       }
 357     }
 358     pTbkg = (externalBackground) ? externalBackground->GetBackground(2)*jet->EffectiveAreaCharged() : 0;
 359     fRndConePt->Fill(ptSum-pTbkg);
 360     rndJet->Delete();
 361     break;
 362   case 1:
 363     jet = (AliAODJet*)jets->At(iPerm[iJetSel]);
 364     fJetsPt->Fill(jetPt[iPerm[iJetSel]]);
 365     phi = jet->Phi()>pi ? jet->Phi()-2*pi : jet->Phi();
 366     fJetsEtaPhi->Fill(jet->Eta(),phi);
 367     jetTracks = (TRefArray*)jet->GetRefTracks();
 368     for (Int_t iTr=0; iTr<jetTracks->GetEntriesFast(); iTr++){
 369       AliAODTrack *tr = (AliAODTrack*)jetTracks->At(iTr);
 370       AliAODPid *pid = tr->GetDetPid();
 371       if (!pid) continue;
 372       if (pid->GetHMPIDsignal() > 0.){
 373         phi = tr->Phi()>pi ? tr->Phi()-2*pi : tr->Phi();
 374         fThetaChJet->Fill(tr->Pt(),pid->GetHMPIDsignal());
 375         fTrksEtaPhiJet->Fill(tr->Eta(),phi);
 376         pid->GetHMPIDprobs(probs);
 377         Short_t charge = tr->Charge();
 378         fTrackPt = tr->Pt();
 379         fJetPt   = jetPt[iPerm[iJetSel]];
 380         fPionBkg = charge*1.2;
 381         fKaonBkg = charge*1.2;
 382         fProtBkg = charge*1.2;
 383         fPionJet = charge*(probs[0] + probs[1] + probs[2]);
 384         fKaonJet = charge*probs[3];
 385         fProtJet = charge*probs[4];
 386         fTree->Fill();
 387       }
 388     }
 389     break;
 390   case 2:
 391     fAwayJetPt->Fill(jetPt[iPerm[iJetSel]]);
 392     for (Int_t iTr = 0; iTr < fAOD->GetNumberOfTracks(); iTr++){
 393       AliAODTrack *tr=(AliAODTrack*)fAOD->GetTrack(iTr);
 394       AliAODPid *pid = tr->GetDetPid();
 395       if (!pid) continue;
 396       if (pid->GetHMPIDsignal() > 0.){
 397         fThetaChBkg->Fill(tr->Pt(),pid->GetHMPIDsignal());
 398         pid->GetHMPIDprobs(probs);
 399         Short_t charge = tr->Charge();
 400         fTrackPt = tr->Pt();
 401         fJetPt   = jetPt[iPerm[iJetSel]];
 402         fPionBkg = charge*(2 + probs[0] + probs[1] + probs[2]);
 403         fKaonBkg = charge*(2 + probs[3]);
 404         fProtBkg = charge*(2 + probs[4]);
 405         fPionJet = charge*1.2;
 406         fKaonJet = charge*1.2;
 407         fProtJet = charge*1.2;
 408         fTree->Fill();
 409       }
 410     }
 411     break;
 412   }
 413
 414   if (iEvSel > -1){
 415     if (iJetSel > -1){
 416       jet = (AliAODJet*)jets->At(iPerm[iJetSel]);
 417       fEvSelDPhi->Fill(jet->Phi()/pi+2*iEvSel);
 418     } else {
 419       for (Int_t iJet=0; iJet<jets->GetEntriesFast(); iJet++){
 420         if (jetPt[iPerm[iJet]] < fJetPtCut) break;
 421         jet = (AliAODJet*)jets->At(iPerm[iJet]);
 422         fEvSelDPhi->Fill(jet->Phi()/pi+2*iEvSel);
 423       }
 424     }
 425   }
 426
 427   delete [] jetPt;
 428   delete [] iPerm;
 429
 430   PostData(1,fHistList);
 431   PostData(2,fTree);
 432   return;
 433 }
 434
 435 //___________________________________________________________________________
 436 void AliAnalysisTaskJetsHMPID::Terminate(Option_t*)
 437 {
 438   AliAnalysisTaskSE::Terminate();
 439
 440 }