PWGLF/SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskJetSpectraAOD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 //         AliAnalysisTaskJetSpectraAOD class
  18 //-----------------------------------------------------------------
  19
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TLegend.h"
  23 #include "TH1F.h"
  24 #include "TH2F.h"
  25 #include "THnSparse.h"
  26 #include "TCanvas.h"
  27 #include "AliAnalysisTask.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliVTrack.h"
  30 #include "AliAODMCParticle.h"
  31 #include "AliVParticle.h"
  32 #include "AliAODEvent.h"
  33 #include "AliAODTrack.h"
  34 #include "AliAODInputHandler.h"
  35 #include "AliAnalysisTaskJetSpectraAOD.h"
  36 #include "AliAnalysisTaskESDfilter.h"
  37 #include "AliAnalysisDataContainer.h"
  38 #include "AliCentrality.h"
  39 #include "TProof.h"
  40 #include "AliVEvent.h"
  41 #include "AliStack.h"
  42 #include <TMCProcess.h>
  43
  44 #include "AliSpectraAODTrackCuts.h"
  45 #include "AliSpectraAODEventCuts.h"
  46
  47 //jet include
  48 #include "AliAODHandler.h"
  49 #include "AliAODJetEventBackground.h"
  50 #include "AliAODJet.h"
  51 #include "AliAnalysisHelperJetTasks.h"
  52
  53 #include <iostream>
  54
  55 using namespace std;
  56
  57 ClassImp(AliAnalysisTaskJetSpectraAOD)
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskJetSpectraAOD::AliAnalysisTaskJetSpectraAOD(const char *name) : AliAnalysisTaskSE(name),
  61   fAOD(0),
  62   fIsMC(0),
  63   fEventCuts(0),
  64   fTrackCuts(0),
  65   fVZEROside(0),
  66   fOutput(0),
  67   fAODJets(0),
  68   fJetBranchName(""),
  69   fListJets(0),
  70   fBackgroundBranch(""),
  71   fOfflineTrgMask(AliVEvent::kMB),
  72   fFilterMask(0),
  73   fJetPtMin(0x0),
  74   fJetEtaMin(0x0),
  75   fJetEtaMax(0x0),
  76   fLeadPtMin(0x0),
  77   fnCentBins(20),
  78   fnQvecBins(20),
  79   fnptLeadBins(4),
  80   fIsQvecCalibMode(0),
  81   fQvecUpperLim(100),
  82   fIsQvecCut(0),
  83   fQvecMin(0),
  84   fQvecMax(100),
  85   fHistEvtSelection(0x0),
  86   fDebug(0),
  87   fMinNcontributors(0),
  88   fRejectPileup(0),
  89   fR(0.4),
  90   fZvertexDiff(1),
  91   fZvertex(10.)
  92 {
  93   // Default constructor
  94
  95   DefineInput(0, TChain::Class());
  96   DefineOutput(1, TList::Class());
  97   DefineOutput(2, AliSpectraAODEventCuts::Class());
  98   DefineOutput(3, AliSpectraAODTrackCuts::Class());
  99
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
 104 {
 105    delete fListJets;
 106 }
 107 //________________________________________________________________________
 108 //________________________________________________________________________
 109 void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
 110 {
 111   // create output objects
 112   fOutput = new TList();
 113   fOutput->SetOwner();
 114   fOutput->SetName("chistpt");
 115
 116   fListJets = new TList;
 117
 118   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
 119   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
 120
 121   // binning common to all the THn
 122   const Double_t ptBins[] = {-50.,-45.,-40.,-35.,-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,65.,70.,75.,80.,85.,90.,100.,150.,200.};
 123   const Int_t nptBins=31;
 124
 125   //dimensions of THnSparse for jets
 126   const Int_t nvarjet=5;
 127   //                                        pt_raw    pt_corr           cent             Q vec                pt_lead
 128   Int_t    binsHistRealJet[nvarjet] = {    nptBins,   nptBins,       fnCentBins,      fnQvecBins,       fnptLeadBins};
 129   Double_t xminHistRealJet[nvarjet] = {         0.,        0.,             0.,                0.,                  0.};
 130   Double_t xmaxHistRealJet[nvarjet] = {       200.,      200.,           100.,     fQvecUpperLim,                 20.};
 131   THnSparseF* NSparseHistJet = new THnSparseF("NSparseHistJet","NSparseHistJet",nvarjet,binsHistRealJet,xminHistRealJet,xmaxHistRealJet);
 132   NSparseHistJet->GetAxis(0)->SetTitle("#it{p}_{T,raw}");
 133   NSparseHistJet->GetAxis(0)->SetName("pT_raw");
 134   NSparseHistJet->SetBinEdges(0,ptBins);
 135   NSparseHistJet->GetAxis(1)->SetTitle("#it{p}_{T,corr}");
 136   NSparseHistJet->GetAxis(1)->SetName("pT_corr");
 137   NSparseHistJet->SetBinEdges(1,ptBins);
 138   NSparseHistJet->GetAxis(2)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 139   NSparseHistJet->GetAxis(2)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 140   NSparseHistJet->GetAxis(3)->SetTitle("Q vec");
 141   NSparseHistJet->GetAxis(3)->SetName("Q_vec");
 142   NSparseHistJet->GetAxis(4)->SetTitle("#it{p}_{T,lead}");
 143   NSparseHistJet->GetAxis(4)->SetName("pT_lead");
 144   fOutput->Add(NSparseHistJet);
 145
 146   //dimensions of THnSparse for the normalization
 147   const Int_t nvarev=3;
 148   //                                             cent         Q vec         rho
 149   Int_t    binsHistRealEv[nvarev] = {    fnCentBins,      fnQvecBins,      100 };
 150   Double_t xminHistRealEv[nvarev] = {           0.,               0.,        0.};
 151   Double_t xmaxHistRealEv[nvarev] = {         100.,      fQvecUpperLim,    200.};
 152   THnSparseF* NSparseHistEv = new THnSparseF("NSparseHistEv","NSparseHistEv",nvarev,binsHistRealEv,xminHistRealEv,xmaxHistRealEv);
 153   NSparseHistEv->GetAxis(0)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 154   NSparseHistEv->GetAxis(0)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 155   NSparseHistEv->GetAxis(1)->SetTitle("Q vec");
 156   NSparseHistEv->GetAxis(1)->SetName("Q_vec");
 157   NSparseHistEv->GetAxis(2)->SetTitle("rho");
 158   NSparseHistEv->GetAxis(2)->SetName("rho");
 159   fOutput->Add(NSparseHistEv);
 160
 161 //   TH1F* fHistTest = new TH1F("fHistTest","fHistTest",nptBins-1,ptBins);
 162 //   fOutput->Add(fHistTest);
 163
 164   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 7.5);
 165   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 166   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 167   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 168   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"centrality (rejected)");
 169   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"multiplicity (rejected)");
 170   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"ESE (rejected)");
 171   fOutput->Add(fHistEvtSelection);
 172
 173   PostData(1, fOutput  );
 174   PostData(2, fEventCuts);
 175   PostData(3, fTrackCuts);
 176
 177 }
 178
 179 //________________________________________________________________________
 180 void AliAnalysisTaskJetSpectraAOD::UserExec(Option_t *)
 181 {
 182
 183   // check for jet branches
 184   if(!strlen(fJetBranchName.Data())){
 185     AliError("Jet branch name not set.");
 186     return;
 187   }
 188
 189   // main event loop
 190   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);
 191   if (!fAOD) {
 192     AliWarning("ERROR: AliAODEvent not available \n");
 193     return;
 194   }
 195
 196   if (strcmp(fAOD->ClassName(), "AliAODEvent"))
 197     {
 198       AliFatal("Not processing AODs");
 199     }
 200
 201   TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 202   if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
 203     fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
 204   }
 205
 206   /* -- event selection -- */
 207   fHistEvtSelection->Fill(1); // number of events before event selection
 208
 209   //jet service task event selection.
 210   Bool_t selected=kTRUE;
 211   selected = AliAnalysisHelperJetTasks::Selected();
 212   if(!selected){
 213     // no selection by the service task, we continue
 214     PostData(1,fOutput);
 215     PostData(2, fEventCuts);
 216     PostData(3, fTrackCuts);
 217   return;}
 218
 219   // physics selection: this is now redundant, all should appear as accepted after service task selection
 220   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 221   ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 222   //std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
 223   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 224     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 225     fHistEvtSelection->Fill(2);
 226     PostData(1,fOutput);
 227     PostData(2, fEventCuts);
 228     PostData(3, fTrackCuts);
 229     return;
 230   }
 231
 232   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 233   Int_t nTracksPrim = primVtx->GetNContributors();
 234
 235   if (fDebug) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
 236   if(nTracksPrim < fMinNcontributors){
 237     if (fDebug) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
 238     fHistEvtSelection->Fill(3);
 239     PostData(1,fOutput);
 240     PostData(2, fEventCuts);
 241     PostData(3, fTrackCuts);
 242     return;
 243   }
 244
 245   TString primVtxName(primVtx->GetName());
 246
 247   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
 248     if (fDebug) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
 249     fHistEvtSelection->Fill(4);
 250     PostData(1,fOutput);
 251     PostData(2, fEventCuts);
 252     PostData(3, fTrackCuts);
 253     return;
 254   }
 255
 256   if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
 257     if (fDebug) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
 258     fHistEvtSelection->Fill(5);
 259     PostData(1,fOutput);
 260     PostData(2, fEventCuts);
 261     PostData(3, fTrackCuts);
 262     return;
 263   }
 264
 265   //cut on difference of the z-vertex
 266   if (fZvertexDiff || (fZvertex>0)) {
 267
 268     Double_t vzPRI = +999;
 269     Double_t vzSPD = -999;
 270
 271     const AliVVertex *pv = fAOD->GetPrimaryVertex(); // AOD
 272     if (pv && pv->GetNContributors()>0) {
 273       vzPRI = pv->GetZ();
 274     }
 275
 276     const AliVVertex *sv = fAOD->GetPrimaryVertexSPD();
 277     if (sv && sv->GetNContributors()>0) {
 278       vzSPD = sv->GetZ();
 279     }
 280
 281     Double_t  dvertex = TMath::Abs(vzPRI-vzSPD);
 282
 283     if (fZvertexDiff && (dvertex>0.1)) {
 284       if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
 285       return;}
 286     if ((fZvertex>0) && (TMath::Abs(vzPRI)>fZvertex)) {
 287       if (fDebug) Printf("%s:%d ZvertexDiff Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
 288       return;}
 289   }
 290
 291   //ESE event cuts
 292   if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
 293     if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
 294     fHistEvtSelection->Fill(6);
 295     PostData(1,fOutput);
 296     PostData(2, fEventCuts);
 297     PostData(3, fTrackCuts);
 298     return;
 299   }
 300
 301   if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
 302   fHistEvtSelection->Fill(0.);
 303   // accepted events
 304   // -- end event selection --
 305
 306
 307   Double_t Qvec=0.;//in case of MC we save space in the memory
 308   if(!fIsMC){
 309     if(fIsQvecCalibMode){
 310       if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
 311       else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
 312     }
 313     else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
 314   }
 315
 316   if(fIsQvecCut && (Qvec<fQvecMin || Qvec>fQvecMax) ) return;
 317
 318   Double_t Cent=fEventCuts->GetCent();
 319
 320   // get background
 321   AliAODJetEventBackground* externalBackground = 0;
 322   if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
 323     externalBackground =  (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
 324     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 325   }
 326
 327   Float_t rho = 0;
 328   if(externalBackground)rho = externalBackground->GetBackground(0);  //default schema
 329     if(rho==0) rho=-9999.; //rho value = 0 are non-physical -> removed from the distribution
 330
 331   // fetch jets
 332   TClonesArray *aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fJetBranchName.Data()));
 333   if(!aodJets){
 334     AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
 335     if(fAOD){
 336       Printf("Input AOD >>>>");
 337       fAOD->Print();
 338     }
 339     return;
 340   }
 341
 342   Bool_t kDeltaPt = kFALSE;
 343   if(fJetBranchName.Contains("RandomCone")){
 344 //     cout<<"RANDOM CONES BRANCH!"<<endl;
 345     kDeltaPt = kTRUE;
 346   }
 347
 348
 349   fListJets->Clear();
 350   for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
 351     AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
 352     if (jet) fListJets->Add(jet);
 353   }
 354
 355
 356
 357   for(Int_t i=0; i<fListJets->GetEntries(); ++i){
 358     AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
 359
 360     if((jet->Eta()<fJetEtaMin)||(jet->Eta()>fJetEtaMax)) continue;
 361
 362     Double_t ptJet   = jet->Pt();
 363
 364     Double_t areaJet = jet->EffectiveAreaCharged();
 365
 366     Double_t ptcorr = -999.;
 367     if(kDeltaPt) ptcorr = ptJet - (rho*TMath::Pi()*fR*fR);
 368     else ptcorr = ptJet-(rho*areaJet);
 369
 370     Double_t varJet[5];
 371     varJet[0]=jet->Pt();
 372     varJet[1]=(Double_t)ptcorr;
 373     varJet[2]=(Double_t)Cent;
 374     varJet[3]=(Double_t)Qvec;
 375
 376     if(!kDeltaPt){
 377       AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
 378       if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
 379
 380       varJet[4]=(Double_t)leadTrack->Pt();
 381     }
 382     else varJet[4] = -1;
 383
 384     ((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
 385   }
 386
 387
 388 //   //track loop
 389 //   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {  //FIXME loop on aod track... should be on jet tracks???
 390 //     AliAODTrack* track = fAOD->GetTrack(iTracks);
 391 // //     if(!(track->TestFilterBit(1024)))continue;
 392 //     if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
 393 //
 394 //     TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
 395 //
 396 //   } // end loop on tracks
 397
 398   Double_t varEv[3];
 399   varEv[0]=Cent;
 400   varEv[1]=Qvec;
 401   varEv[2]=rho;
 402   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
 403
 404   PostData(1,fOutput);
 405   PostData(2, fEventCuts);
 406   PostData(3, fTrackCuts);
 407   //Printf("............. end of Exec");
 408
 409 }
 410
 411 //_________________________________________________________________
 412 AliVParticle *AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet* jet){
 413   if(!jet)return 0;
 414   TRefArray *refs = jet->GetRefTracks();
 415   if(!refs) return 0;
 416   AliVParticle *leading = 0;
 417   Float_t fMaxPt = 0;
 418   for(int i = 0;i<refs->GetEntriesFast();i++){
 419     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
 420     if(!tmp)continue;
 421     if(tmp->Pt()>fMaxPt){
 422       leading = tmp;
 423       fMaxPt = tmp->Pt();
 424     }
 425   }
 426   return leading;
 427 }
 428
 429 //_________________________________________________________________
 430 void   AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
 431 {
 432   // Terminate
 433 }
 434
 435 //jet
 436 void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
 437 {
 438    fJetBranchName = branch;
 439 }
 440
 441 void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
 442 {
 443    fBackgroundBranch = bckbranch;
 444 }