PWG4/JetTasks/AliAnalysisTaskJetServices.cxx

   1 // **************************************
   2 // Task used for the correction of determiantion of reconstructed jet spectra
   3 // Compares input (gen) and output (rec) jets
   4 // *******************************************
   5
   6
   7 /**************************************************************************
   8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   9  *                                                                        *
  10  * Author: The ALICE Off-line Project.                                    *
  11  * Contributors are mentioned in the code where appropriate.              *
  12  *                                                                        *
  13  * Permission to use, copy, modify and distribute this software and its   *
  14  * documentation strictly for non-commercial purposes is hereby granted   *
  15  * without fee, provided that the above copyright notice appears in all   *
  16  * copies and that both the copyright notice and this permission notice   *
  17  * appear in the supporting documentation. The authors make no claims     *
  18  * about the suitability of this software for any purpose. It is          *
  19  * provided "as is" without express or implied warranty.                  *
  20  **************************************************************************/
  21
  22
  23 #include <TROOT.h>
  24 #include <TRandom.h>
  25 #include <TString.h>
  26 #include <TSystem.h>
  27 #include <TInterpreter.h>
  28 #include <TChain.h>
  29 #include <TFile.h>
  30 #include <TKey.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TH3F.h>
  34 #include <TProfile.h>
  35 #include <TList.h>
  36 #include <TLorentzVector.h>
  37 #include <TClonesArray.h>
  38 #include  "TDatabasePDG.h"
  39
  40 #include "AliAnalysisTaskJetServices.h"
  41 #include "AliAnalysisDataContainer.h"
  42 #include "AliAnalysisDataSlot.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliJetFinder.h"
  45 #include "AliJetHeader.h"
  46 #include "AliJetReader.h"
  47 #include "AliJetReaderHeader.h"
  48 #include "AliUA1JetHeaderV1.h"
  49 #include "AliESDEvent.h"
  50 #include "AliAODEvent.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODTrack.h"
  53 #include "AliAODJet.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliMCEvent.h"
  57 #include "AliStack.h"
  58 #include "AliGenPythiaEventHeader.h"
  59 #include "AliJetKineReaderHeader.h"
  60 #include "AliGenCocktailEventHeader.h"
  61 #include "AliInputEventHandler.h"
  62 #include "AliPhysicsSelection.h"
  63
  64
  65 #include "AliAnalysisHelperJetTasks.h"
  66
  67 ClassImp(AliAnalysisTaskJetServices)
  68
  69 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
  70   fUseAODInput(kFALSE),
  71   fUsePhysicsSelection(kFALSE),
  72   fRealData(kFALSE),
  73   fAvgTrials(1),
  74   fZVtxCut(8.),
  75   fh1Xsec(0x0),
  76   fh1Trials(0x0),
  77   fh1PtHard(0x0),
  78   fh1PtHardTrials(0x0),
  79   fh2TriggerCount(0x0),
  80   fh2ESDTriggerCount(0x0),
  81   fh2TriggerVtx(0x0),
  82   fh2ESDTriggerVtx(0x0),
  83   fh2ESDTriggerRun(0x0),
  84   fh2VtxXY(0x0),
  85   fHistList(0x0)
  86 {
  87   fRunRange[0] = fRunRange[1] = 0;
  88 }
  89
  90 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
  91   AliAnalysisTaskSE(name),
  92   fUseAODInput(kFALSE),
  93   fUsePhysicsSelection(kFALSE),
  94   fRealData(kFALSE),
  95   fAvgTrials(1),
  96   fZVtxCut(8.),
  97   fh1Xsec(0x0),
  98   fh1Trials(0x0),
  99   fh1PtHard(0x0),
 100   fh1PtHardTrials(0x0),
 101   fh2TriggerCount(0x0),
 102   fh2ESDTriggerCount(0x0),
 103   fh2TriggerVtx(0x0),
 104   fh2ESDTriggerVtx(0x0),
 105   fh2ESDTriggerRun(0x0),
 106   fh2VtxXY(0x0),
 107   fHistList(0x0)
 108 {
 109   fRunRange[0] = fRunRange[1] = 0;
 110   DefineOutput(1,TList::Class());
 111 }
 112
 113
 114
 115 Bool_t AliAnalysisTaskJetServices::Notify()
 116 {
 117   //
 118   // Implemented Notify() to read the cross sections
 119   // and number of trials from pyxsec.root
 120   //
 121
 122   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 123   Float_t xsection = 0;
 124   Float_t ftrials  = 1;
 125
 126   fAvgTrials = 1;
 127   if(tree){
 128     TFile *curfile = tree->GetCurrentFile();
 129     if (!curfile) {
 130       Error("Notify","No current file");
 131       return kFALSE;
 132     }
 133     if(!fh1Xsec||!fh1Trials){
 134       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 135       return kFALSE;
 136     }
 137     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 138     fh1Xsec->Fill("<#sigma>",xsection);
 139     // construct a poor man average trials
 140     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 141     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 142   }
 143   return kTRUE;
 144 }
 145
 146 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 147 {
 148
 149   //
 150   // Create the output container
 151   //
 152
 153
 154   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 155
 156   OpenFile(1);
 157   if(!fHistList)fHistList = new TList();
 158
 159   Bool_t oldStatus = TH1::AddDirectoryStatus();
 160   TH1::AddDirectory(kFALSE);
 161   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 162   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 163   fHistList->Add(fh1Xsec);
 164
 165   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 166   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 167   fHistList->Add(fh1Trials);
 168
 169   const Int_t nBinPt = 100;
 170   Double_t binLimitsPt[nBinPt+1];
 171   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 172     if(iPt == 0){
 173       binLimitsPt[iPt] = 0.0;
 174     }
 175     else {// 1.0
 176       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 177     }
 178   }
 179
 180   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 181   fHistList->Add(fh2TriggerCount);
 182
 183   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 184   fHistList->Add(fh2ESDTriggerCount);
 185
 186   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 187   fHistList->Add(fh2TriggerVtx);
 188
 189   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 190   fHistList->Add(fh2ESDTriggerVtx);
 191
 192
 193   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 194   fHistList->Add(fh1PtHard);
 195   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 196   fHistList->Add(fh1PtHardTrials);
 197
 198   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 199   // 3 triggers BB BE/EB EE
 200
 201   fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
 202   fHistList->Add(fh2ESDTriggerRun);
 203
 204   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 205   fHistList->Add(fh2VtxXY);
 206   // =========== Switch on Sumw2 for all histos ===========
 207   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 208     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 209     if (h1){
 210       h1->Sumw2();
 211       continue;
 212     }
 213     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 214     if(hn)hn->Sumw2();
 215   }
 216
 217
 218   TH1::AddDirectory(oldStatus);
 219 }
 220
 221 void AliAnalysisTaskJetServices::Init()
 222 {
 223   //
 224   // Initialization
 225   //
 226   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 227
 228 }
 229
 230 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 231 {
 232
 233   //
 234   // Execute analysis for current event
 235   //
 236
 237   AliAODEvent *aod = 0;
 238   AliESDEvent *esd = 0;
 239
 240   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 241
 242   if(fUseAODInput){
 243     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 244     if(!aod){
 245       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 246       return;
 247     }
 248     // fethc the header
 249   }
 250   else{
 251     //  assume that the AOD is in the general output...
 252     aod  = AODEvent();
 253     if(!aod){
 254       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 255       return;
 256     }
 257     esd = dynamic_cast<AliESDEvent*>(InputEvent());
 258   }
 259
 260   if(esd){
 261     Float_t run = (Float_t)esd->GetRunNumber();
 262     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 263     Float_t zvtx = -999;
 264     Float_t xvtx = -999;
 265     Float_t yvtx = -999;
 266
 267     if(vtxESD->GetNContributors()>0){
 268       zvtx = vtxESD->GetZ();
 269       yvtx = vtxESD->GetY();
 270       xvtx = vtxESD->GetX();
 271     }
 272
 273     Int_t iTrig = -1;
 274     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
 275        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
 276        ||esd->GetFiredTriggerClasses().Contains("MB1")
 277        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
 278       iTrig = 0;
 279     }
 280     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
 281             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
 282             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
 283             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
 284             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
 285             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
 286       // empty bunch or bunch empty
 287       iTrig = 1;
 288     }
 289     if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
 290        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
 291       iTrig = 2;
 292     }
 293
 294
 295     if(iTrig>=0){
 296       iTrig *= 3;
 297       fh2ESDTriggerRun->Fill(run,iTrig+1);
 298       if(vtxESD->GetNContributors()>0){
 299         fh2ESDTriggerRun->Fill(run,iTrig+2);
 300         fh2VtxXY->Fill(xvtx,yvtx);
 301       }
 302       Float_t r2 = xvtx *xvtx + yvtx *yvtx;
 303       if(TMath::Abs(zvtx)<fZVtxCut&&r2<1)fh2ESDTriggerRun->Fill(run,iTrig+3);
 304     }
 305     else{
 306       fh2ESDTriggerRun->Fill(run,0);
 307     }
 308
 309
 310   }
 311
 312   // loop over all possible trigger and
 313   for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 314     Bool_t esdTrig = kFALSE;
 315     if(esd){
 316       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
 317       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
 318     }
 319     Bool_t aodTrig = kFALSE;
 320     if(aod){
 321       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 322       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 323     }
 324
 325     // Check wether we have also an SPD vertex
 326
 327     if(aod){
 328       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 329       //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
 330
 331       if(vtxAOD->GetNContributors()>0){
 332         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
 333         Float_t zvtx = vtxAOD->GetZ();
 334         Float_t yvtx = vtxAOD->GetY();
 335         Float_t xvtx = vtxAOD->GetX();
 336         fh2TriggerVtx->Fill(it,zvtx);
 337         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 338         if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&r2<1){
 339           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 340         }
 341       }
 342     }
 343     if(esd){
 344       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 345       //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
 346       Bool_t cand = fInputHandler->IsEventSelected();
 347       if(cand){
 348         fh2ESDTriggerCount->Fill(it,kSelectedALICE);
 349       }
 350       if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
 351
 352       if(vtxESD->GetNContributors()>0){
 353         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
 354         Float_t zvtx = vtxESD->GetZ();
 355         Float_t yvtx = vtxESD->GetY();
 356         Float_t xvtx = vtxESD->GetX();
 357         Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 358         fh2ESDTriggerVtx->Fill(it,zvtx);
 359         Bool_t vtxIn = TMath::Abs(zvtx)<fZVtxCut&&r2<1;
 360
 361         if(cand&&vtxIn){
 362           fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
 363           fh2ESDTriggerCount->Fill(it,kSelected);
 364           AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 365         }
 366         if(vtxIn&&esdTrig){
 367           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
 368           // here we select based on ESD info...
 369         }
 370       }
 371     }
 372   }
 373
 374
 375
 376   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 377
 378
 379   Double_t ptHard = 0;
 380   Double_t nTrials = 1; // Trials for MC trigger
 381
 382   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 383   AliMCEvent* mcEvent = MCEvent();
 384   //    AliStack *pStack = 0;
 385   if(mcEvent){
 386     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 387     if(pythiaGenHeader){
 388       nTrials = pythiaGenHeader->Trials();
 389       ptHard  = pythiaGenHeader->GetPtHard();
 390       int iProcessType = pythiaGenHeader->ProcessType();
 391       // 11 f+f -> f+f
 392       // 12 f+barf -> f+barf
 393       // 13 f+barf -> g+g
 394       // 28 f+g -> f+g
 395       // 53 g+g -> f+barf
 396       // 68 g+g -> g+g
 397       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 398       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 399       fh1PtHard->Fill(ptHard);
 400       fh1PtHardTrials->Fill(ptHard,nTrials);
 401
 402     }// if pythia gen header
 403   }
 404
 405   // trigger selection
 406
 407
 408   PostData(1, fHistList);
 409 }
 410
 411
 412 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 413 {
 414   // Terminate analysis
 415   //
 416 }