PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetPatchTriggerQA.cxx

   1 // ChristineQA task
   2 //
   3 // Author: J Mazer
   4
   5 #include "AliAnalysisTaskEmcalJetPatchTriggerQA.h"
   6
   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
  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 #include "AliRhoParameter.h"
  30 #include "AliEmcalParticle.h"
  31 #include "AliLocalRhoParameter.h"
  32 #include "AliAnalysisTaskLocalRho.h"
  33 #include <AliInputEventHandler.h>
  34 #include <AliVEventHandler.h>
  35
  36 ClassImp(AliAnalysisTaskEmcalJetPatchTriggerQA)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskEmcalJetPatchTriggerQA::AliAnalysisTaskEmcalJetPatchTriggerQA() :
  40   AliAnalysisTaskEmcalJet("ChristineQA",kFALSE),
  41   fPhimin(-10), fPhimax(10),
  42   fEtamin(-0.9), fEtamax(0.9),
  43   fAreacut(0.0),
  44   fLocalRhoVal(0),
  45   fHistNjetvsCent(0),
  46   fhnJetTriggerQA(0x0)
  47 {
  48   // Default constructor.
  49   SetMakeGeneralHistograms(kTRUE);
  50
  51 }
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskEmcalJetPatchTriggerQA::AliAnalysisTaskEmcalJetPatchTriggerQA(const char *name) :
  55   AliAnalysisTaskEmcalJet(name,kTRUE),
  56   fPhimin(-10), fPhimax(10),
  57   fEtamin(-0.9), fEtamax(0.9),
  58   fAreacut(0.0),
  59   fLocalRhoVal(0),
  60   fHistNjetvsCent(0),
  61   fhnJetTriggerQA(0x0)
  62 {
  63
  64    SetMakeGeneralHistograms(kTRUE);
  65
  66    // DefineInput(0,TChain::Class());
  67    // DefineOutput(1, TList::Class());
  68 }
  69
  70 //_______________________________________________________________________
  71 AliAnalysisTaskEmcalJetPatchTriggerQA::~AliAnalysisTaskEmcalJetPatchTriggerQA()
  72 {
  73   // destructor
  74   //
  75   if (fOutput) {
  76     delete fOutput;
  77     fOutput = 0;
  78   }
  79 }
  80
  81 //________________________________________________________________________
  82 void AliAnalysisTaskEmcalJetPatchTriggerQA::UserCreateOutputObjects()
  83 {
  84   if (! fCreateHisto)
  85     return;
  86   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  87
  88   fHistNjetvsCent            = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
  89   fOutput->Add(fHistNjetvsCent);
  90
  91   UInt_t bitcode = 0; // bit codes, see GetDimParams() below
  92   //bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9;
  93   bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 ;
  94   fhnJetTriggerQA = NewTHnSparseF("fhnJetTriggerQA", bitcode);
  95   fhnJetTriggerQA->Sumw2();
  96   fOutput->Add(fhnJetTriggerQA);
  97
  98   PostData(1, fOutput);
  99 }
 100
 101 //________________________________________________________
 102 void AliAnalysisTaskEmcalJetPatchTriggerQA::ExecOnce()
 103 {
 104   //  Initialize the analysis
 105   AliAnalysisTaskEmcalJet::ExecOnce();
 106
 107 } // end of ExecOnce
 108
 109 //________________________________________________________________________
 110 Bool_t AliAnalysisTaskEmcalJetPatchTriggerQA::Run()
 111 {
 112   // TEST TEST TEST TEST for OBJECTS!!
 113   if(!fLocalRho) {
 114     AliWarning(Form("%s: Could not retrieve LocalRho object!", GetName()));
 115     fLocalRho = GetLocalRhoFromEvent(fLocalRhoName);
 116   }
 117
 118   // check to see if we have jet object
 119   if (!fJets)  return kTRUE;
 120
 121   // find NUMBER of jets
 122   const Int_t Njets = fJets->GetEntries();
 123   Int_t NjetAcc = 0;
 124
 125   // loop over jets in the event and make appropriate cuts
 126   //cout<<"I at least get in the event and I have "<<Njets<<" jets"<<endl;
 127   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 128      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 129      if (!jet)  // see if we have a jet
 130        continue;
 131      if (!AcceptMyJet(jet))
 132        continue;
 133
 134      //cout<<"jet pt "<<jet->Pt()<<" area "<<jet->Area()<<" maxtrackpt "<<jet->MaxTrackPt()<<endl;
 135      //This somehow needs to be fixed but I'm not sure what it does yet.  It seems the defaults are wacky.
 136 //     if (! AcceptJet(jet)) // sees if jet is accepted
 137 //     continue;
 138      //  jets.push_back(jet);
 139
 140      NjetAcc++;
 141
 142      // Initializations and Calculations
 143      // Double_t jetphi = jet->Phi();
 144      //Double_t jeteta = jet->Eta();                                    // ETA of jet
 145
 146      Double_t jetPtraw = jet->Pt();                             // raw pT of jet
 147      Double_t jetPtLocal = -500;                    // initialize corr jet pt LOCAL
 148      Double_t jetPtGlobal = -500;                                       // initialize corr jet pt GLOBAL
 149      Double_t jetarea = -500;                                   // initialize jet area
 150      jetarea = jet->Area();                                     // jet area
 151
 152      Double_t dEP = -500;                                       // initialize angle between jet and event plane
 153      dEP = RelativeEPJET(jet->Phi(),fEPV0);
 154
 155      // get LOCAL rho from event and fill histo's
 156      fRhoVal = fRho->GetVal();
 157      fLocalRhoVal = fLocalRho->GetLocalVal(jet->Phi(), 0.2);  // jet angle, then jet radius
 158
 159      jetPtLocal = jet->Pt() - jetarea*fLocalRhoVal;              // corrected pT of jet from LOCAL rho value
 160      jetPtGlobal = jet->Pt() - jetarea*fRhoVal;                  // corrected pT of jet from GLOBAL rho value
 161
 162      // need to update entries soon
 163      //Double_t Entries[10] = {centbin, jetPtLocal, jetPtGlobal, jetPtraw, jet->Eta(), jet->Phi(), dEP, fLocalRhoVal, fRhoVal, fEPV0};
 164      Double_t Entries[9] = {fCent, jetPtLocal, jetPtGlobal, jetPtraw, jet->Phi(), dEP, fLocalRhoVal, fRhoVal, fEPV0};
 165      fhnJetTriggerQA->Fill(Entries);        // fill Sparse Histo with entries
 166
 167      // in plane and out of plane histo's
 168      if( dEP>0 && dEP<=(TMath::Pi()/6) ){
 169         // we are IN plane, lets fill some histo's
 170      }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
 171         // we are OUT of PLANE, lets fill some histo's
 172          }
 173
 174     fHistNjetvsCent->Fill(fCent,NjetAcc);
 175
 176   } // LOOP over JETS in event
 177
 178   return kTRUE;
 179 }
 180
 181 //________________________________________________________________________
 182 void AliAnalysisTaskEmcalJetPatchTriggerQA::Terminate(Option_t *)
 183 {
 184
 185 } // end of terminate
 186
 187 //________________________________________________________________________
 188 Int_t AliAnalysisTaskEmcalJetPatchTriggerQA::GetCentBin(Double_t cent) const
 189 {  // Get centrality bin.
 190
 191   Int_t centbin = -1;
 192   if (cent>=0 && cent<10)
 193     centbin = 0;
 194   else if (cent>=10 && cent<20)
 195     centbin = 1;
 196   else if (cent>=20 && cent<30)
 197     centbin = 2;
 198   else if (cent>=30 && cent<40)
 199     centbin = 3;
 200   else if (cent>=40 && cent<50)
 201     centbin = 4;
 202   else if (cent>=50 && cent<90)
 203     centbin = 5;
 204   return centbin;
 205 }
 206
 207 //_________________________________________________________________________
 208 Double_t AliAnalysisTaskEmcalJetPatchTriggerQA:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
 209 { // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
 210   Double_t dphi = (EPAng - jetAng);
 211
 212   // ran into trouble with a few dEP<-Pi so trying this...
 213   if( dphi<-1*TMath::Pi() ){
 214     dphi = dphi + 1*TMath::Pi();
 215   }
 216
 217   if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
 218     // Do nothing! we are in quadrant 1
 219   }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
 220     dphi = 1*TMath::Pi() - dphi;
 221   }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
 222     dphi = fabs(dphi);
 223   }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
 224     dphi = dphi + 1*TMath::Pi();
 225   }
 226
 227   // test
 228   if( dphi < 0 || dphi > TMath::Pi()/2 )
 229     AliWarning(Form("%s: dPHI not constrained from [0, Pi/2]", GetName()));
 230
 231   return dphi;   // dphi in [0, Pi/2]
 232 }
 233
 234 //______________________________________________________________________
 235 THnSparse* AliAnalysisTaskEmcalJetPatchTriggerQA::NewTHnSparseF(const char* name, UInt_t entries)
 236 {
 237    // generate new THnSparseF, axes are defined in GetDimParams()
 238    Int_t count = 0;
 239    UInt_t tmp = entries;
 240    while(tmp!=0){
 241       count++;
 242       tmp = tmp &~ -tmp;  // clear lowest bit
 243    }
 244
 245    TString hnTitle(name);
 246    const Int_t dim = count;
 247    Int_t nbins[dim];
 248    Double_t xmin[dim];
 249    Double_t xmax[dim];
 250
 251    Int_t i=0;
 252    Int_t c=0;
 253    while(c<dim && i<32){
 254       if(entries&(1<<i)){
 255
 256          TString label("");
 257          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
 258          hnTitle += Form(";%s",label.Data());
 259          c++;
 260       }
 261
 262       i++;
 263    }
 264    hnTitle += ";";
 265
 266    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 267 } // end of NewTHnSparseF
 268
 269 void AliAnalysisTaskEmcalJetPatchTriggerQA::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 270 {
 271    // stores label and binning of axis for THnSparse
 272    const Double_t pi = TMath::Pi();
 273
 274    switch(iEntry){
 275
 276    case 0:
 277       label = "V0 centrality (%)";
 278       nbins = 20;
 279       xmin = 0.;
 280       xmax = 100.;
 281       break;
 282
 283    case 1:
 284       label = "Corrected jet p_{T}: Local #rho";
 285       nbins = 500;
 286       xmin = -250.;
 287       xmax = 250.;
 288       break;
 289
 290    case 2:
 291       label = "Corrected jet p_{T}: Global #rho";
 292       nbins = 500;
 293       xmin = -250.;
 294       xmax = 250.;
 295       break;
 296
 297    case 3:
 298       label = "Raw jet p_{T}";
 299       nbins = 250;
 300       xmin = 0;
 301       xmax = 250.;
 302       break;
 303
 304    case 4:
 305       label = "Jet Phi";
 306       nbins = 144;
 307       xmin = 1.4;//-1.0*pi;
 308       xmax = 3.4;//2.0*pi;
 309       break;
 310
 311   case 5:
 312       label = "Relative angle between jet and Reaction plane";
 313       nbins = 72;
 314       xmin = 0;
 315       xmax = 0.5*pi;
 316       break;
 317
 318   case 6:
 319       label = "Local #rho value";
 320       nbins = 120;
 321       xmin = 0.0;
 322       xmax = 300.0;
 323       break;
 324
 325   case 7:
 326       label = "Global #rho value";
 327       nbins = 120;
 328       xmin = 0.0;
 329       xmax = 300.0;
 330       break;
 331
 332   case 8:
 333       label = "fEPV0 event plane";
 334       nbins = 72;
 335       xmin = -TMath::Pi()/2.0;
 336       xmax = TMath::Pi()/2.0;
 337       break;
 338
 339    } // end of switch
 340 } // end of getting dim-params
 341
 342 //________________________________________________________________________
 343 Int_t AliAnalysisTaskEmcalJetPatchTriggerQA::AcceptMyJet(AliEmcalJet *jet) {
 344   //applies all jet cuts except pt
 345   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
 346   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
 347   if (jet->Area()<fAreacut) return 0;
 348   // prevents 0 area jets from sneaking by when area cut == 0
 349   if (jet->Area()==0) return 0;
 350   //exclude jets with extremely high pt tracks which are likely misreconstructed
 351   if(jet->MaxTrackPt()>100) return 0;
 352
 353   //passed all above cuts
 354   return 1;
 355 }