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 {
  90   // Default constructor
  91
  92   DefineInput(0, TChain::Class());
  93   DefineOutput(1, TList::Class());
  94   DefineOutput(2, AliSpectraAODEventCuts::Class());
  95   DefineOutput(3, AliSpectraAODTrackCuts::Class());
  96
  97 }
  98
  99 //________________________________________________________________________
 100 AliAnalysisTaskJetSpectraAOD::~AliAnalysisTaskJetSpectraAOD()
 101 {
 102    delete fListJets;
 103 }
 104 //________________________________________________________________________
 105 //________________________________________________________________________
 106 void AliAnalysisTaskJetSpectraAOD::UserCreateOutputObjects()
 107 {
 108   Printf("\n\n\n\n\n\n In CreateOutput Object:");
 109
 110   fOutput = new TList();
 111   fOutput->SetOwner();
 112   fOutput->SetName("chistpt");
 113
 114   fListJets = new TList;
 115
 116   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");
 117   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");
 118
 119   // binning common to all the THn
 120   const Double_t ptBins[] = {0.15,5.,10.,15.,20.,25.,30.,35.,40.,50.,75.,100.,150.,200.};
 121   const Int_t nptBins=13;
 122
 123   //dimensions of THnSparse for jets
 124   const Int_t nvarjet=6;
 125   //                                        pt_raw    pt_corr           cent             Q vec            rho          pt_lead
 126   Int_t    binsHistRealJet[nvarjet] = {    nptBins,   nptBins,       fnCentBins,      fnQvecBins,          40.,    fnptLeadBins};
 127   Double_t xminHistRealJet[nvarjet] = {         0.,        0.,             0.,                0.,           0.,              0.};
 128   Double_t xmaxHistRealJet[nvarjet] = {       200.,      200.,           100.,     fQvecUpperLim,         200.,             20.};
 129   THnSparseF* NSparseHistJet = new THnSparseF("NSparseHistJet","NSparseHistJet",nvarjet,binsHistRealJet,xminHistRealJet,xmaxHistRealJet);
 130   NSparseHistJet->GetAxis(0)->SetTitle("#it{p}_{T,raw}");
 131   NSparseHistJet->GetAxis(0)->SetName("pT_raw");
 132   NSparseHistJet->SetBinEdges(0,ptBins);
 133   NSparseHistJet->GetAxis(1)->SetTitle("#it{p}_{T,corr}");
 134   NSparseHistJet->GetAxis(1)->SetName("pT_corr");
 135   NSparseHistJet->SetBinEdges(1,ptBins);
 136   NSparseHistJet->GetAxis(2)->SetTitle(Form("%s cent",fEventCuts->GetCentralityMethod().Data()));
 137   NSparseHistJet->GetAxis(2)->SetName(Form("%s_cent",fEventCuts->GetCentralityMethod().Data()));
 138   NSparseHistJet->GetAxis(3)->SetTitle("Q vec");
 139   NSparseHistJet->GetAxis(3)->SetName("Q_vec");
 140   NSparseHistJet->GetAxis(4)->SetTitle("rho");
 141   NSparseHistJet->GetAxis(4)->SetName("rho");
 142   NSparseHistJet->GetAxis(5)->SetTitle("#it{p}_{T,lead}");
 143   NSparseHistJet->GetAxis(5)->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,       40.};
 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   //ESE event cuts
 266   if(!fEventCuts->IsSelected(fAOD,fTrackCuts)){
 267     if (fDebug) Printf("%s:%d ESE Event selection: event REJECTED...",(char*)__FILE__,__LINE__);
 268     fHistEvtSelection->Fill(6);
 269     PostData(1,fOutput);
 270     PostData(2, fEventCuts);
 271     PostData(3, fTrackCuts);
 272     return;
 273   }
 274
 275   if (fDebug) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
 276   fHistEvtSelection->Fill(0.);
 277   // accepted events
 278   // -- end event selection --
 279
 280
 281   Double_t Qvec=0.;//in case of MC we save space in the memory
 282   if(!fIsMC){
 283     if(fIsQvecCalibMode){
 284       if(fVZEROside==0)Qvec=fEventCuts->GetqV0A();
 285       else if (fVZEROside==1)Qvec=fEventCuts->GetqV0C();
 286     }
 287     else Qvec=fEventCuts->GetQvecPercentile(fVZEROside);
 288   }
 289
 290   if(fIsQvecCut && (Qvec<fQvecMin || Qvec>fQvecMax) ) return;
 291
 292   Double_t Cent=fEventCuts->GetCent();
 293
 294   // get background
 295   AliAODJetEventBackground* externalBackground = 0;
 296   if(fAODJets && !externalBackground && fBackgroundBranch.Length()){
 297     externalBackground =  (AliAODJetEventBackground*)(fAODJets->FindListObject(fBackgroundBranch.Data()));
 298     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 299   }
 300
 301   Float_t rho = 0;
 302   if(externalBackground)rho = externalBackground->GetBackground(0);  //default schema
 303
 304   // fetch jets
 305   TClonesArray *aodJets = dynamic_cast<TClonesArray*>(fAODJets->FindListObject(fJetBranchName.Data()));
 306   if(!aodJets){
 307     AliError(Form("no jet branch \"%s\" found, in the AODs are:", fJetBranchName.Data()));
 308     if(fAOD){
 309       Printf("Input AOD >>>>");
 310       fAOD->Print();
 311     }
 312     return;
 313   }
 314
 315   fListJets->Clear();
 316   for (Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
 317     AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
 318     if (jet) fListJets->Add(jet);
 319   }
 320
 321
 322
 323   for(Int_t i=0; i<fListJets->GetEntries(); ++i){
 324     AliAODJet* jet = (AliAODJet*)(fListJets->At(i));
 325
 326     if((jet->Eta()<fJetEtaMin)||(jet->Eta()>fJetEtaMax)) continue;
 327
 328     Double_t ptJet   = jet->Pt();
 329
 330     Double_t areaJet = jet->EffectiveAreaCharged();
 331
 332     Double_t ptcorr = ptJet-rho*areaJet;
 333
 334     Double_t varJet[6];
 335     varJet[0]=jet->Pt();
 336     varJet[1]=(Double_t)ptcorr;
 337     varJet[2]=(Double_t)Cent;
 338     varJet[3]=(Double_t)Qvec;
 339     varJet[4]=(Double_t)rho;
 340
 341     AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
 342     if(fLeadPtMin>0 && leadTrack->Pt()<fLeadPtMin)continue;
 343
 344     varJet[5]=(Double_t)leadTrack->Pt();
 345
 346     ((THnSparseF*)fOutput->FindObject("NSparseHistJet"))->Fill(varJet);//jet loop
 347   }
 348
 349
 350 //   //track loop
 351 //   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {  //FIXME loop on aod track... should be on jet tracks???
 352 //     AliAODTrack* track = fAOD->GetTrack(iTracks);
 353 // //     if(!(track->TestFilterBit(1024)))continue;
 354 //     if (!fTrackCuts->IsSelected(track,kTRUE)) continue;
 355 //
 356 //     TH1F* h=(TH1F*)fOutput->FindObject("fHistTest");h->Fill(track->Pt());
 357 //
 358 //   } // end loop on tracks
 359
 360   Double_t varEv[3];
 361   varEv[0]=Cent;
 362   varEv[1]=Qvec;
 363   varEv[2]=rho;
 364   ((THnSparseF*)fOutput->FindObject("NSparseHistEv"))->Fill(varEv);//event loop
 365
 366   PostData(1,fOutput);
 367   PostData(2, fEventCuts);
 368   PostData(3, fTrackCuts);
 369   //Printf("............. end of Exec");
 370
 371 }
 372
 373 //_________________________________________________________________
 374 AliVParticle *AliAnalysisTaskJetSpectraAOD::LeadingTrackFromJetRefs(AliAODJet* jet){
 375   if(!jet)return 0;
 376   TRefArray *refs = jet->GetRefTracks();
 377   if(!refs) return 0;
 378   AliVParticle *leading = 0;
 379   Float_t fMaxPt = 0;
 380   for(int i = 0;i<refs->GetEntriesFast();i++){
 381     AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
 382     if(!tmp)continue;
 383     if(tmp->Pt()>fMaxPt){
 384       leading = tmp;
 385       fMaxPt = tmp->Pt();
 386     }
 387   }
 388   return leading;
 389 }
 390
 391 //_________________________________________________________________
 392 void   AliAnalysisTaskJetSpectraAOD::Terminate(Option_t *)
 393 {
 394   // Terminate
 395   printf("AliAnalysisTaskJetSpectraAOD: Terminate() \n");
 396 }
 397
 398 //jet
 399 void AliAnalysisTaskJetSpectraAOD::SetBranchNames(const TString &branch)
 400 {
 401    fJetBranchName = branch;
 402 }
 403
 404 void AliAnalysisTaskJetSpectraAOD::SetRecBackgroundBranch(const TString &bckbranch)
 405 {
 406    fBackgroundBranch = bckbranch;
 407 }