PWG4/JetTasks/AliPWG4HighPtSpectra.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // Example of task (running locally, on AliEn and CAF),
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 //-----------------------------------------------------------------------
  23 // Author : Marta Verweij - UU
  24 //-----------------------------------------------------------------------
  25
  26
  27 #ifndef ALIPWG4HIGHPTSPECTRA_CXX
  28 #define ALIPWG4HIGHPTSPECTRA_CXX
  29
  30 #include "AliPWG4HighPtSpectra.h"
  31
  32 #include "TVector3.h"
  33 #include <iostream>
  34 #include "TH1.h"
  35 #include "TH2.h"
  36 #include "TH3.h"
  37 #include "TProfile.h"
  38 #include "TList.h"
  39 #include "TChain.h"
  40 #include "TKey.h"
  41 #include "TSystem.h"
  42 #include "TFile.h"
  43
  44 #include "AliAnalysisManager.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliExternalTrackParam.h"
  49
  50 #include "AliLog.h"
  51
  52 #include "AliStack.h"
  53 #include "TParticle.h"
  54 #include "AliMCEvent.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliCFContainer.h"
  57 #include "AliGenPythiaEventHeader.h"
  58 #include "AliGenCocktailEventHeader.h"
  59 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
  60
  61 //#include <iostream>
  62 using namespace std; //required for resolving the 'cout' symbol
  63
  64 ClassImp(AliPWG4HighPtSpectra)
  65
  66 //__________________________________________________________________________
  67 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
  68   fReadAODData(0),
  69   fCFManagerPos(0x0),
  70   fCFManagerNeg(0x0),
  71   fESD(0x0),
  72   fMC(0x0),
  73   fStack(0x0),
  74   fTrackCuts(0x0),
  75   fTrackCutsTPConly(0x0),
  76   fAvgTrials(1),
  77   fHistList(0),
  78   fNEventAll(0),
  79   fNEventSel(0),
  80   fNEventReject(0),
  81   fh1Xsec(0),
  82   fh1Trials(0),
  83   fh1PtHard(0),
  84   fh1PtHardTrials(0)
  85 {
  86   //
  87   //Default ctor
  88   //
  89 }
  90 //___________________________________________________________________________
  91 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
  92   AliAnalysisTask(name,""),
  93   fReadAODData(0),
  94   fCFManagerPos(0x0),
  95   fCFManagerNeg(0x0),
  96   fESD(0x0),
  97   fMC(0x0),
  98   fStack(0x0),
  99   fTrackCuts(0x0),
 100   fTrackCutsTPConly(0x0),
 101   fAvgTrials(1),
 102   fHistList(0),
 103   fNEventAll(0),
 104   fNEventSel(0),
 105   fNEventReject(0),
 106   fh1Xsec(0),
 107   fh1Trials(0),
 108   fh1PtHard(0),
 109   fh1PtHardTrials(0)
 110 {
 111   //
 112   // Constructor. Initialization of Inputs and Outputs
 113   //
 114   AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
 115   // Input slot #0 works with a TChain ESD
 116   DefineInput(0, TChain::Class());
 117   // Output slot #0 writes into a TList
 118   DefineOutput(0,TList::Class());
 119   // Output slot #1, #2 writes into a AliCFContainer
 120   DefineOutput(1,AliCFContainer::Class());
 121   DefineOutput(2,AliCFContainer::Class());
 122   // Output slot #3 writes into a AliESDtrackCuts
 123   DefineOutput(3, AliESDtrackCuts::Class());
 124   DefineOutput(4, AliESDtrackCuts::Class());
 125 }
 126
 127 //________________________________________________________________________
 128 void AliPWG4HighPtSpectra::LocalInit()
 129 {
 130   //
 131   // Only called once at beginning
 132   //
 133   PostData(3,fTrackCuts);
 134   PostData(4,fTrackCutsTPConly);
 135 }
 136
 137 //________________________________________________________________________
 138 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
 139 {
 140   // Connect ESD here
 141   // Called once
 142   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 143
 144   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 145   if (!tree) {
 146     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 147     return;
 148   }
 149
 150   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 151
 152   if (!esdH) {
 153     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 154     return;
 155   } else
 156     fESD = esdH->GetEvent();
 157
 158   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 159   if (!eventHandler) {
 160     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 161   }
 162   else
 163     fMC = eventHandler->MCEvent();
 164
 165 }
 166
 167 //________________________________________________________________________
 168 Bool_t AliPWG4HighPtSpectra::SelectEvent() {
 169   //
 170   // Decide if event should be selected for analysis
 171   //
 172
 173   // Checks following requirements:
 174   // - fESD available
 175   // - trigger info from AliPhysicsSelection
 176   // - number of reconstructed tracks > 1
 177   // - primary vertex reconstructed
 178   // - z-vertex < 10 cm
 179
 180   Bool_t selectEvent = kTRUE;
 181
 182   //fESD object available?
 183   if (!fESD) {
 184     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 185     fNEventReject->Fill("noESD",1);
 186     selectEvent = kFALSE;
 187     return selectEvent;
 188   }
 189
 190   //Trigger
 191   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 192   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 193     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 194     fNEventReject->Fill("Trigger",1);
 195     selectEvent = kFALSE;
 196     return selectEvent;
 197   }
 198
 199   //Check if number of reconstructed tracks is larger than 1
 200   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 201     fNEventReject->Fill("NTracks<2",1);
 202     selectEvent = kFALSE;
 203     return selectEvent;
 204   }
 205
 206   //Check if vertex is reconstructed
 207   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 208   // Need vertex cut
 209   TString vtxName(vtx->GetName());
 210   if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 211     // SPD vertex
 212     vtx = fESD->GetPrimaryVertexSPD();
 213     if(vtx->GetNContributors()<2) {
 214       vtx = 0x0;
 215       fNEventReject->Fill("noVTX",1);
 216       selectEvent = kFALSE;
 217       return selectEvent;
 218     }
 219   }
 220
 221   //Check if z-vertex < 10 cm
 222   double primVtx[3];
 223   vtx->GetXYZ(primVtx);
 224   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 225     fNEventReject->Fill("ZVTX>10",1);
 226     selectEvent = kFALSE;
 227     return selectEvent;
 228   }
 229
 230   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 231
 232   return selectEvent;
 233
 234 }
 235
 236 //_________________________________________________
 237 void AliPWG4HighPtSpectra::Exec(Option_t *)
 238 {
 239   //
 240   // Main loop function
 241   //
 242   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
 243
 244   // All events without selection
 245   fNEventAll->Fill(0.);
 246
 247   if(!SelectEvent()) {
 248     fNEventReject->Fill("NTracks<2",1);
 249     // Post output data
 250     PostData(0,fHistList);
 251     PostData(1,fCFManagerPos->GetParticleContainer());
 252     PostData(2,fCFManagerNeg->GetParticleContainer());
 253     return;
 254   }
 255
 256   //MCEvent available?
 257   //if yes: get stack
 258   if(fMC) {
 259     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 260     fStack = fMC->Stack();                //Particles Stack
 261     AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
 262   }
 263
 264   // ---- Get MC Header information (for MC productions in pThard bins) ----
 265   Double_t ptHard = 0.;
 266   Double_t nTrials = 1; // trials for MC trigger weight for real data
 267
 268   if(fMC){
 269     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 270      if(pythiaGenHeader){
 271        nTrials = pythiaGenHeader->Trials();
 272        ptHard  = pythiaGenHeader->GetPtHard();
 273
 274        fh1PtHard->Fill(ptHard);
 275        fh1PtHardTrials->Fill(ptHard,nTrials);
 276
 277        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 278      }
 279   }
 280
 281   Int_t nTracks = fESD->GetNumberOfTracks();
 282   AliDebug(2,Form("nTracks %d", nTracks));
 283
 284   if(!fTrackCuts) {
 285     fNEventReject->Fill("noTrackCuts",1);
 286     // Post output data
 287     PostData(0,fHistList);
 288     PostData(1,fCFManagerPos->GetParticleContainer());
 289     PostData(2,fCFManagerNeg->GetParticleContainer());
 290     return;
 291   }
 292
 293   // Selected events for analysis
 294   fNEventSel->Fill(0.);
 295
 296
 297   Double_t containerInputRec[3]       = {0.,0.,0.};
 298   Double_t containerInputTPConly[3]   = {0.,0.,0.};
 299   Double_t containerInputMC[3]        = {0.,0.,0.};
 300   Double_t containerInputRecMC[3]     = {0.,0.,0.};
 301   Double_t containerInputTPConlyMC[3] = {0.,0.,0.};
 302
 303   //Now go to rec level
 304   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 305     {
 306       if(!fESD->GetTrack(iTrack) ) continue;
 307       AliESDtrack* track = fESD->GetTrack(iTrack);
 308       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
 309       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
 310       if(!track || !trackTPC) continue;
 311
 312       //fill the container
 313       containerInputRec[0] = track->Pt();
 314       containerInputRec[1] = track->Phi();
 315       containerInputRec[2] = track->Eta();
 316
 317       //Store TPC Inner Params for TPConly tracks
 318       containerInputTPConly[0] = trackTPC->Pt();
 319       containerInputTPConly[1] = trackTPC->Phi();
 320       containerInputTPConly[2] = trackTPC->Eta();
 321
 322       AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
 323       if(trackTPCESD) {
 324         if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
 325           if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 326           if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 327
 328           //Only fill the MC containers if MC information is available
 329           if(fMC) {
 330             Int_t label = TMath::Abs(track->GetLabel());
 331             TParticle *particle = fStack->Particle(label) ;
 332             if(!particle) continue;
 333
 334             containerInputTPConlyMC[0] = particle->Pt();
 335             containerInputTPConlyMC[1] = particle->Phi();
 336             containerInputTPConlyMC[2] = particle->Eta();
 337
 338             //Container with primaries
 339             if(fStack->IsPhysicalPrimary(label)) {
 340               if(particle->GetPDG()->Charge()>0.) {
 341                 fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConlyMC,kStepReconstructedTPCOnlyMC);
 342               }
 343               if(particle->GetPDG()->Charge()<0.) {
 344                 fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConlyMC,kStepReconstructedTPCOnlyMC);
 345               }
 346             }
 347           }
 348
 349         }
 350       }
 351
 352       if (fTrackCuts->AcceptTrack(track)) {
 353         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 354         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 355
 356
 357         //Only fill the MC containers if MC information is available
 358         if(fMC) {
 359           Int_t label = TMath::Abs(track->GetLabel());
 360           TParticle *particle = fStack->Particle(label) ;
 361           if(!particle) continue;
 362
 363           containerInputRecMC[0] = particle->Pt();
 364           containerInputRecMC[1] = particle->Phi();
 365           containerInputRecMC[2] = particle->Eta();
 366
 367           //Container with primaries
 368           if(fStack->IsPhysicalPrimary(label)) {
 369             if(particle->GetPDG()->Charge()>0.) {
 370               fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 371             }
 372             if(particle->GetPDG()->Charge()<0.) {
 373               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 374             }
 375           }
 376
 377           //Container with secondaries
 378           if (!fStack->IsPhysicalPrimary(label) ) {
 379             if(particle->GetPDG()->Charge()>0.) {
 380               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 381             }
 382             if(particle->GetPDG()->Charge()<0.) {
 383               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 384             }
 385           }
 386         }
 387
 388       }//trackCuts
 389
 390       delete trackTPCESD;
 391     }//track loop
 392
 393
 394   //Fill MC containters if particles are findable
 395   if(fMC) {
 396     for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
 397       AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
 398       if(!mcPart) continue;
 399       //fill the container
 400       containerInputMC[0] = mcPart->Pt();
 401       containerInputMC[1] = mcPart->Phi();
 402       containerInputMC[2] = mcPart->Eta();
 403
 404       if(fStack->IsPhysicalPrimary(iPart)) {
 405         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 406         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 407       }
 408     }
 409   }
 410
 411   PostData(0,fHistList);
 412   PostData(1,fCFManagerPos->GetParticleContainer());
 413   PostData(2,fCFManagerNeg->GetParticleContainer());
 414
 415 }
 416 //________________________________________________________________________
 417 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 418   //
 419   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 420   // This is to called in Notify and should provide the path to the AOD/ESD file
 421   // Copied from AliAnalysisTaskJetSpectrum2
 422   //
 423
 424   TString file(currFile);
 425   fXsec = 0;
 426   fTrials = 1;
 427
 428   if(file.Contains("root_archive.zip#")){
 429     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 430     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 431     file.Replace(pos+1,20,"");
 432   }
 433   else {
 434     // not an archive take the basename....
 435     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 436   }
 437   //  Printf("%s",file.Data());
 438
 439   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 440   if(!fxsec){
 441     // next trial fetch the histgram file
 442     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 443     if(!fxsec){
 444         // not a severe condition but inciate that we have no information
 445       return kFALSE;
 446     }
 447     else{
 448       // find the tlist we want to be independtent of the name so use the Tkey
 449       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 450       if(!key){
 451         fxsec->Close();
 452         return kFALSE;
 453       }
 454       TList *list = dynamic_cast<TList*>(key->ReadObj());
 455       if(!list){
 456         fxsec->Close();
 457         return kFALSE;
 458       }
 459       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 460       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 461       fxsec->Close();
 462     }
 463   } // no tree pyxsec.root
 464   else {
 465     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 466     if(!xtree){
 467       fxsec->Close();
 468       return kFALSE;
 469     }
 470     UInt_t   ntrials  = 0;
 471     Double_t  xsection  = 0;
 472     xtree->SetBranchAddress("xsection",&xsection);
 473     xtree->SetBranchAddress("ntrials",&ntrials);
 474     xtree->GetEntry(0);
 475     fTrials = ntrials;
 476     fXsec = xsection;
 477     fxsec->Close();
 478   }
 479   return kTRUE;
 480 }
 481 //________________________________________________________________________
 482 Bool_t AliPWG4HighPtSpectra::Notify()
 483 {
 484   //
 485   // Implemented Notify() to read the cross sections
 486   // and number of trials from pyxsec.root
 487   // Copied from AliAnalysisTaskJetSpectrum2
 488   //
 489
 490   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 491   Float_t xsection = 0;
 492   Float_t ftrials  = 1;
 493
 494   fAvgTrials = 1;
 495   if(tree){
 496     TFile *curfile = tree->GetCurrentFile();
 497     if (!curfile) {
 498       Error("Notify","No current file");
 499       return kFALSE;
 500     }
 501     if(!fh1Xsec||!fh1Trials){
 502       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 503       return kFALSE;
 504     }
 505      PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 506     fh1Xsec->Fill("<#sigma>",xsection);
 507     // construct a poor man average trials
 508     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 509     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 510   }
 511   return kTRUE;
 512 }
 513
 514 //________________________________________________________________________
 515 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
 516
 517   if(!mcEvent)return 0;
 518   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 519   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 520   if(!pythiaGenHeader){
 521     // cocktail ??
 522     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 523
 524     if (!genCocktailHeader) {
 525       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 526       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 527       return 0;
 528     }
 529     TList* headerList = genCocktailHeader->GetHeaders();
 530     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 531       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 532       if (pythiaGenHeader)
 533         break;
 534     }
 535     if(!pythiaGenHeader){
 536       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 537       return 0;
 538     }
 539   }
 540   return pythiaGenHeader;
 541
 542 }
 543
 544
 545 //___________________________________________________________________________
 546 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 547 {
 548   // The Terminate() function is the last function to be called during
 549   // a query. It always runs on the client, it can be used to present
 550   // the results graphically or save the results to file.
 551
 552 }
 553
 554 //___________________________________________________________________________
 555 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 556   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 557   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 558   //
 559   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 560
 561   Bool_t oldStatus = TH1::AddDirectoryStatus();
 562   TH1::AddDirectory(kFALSE);
 563
 564   //slot #1
 565   OpenFile(0);
 566   fHistList = new TList();
 567   fHistList->SetOwner(kTRUE);
 568   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 569   fHistList->Add(fNEventAll);
 570   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 571   fHistList->Add(fNEventSel);
 572
 573   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 574   fHistList->Add(fNEventReject);
 575
 576   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 577   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 578   fHistList->Add(fh1Xsec);
 579
 580   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 581   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 582   fHistList->Add(fh1Trials);
 583
 584   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 585   fHistList->Add(fh1PtHard);
 586   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 587   fHistList->Add(fh1PtHardTrials);
 588
 589   TH1::AddDirectory(oldStatus);
 590
 591 }
 592
 593 #endif