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