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 <TSystem.h>
  26 #include <TInterpreter.h>
  27 #include <TChain.h>
  28 #include <TFile.h>
  29 #include <TKey.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TH3F.h>
  33 #include <TProfile.h>
  34 #include <TList.h>
  35 #include <TLorentzVector.h>
  36 #include <TClonesArray.h>
  37 #include  "TDatabasePDG.h"
  38
  39 #include "AliAnalysisTaskJetServices.h"
  40 #include "AliAnalysisDataContainer.h"
  41 #include "AliAnalysisDataSlot.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliJetFinder.h"
  44 #include "AliJetHeader.h"
  45 #include "AliJetReader.h"
  46 #include "AliJetReaderHeader.h"
  47 #include "AliUA1JetHeaderV1.h"
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 #include "AliAODHandler.h"
  51 #include "AliAODTrack.h"
  52 #include "AliAODJet.h"
  53 #include "AliAODMCParticle.h"
  54 #include "AliMCEventHandler.h"
  55 #include "AliMCEvent.h"
  56 #include "AliStack.h"
  57 #include "AliGenPythiaEventHeader.h"
  58 #include "AliJetKineReaderHeader.h"
  59 #include "AliGenCocktailEventHeader.h"
  60 #include "AliInputEventHandler.h"
  61 #include "AliPhysicsSelection.h"
  62
  63
  64 #include "AliAnalysisHelperJetTasks.h"
  65
  66 ClassImp(AliAnalysisTaskJetServices)
  67
  68 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
  69   fUseAODInput(kFALSE),
  70   fAvgTrials(1),
  71   fZVtxCut(8.),
  72   fRealData(kFALSE),
  73   fPhysicsSelection(0x0),
  74   fh1Xsec(0x0),
  75   fh1Trials(0x0),
  76   fh1PtHard(0x0),
  77   fh1PtHardTrials(0x0),
  78   fh2TriggerCount(0x0),
  79   fh2ESDTriggerCount(0x0),
  80   fh2TriggerVtx(0x0),
  81   fh2ESDTriggerVtx(0x0),
  82   fh2ESDTriggerRun(0x0),
  83   fh2VtxXY(0x0),
  84   fHistList(0x0)
  85 {
  86   fRunRange[0] = fRunRange[1] = 0;
  87 }
  88
  89 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
  90   AliAnalysisTaskSE(name),
  91   fUseAODInput(kFALSE),
  92   fAvgTrials(1),
  93   fZVtxCut(8.),
  94   fRealData(kFALSE),
  95   fPhysicsSelection(0x0),
  96   fh1Xsec(0x0),
  97   fh1Trials(0x0),
  98   fh1PtHard(0x0),
  99   fh1PtHardTrials(0x0),
 100   fh2TriggerCount(0x0),
 101   fh2ESDTriggerCount(0x0),
 102   fh2TriggerVtx(0x0),
 103   fh2ESDTriggerVtx(0x0),
 104   fh2ESDTriggerRun(0x0),
 105   fh2VtxXY(0x0),
 106   fHistList(0x0)
 107 {
 108   fRunRange[0] = fRunRange[1] = 0;
 109   DefineOutput(1,TList::Class());
 110 }
 111
 112
 113
 114 Bool_t AliAnalysisTaskJetServices::Notify()
 115 {
 116   //
 117   // Implemented Notify() to read the cross sections
 118   // and number of trials from pyxsec.root
 119   //
 120
 121   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 122   Float_t xsection = 0;
 123   Float_t ftrials  = 1;
 124
 125   fAvgTrials = 1;
 126   if(tree){
 127     TFile *curfile = tree->GetCurrentFile();
 128     if (!curfile) {
 129       Error("Notify","No current file");
 130       return kFALSE;
 131     }
 132     if(!fh1Xsec||!fh1Trials){
 133       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 134       return kFALSE;
 135     }
 136     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 137     fh1Xsec->Fill("<#sigma>",xsection);
 138     // construct a poor man average trials
 139     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 140     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 141   }
 142   return kTRUE;
 143 }
 144
 145 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 146 {
 147
 148   //
 149   // Create the output container
 150   //
 151
 152
 153   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 154
 155   OpenFile(1);
 156   if(!fHistList)fHistList = new TList();
 157
 158   Bool_t oldStatus = TH1::AddDirectoryStatus();
 159   TH1::AddDirectory(kFALSE);
 160
 161   if(!fPhysicsSelection)
 162     fPhysicsSelection = new AliPhysicsSelection();
 163   fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back
 164   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
 165   fHistList->Add(fPhysicsSelection);
 166
 167   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 168   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 169   fHistList->Add(fh1Xsec);
 170
 171   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 172   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 173   fHistList->Add(fh1Trials);
 174
 175   const Int_t nBinPt = 100;
 176   Double_t binLimitsPt[nBinPt+1];
 177   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 178     if(iPt == 0){
 179       binLimitsPt[iPt] = 0.0;
 180     }
 181     else {// 1.0
 182       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 183     }
 184   }
 185
 186   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 187   fHistList->Add(fh2TriggerCount);
 188
 189   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 190   fHistList->Add(fh2ESDTriggerCount);
 191
 192   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 193   fHistList->Add(fh2TriggerVtx);
 194
 195   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 196   fHistList->Add(fh2ESDTriggerVtx);
 197
 198
 199   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 200   fHistList->Add(fh1PtHard);
 201   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 202   fHistList->Add(fh1PtHardTrials);
 203
 204   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 205   // 3 triggers BB BE/EB EE
 206
 207   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);
 208   fHistList->Add(fh2ESDTriggerRun);
 209
 210   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 211   fHistList->Add(fh2VtxXY);
 212   // =========== Switch on Sumw2 for all histos ===========
 213   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 214     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 215     if (h1){
 216       h1->Sumw2();
 217       continue;
 218     }
 219     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 220     if(hn)hn->Sumw2();
 221   }
 222
 223
 224   TH1::AddDirectory(oldStatus);
 225 }
 226
 227 void AliAnalysisTaskJetServices::Init()
 228 {
 229   //
 230   // Initialization
 231   //
 232
 233   Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
 234   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 235
 236 }
 237
 238 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 239 {
 240
 241   //
 242   // Execute analysis for current event
 243   //
 244
 245   AliAODEvent *aod = 0;
 246   AliESDEvent *esd = 0;
 247
 248   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
 249
 250   if(fUseAODInput){
 251     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 252     if(!aod){
 253       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 254       return;
 255     }
 256     // fethc the header
 257   }
 258   else{
 259     //  assume that the AOD is in the general output...
 260     aod  = AODEvent();
 261     if(!aod){
 262       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 263       return;
 264     }
 265     esd = dynamic_cast<AliESDEvent*>(InputEvent());
 266   }
 267
 268   if(esd){
 269     Float_t run = (Float_t)esd->GetRunNumber();
 270     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 271     Float_t zvtx = -999;
 272     Float_t xvtx = -999;
 273     Float_t yvtx = -999;
 274
 275     if(vtxESD->GetNContributors()>0){
 276       zvtx = vtxESD->GetZ();
 277       yvtx = vtxESD->GetY();
 278       xvtx = vtxESD->GetX();
 279     }
 280
 281     Int_t iTrig = -1;
 282     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
 283        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
 284        ||esd->GetFiredTriggerClasses().Contains("MB1")
 285        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
 286       iTrig = 0;
 287     }
 288     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
 289             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
 290             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
 291             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
 292             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
 293             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
 294       // empty bunch or bunch empty
 295       iTrig = 1;
 296     }
 297     if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
 298        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
 299       iTrig = 2;
 300     }
 301
 302
 303     if(iTrig>=0){
 304       iTrig *= 3;
 305       fh2ESDTriggerRun->Fill(run,iTrig+1);
 306       if(vtxESD->GetNContributors()>0){
 307         fh2ESDTriggerRun->Fill(run,iTrig+2);
 308         fh2VtxXY->Fill(xvtx,yvtx);
 309       }
 310       if(TMath::Abs(zvtx)<fZVtxCut&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5)fh2ESDTriggerRun->Fill(run,iTrig+3);
 311     }
 312     else{
 313       fh2ESDTriggerRun->Fill(run,0);
 314     }
 315
 316
 317   }
 318
 319   // loop over all possible trigger and
 320   for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 321     Bool_t esdTrig = kFALSE;
 322     if(esd){
 323       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
 324       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
 325     }
 326     Bool_t aodTrig = kFALSE;
 327     if(aod){
 328       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 329       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 330     }
 331
 332     // Check wether we have also an SPD vertex
 333
 334     if(aod){
 335       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 336       //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
 337
 338       if(vtxAOD->GetNContributors()>0){
 339         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
 340         Float_t zvtx = vtxAOD->GetZ();
 341         Float_t yvtx = vtxAOD->GetY();
 342         Float_t xvtx = vtxAOD->GetX();
 343         fh2TriggerVtx->Fill(it,zvtx);
 344         if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
 345           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 346         }
 347       }
 348     }
 349     if(esd){
 350       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 351       //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
 352       Bool_t cand = fPhysicsSelection->IsCollisionCandidate(esd);
 353       if(cand) fh2ESDTriggerCount->Fill(it,kSelectedALICE);
 354       if(vtxESD->GetNContributors()>0){
 355         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
 356         Float_t zvtx = vtxESD->GetZ();
 357         Float_t yvtx = vtxESD->GetY();
 358         Float_t xvtx = vtxESD->GetX();
 359         fh2ESDTriggerVtx->Fill(it,zvtx);
 360         if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
 361           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
 362           // here we select based on ESD info...
 363           if(fRealData){
 364             if(cand){
 365               fh2ESDTriggerCount->Fill(it,kSelected);
 366               AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 367             }
 368           }
 369           else{
 370             fh2ESDTriggerCount->Fill(it,kSelected);
 371             AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 372           }
 373         }
 374
 375
 376       }
 377
 378     }
 379
 380   }
 381
 382
 383
 384   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 385
 386
 387   Double_t ptHard = 0;
 388   Double_t nTrials = 1; // Trials for MC trigger
 389
 390   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 391   AliMCEvent* mcEvent = MCEvent();
 392   //    AliStack *pStack = 0;
 393   if(mcEvent){
 394     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 395     if(pythiaGenHeader){
 396       nTrials = pythiaGenHeader->Trials();
 397       ptHard  = pythiaGenHeader->GetPtHard();
 398       int iProcessType = pythiaGenHeader->ProcessType();
 399       // 11 f+f -> f+f
 400       // 12 f+barf -> f+barf
 401       // 13 f+barf -> g+g
 402       // 28 f+g -> f+g
 403       // 53 g+g -> f+barf
 404       // 68 g+g -> g+g
 405       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 406       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 407       fh1PtHard->Fill(ptHard);
 408       fh1PtHardTrials->Fill(ptHard,nTrials);
 409
 410     }// if pythia gen header
 411   }
 412
 413   // trigger selection
 414
 415
 416   PostData(1, fHistList);
 417 }
 418
 419
 420 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 421 {
 422   // Terminate analysis
 423   //
 424
 425   TDirectory* owd = gDirectory;
 426
 427   if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
 428
 429   fHistList = dynamic_cast<TList*> (GetOutputData(1));
 430   if (!fHistList)
 431     Printf("ERROR: fOutput not available");
 432
 433
 434
 435   AliAnalysisDataContainer *cont = GetOutputSlot(1)->GetContainer();
 436   TString filename = cont->GetFileName();
 437   TFile *f = NULL;
 438   // Check first if the file is already opened
 439   f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
 440   if (f) {
 441     // Cd to file
 442     f->cd();
 443     // Check for a folder request
 444     TString dir = cont->GetFolderName();
 445     if (!dir.IsNull()) {
 446       if (!f->GetDirectory(dir)) f->mkdir(dir);
 447       f->cd(dir);
 448     }
 449   }
 450
 451   if (fHistList)
 452   {
 453     fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (fHistList->FindObject("AliPhysicsSelection_outputlist"));
 454   }
 455
 456   if (fPhysicsSelection)
 457     {
 458       fPhysicsSelection->SaveHistograms("physics_selection");
 459       fPhysicsSelection->Print();
 460     }
 461
 462   owd->cd();
 463
 464 }