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