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   fVtxXMean(0),
  75   fVtxYMean(0),
  76   fVtxZMean(0),
  77   fVtxRCut(1.),
  78   fVtxZCut(8.),
  79   fh1Xsec(0x0),
  80   fh1Trials(0x0),
  81   fh1PtHard(0x0),
  82   fh1PtHardTrials(0x0),
  83   fh2TriggerCount(0x0),
  84   fh2ESDTriggerCount(0x0),
  85   fh2TriggerVtx(0x0),
  86   fh2ESDTriggerVtx(0x0),
  87   fh2ESDTriggerRun(0x0),
  88   fh2VtxXY(0x0),
  89   fHistList(0x0)
  90 {
  91   fRunRange[0] = fRunRange[1] = 0;
  92 }
  93
  94 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
  95   AliAnalysisTaskSE(name),
  96   fUseAODInput(kFALSE),
  97   fUsePhysicsSelection(kFALSE),
  98   fRealData(kFALSE),
  99   fAvgTrials(1),
 100   fVtxXMean(0),
 101   fVtxYMean(0),
 102   fVtxZMean(0),
 103   fVtxRCut(1.),
 104   fVtxZCut(8.),
 105   fh1Xsec(0x0),
 106   fh1Trials(0x0),
 107   fh1PtHard(0x0),
 108   fh1PtHardTrials(0x0),
 109   fh2TriggerCount(0x0),
 110   fh2ESDTriggerCount(0x0),
 111   fh2TriggerVtx(0x0),
 112   fh2ESDTriggerVtx(0x0),
 113   fh2ESDTriggerRun(0x0),
 114   fh2VtxXY(0x0),
 115   fHistList(0x0)
 116 {
 117   fRunRange[0] = fRunRange[1] = 0;
 118   DefineOutput(1,TList::Class());
 119 }
 120
 121
 122
 123 Bool_t AliAnalysisTaskJetServices::Notify()
 124 {
 125   //
 126   // Implemented Notify() to read the cross sections
 127   // and number of trials from pyxsec.root
 128   //
 129
 130   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 131   Float_t xsection = 0;
 132   Float_t ftrials  = 1;
 133
 134   fAvgTrials = 1;
 135   if(tree){
 136     TFile *curfile = tree->GetCurrentFile();
 137     if (!curfile) {
 138       Error("Notify","No current file");
 139       return kFALSE;
 140     }
 141     if(!fh1Xsec||!fh1Trials){
 142       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 143       return kFALSE;
 144     }
 145     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 146     fh1Xsec->Fill("<#sigma>",xsection);
 147     // construct a poor man average trials
 148     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 149     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 150   }
 151   return kTRUE;
 152 }
 153
 154 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 155 {
 156
 157   //
 158   // Create the output container
 159   //
 160
 161
 162   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
 163
 164   OpenFile(1);
 165   if(!fHistList)fHistList = new TList();
 166
 167   Bool_t oldStatus = TH1::AddDirectoryStatus();
 168   TH1::AddDirectory(kFALSE);
 169   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 170   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 171   fHistList->Add(fh1Xsec);
 172
 173   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 174   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 175   fHistList->Add(fh1Trials);
 176
 177   const Int_t nBinPt = 100;
 178   Double_t binLimitsPt[nBinPt+1];
 179   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 180     if(iPt == 0){
 181       binLimitsPt[iPt] = 0.0;
 182     }
 183     else {// 1.0
 184       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
 185     }
 186   }
 187
 188   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 189   fHistList->Add(fh2TriggerCount);
 190
 191   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5);
 192   fHistList->Add(fh2ESDTriggerCount);
 193
 194   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 195   fHistList->Add(fh2TriggerVtx);
 196
 197   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.);
 198   fHistList->Add(fh2ESDTriggerVtx);
 199
 200
 201   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 202   fHistList->Add(fh1PtHard);
 203   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 204   fHistList->Add(fh1PtHardTrials);
 205
 206   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range
 207   // 3 triggers BB BE/EB EE
 208
 209   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);
 210   fHistList->Add(fh2ESDTriggerRun);
 211
 212   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
 213   fHistList->Add(fh2VtxXY);
 214   // =========== Switch on Sumw2 for all histos ===========
 215   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
 216     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
 217     if (h1){
 218       h1->Sumw2();
 219       continue;
 220     }
 221     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
 222     if(hn)hn->Sumw2();
 223   }
 224
 225
 226   TH1::AddDirectory(oldStatus);
 227 }
 228
 229 void AliAnalysisTaskJetServices::Init()
 230 {
 231   //
 232   // Initialization
 233   //
 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()>2){
 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()>2){
 307         fh2ESDTriggerRun->Fill(run,iTrig+2);
 308         fh2VtxXY->Fill(xvtx,yvtx);
 309       }
 310       xvtx -= fVtxXMean;
 311       yvtx -= fVtxYMean;
 312       zvtx -= fVtxZMean;
 313       Float_t r2 = xvtx *xvtx + yvtx *yvtx;
 314       if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
 315     }
 316     else{
 317       fh2ESDTriggerRun->Fill(run,0);
 318     }
 319
 320
 321   }
 322
 323
 324   // Apply additional constraints
 325   Bool_t esdEventSelected = IsEventSelectedESD(esd);
 326   Bool_t aodEventSelected = IsEventSelectedAOD(aod);
 327
 328
 329   if(esd&&aod&&fDebug){
 330     if(esdEventSelected&&!aodEventSelected){
 331       Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
 332       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 333       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 334       Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
 335       vtxESD->Print();
 336       Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
 337       vtxAOD->Print();
 338     }
 339   }
 340
 341   // loop over all possible triggers for esd
 342
 343   if(esd){
 344     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
 345     //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
 346     TString vtxName(vtxESD->GetName());
 347     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 348       Bool_t esdTrig = kFALSE;
 349       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
 350       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
 351       Bool_t cand = fInputHandler->IsEventSelected();
 352       if(cand){
 353         fh2ESDTriggerCount->Fill(it,kSelectedALICE);
 354       }
 355       if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
 356       if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
 357         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
 358         Float_t zvtx = vtxESD->GetZ();
 359         if(esdEventSelected&&esdTrig){
 360           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
 361           fh2ESDTriggerVtx->Fill(it,zvtx);
 362         }
 363       }
 364       if(cand&&esdEventSelected){
 365         fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
 366         fh2ESDTriggerCount->Fill(it,kSelected);
 367         AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
 368       }
 369     }
 370   }
 371
 372   if(aod){
 373     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 374     //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
 375     TString vtxTitle(vtxAOD->GetTitle());
 376     for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
 377       Bool_t aodTrig = kFALSE;
 378       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
 379       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
 380       if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
 381         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
 382         Float_t zvtx = vtxAOD->GetZ();
 383         if(aodTrig&&aodEventSelected){
 384           fh2TriggerVtx->Fill(it,zvtx);
 385           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
 386         }
 387       }
 388     }
 389   }
 390
 391   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 392
 393
 394   Double_t ptHard = 0;
 395   Double_t nTrials = 1; // Trials for MC trigger
 396
 397   fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 398   AliMCEvent* mcEvent = MCEvent();
 399   //    AliStack *pStack = 0;
 400   if(mcEvent){
 401     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 402     if(pythiaGenHeader){
 403       nTrials = pythiaGenHeader->Trials();
 404       ptHard  = pythiaGenHeader->GetPtHard();
 405       int iProcessType = pythiaGenHeader->ProcessType();
 406       // 11 f+f -> f+f
 407       // 12 f+barf -> f+barf
 408       // 13 f+barf -> g+g
 409       // 28 f+g -> f+g
 410       // 53 g+g -> f+barf
 411       // 68 g+g -> g+g
 412       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 413       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 414       fh1PtHard->Fill(ptHard);
 415       fh1PtHardTrials->Fill(ptHard,nTrials);
 416
 417     }// if pythia gen header
 418   }
 419
 420   // trigger selection
 421
 422
 423   PostData(1, fHistList);
 424 }
 425
 426 Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
 427   if(!esd)return kFALSE;
 428   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 429   // We can check the type of the vertex by its name and title
 430   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 431   // vtx       name                  title
 432   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 433   // TPC       TPCVertex             VertexerTracksNoConstraint
 434   // SPD       SPDVertex             vertexer: 3D
 435   // SPD       SPDVertex             vertexer: Z
 436
 437
 438
 439   if(vtx->GetNContributors()<3)return kFALSE;
 440
 441   // do not want tpc only primary vertex
 442   TString vtxName(vtx->GetName());
 443   if(vtxName.Contains("TPCVertex"))return kFALSE;
 444   Float_t zvtx = vtx->GetZ();
 445   Float_t yvtx = vtx->GetY();
 446   Float_t xvtx = vtx->GetX();
 447
 448   // here we may fill the histos for run dependence....
 449
 450   xvtx -= fVtxXMean;
 451   yvtx -= fVtxYMean;
 452   zvtx -= fVtxZMean;
 453
 454   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 455
 456   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 457   return eventSel;
 458 }
 459
 460 Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
 461   if(!aod)return kFALSE;
 462   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 463   // We can check the type of the vertex by its name and title
 464   // if vertexer failed title is empty (default c'tor) and ncontributors is 0
 465   // vtx       name                  title
 466   // Tracks    PrimaryVertex         VertexerTracksNoConstraint
 467   // TPC       TPCVertex             VertexerTracksNoConstraint
 468   // SPD       SPDVertex             vertexer: 3D
 469   // SPD       SPDVertex             vertexer: Z
 470
 471
 472
 473   if(vtx->GetNContributors()<3)return kFALSE;
 474
 475   // do not want tpc only primary vertex
 476   TString vtxName(vtx->GetName());
 477   if(vtxName.Contains("TPCVertex"))return kFALSE;
 478
 479   Float_t zvtx = vtx->GetZ();
 480   Float_t yvtx = vtx->GetY();
 481   Float_t xvtx = vtx->GetX();
 482
 483   // here we may fill the histos for run dependence....
 484
 485   xvtx -= fVtxXMean;
 486   yvtx -= fVtxYMean;
 487   zvtx -= fVtxZMean;
 488
 489   Float_t r2   = yvtx*yvtx+xvtx*xvtx;
 490   Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
 491   return eventSel;
 492 }
 493
 494
 495
 496 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 497 {
 498   // Terminate analysis
 499   //
 500 }