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