daac60f0f496e8f2d8f2a4555dab6e36aa1b542a
[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(!track)continue;
319       AliESDtrack* trackTPCESD = 0;
320       if(fTrackType==0)
321         trackTPCESD = AliESDtrackCuts::GetTPCOnlyTrack(fESD,track->GetID());
322       else if(fTrackType==1) {
323         trackTPCESD = AliESDtrackCuts::GetTPCOnlyTrack(fESD,track->GetID());
324         if(!trackTPCESD) {
325           continue;
326         }
327       AliExternalTrackParam exParam;
328       Bool_t relate = trackTPCESD->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
329       if( !relate ) {
330         delete trackTPCESD;
331         continue;
332       }
333       trackTPCESD->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
334       }
335       else
336         trackTPCESD = track;
337     
338       if(!trackTPCESD) {
339         delete trackTPCESD;
340         continue;
341       }
342       if(!track) {
343         delete trackTPCESD;
344         continue;
345       }
346
347
348       //fill the container
349       containerInputRec[0] = track->Pt();
350       containerInputRec[1] = track->Phi();
351       containerInputRec[2] = track->Eta();
352     
353       //Store TPC Inner Params for TPConly tracks
354       containerInputTPConly[0] = trackTPCESD->Pt();
355       containerInputTPConly[1] = trackTPCESD->Phi();
356       containerInputTPConly[2] = trackTPCESD->Eta();
357
358       if(trackTPCESD) {
359         if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
360           if(trackTPCESD->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
361           if(trackTPCESD->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
362
363           //Only fill the MC containers if MC information is available
364           if(fMC) {
365             Int_t label = TMath::Abs(track->GetLabel());
366             TParticle *particle = fStack->Particle(label) ;
367             if(!particle) {
368               delete trackTPCESD;
369               continue;
370             }
371             
372             containerInputTPConlyMC[0] = particle->Pt();      
373             containerInputTPConlyMC[1] = particle->Phi();      
374             containerInputTPConlyMC[2] = particle->Eta();  
375             
376             //Container with primaries
377             if(fStack->IsPhysicalPrimary(label)) {
378               if(particle->GetPDG()->Charge()>0.) {
379                 fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConlyMC,kStepReconstructedTPCOnlyMC);
380               }
381               if(particle->GetPDG()->Charge()<0.) {
382                 fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConlyMC,kStepReconstructedTPCOnlyMC);
383               }
384             }
385           }
386           
387         }
388       }
389
390       if (fTrackCuts->AcceptTrack(track)) {
391         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
392         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
393
394         
395         //Only fill the MC containers if MC information is available
396         if(fMC) {
397           Int_t label = TMath::Abs(track->GetLabel());
398           TParticle *particle = fStack->Particle(label) ;
399           if(!particle) {
400             delete trackTPCESD;
401             continue;
402           }
403           containerInputRecMC[0] = particle->Pt();      
404           containerInputRecMC[1] = particle->Phi();      
405           containerInputRecMC[2] = particle->Eta();  
406
407           //Container with primaries
408           if(fStack->IsPhysicalPrimary(label)) {
409             if(particle->GetPDG()->Charge()>0.) {
410               fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
411             }
412             if(particle->GetPDG()->Charge()<0.) {
413               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
414             }
415           }
416
417           //Container with secondaries
418           if (!fStack->IsPhysicalPrimary(label) ) {
419             if(particle->GetPDG()->Charge()>0.) {
420               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
421             }
422             if(particle->GetPDG()->Charge()<0.) {
423               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
424             }
425           }
426         }
427         
428       }//trackCuts
429
430       delete trackTPCESD;
431     }//track loop
432   
433
434   //Fill MC containters if particles are findable
435   if(fMC) {
436     for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
437       AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
438       if(!mcPart) continue;
439       //fill the container
440       containerInputMC[0] = mcPart->Pt();
441       containerInputMC[1] = mcPart->Phi();      
442       containerInputMC[2] = mcPart->Eta();  
443       
444       if(fStack->IsPhysicalPrimary(iPart)) {
445         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
446         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
447       }
448     }
449   }
450   
451   PostData(0,fHistList);
452   PostData(1,fCFManagerPos->GetParticleContainer());
453   PostData(2,fCFManagerNeg->GetParticleContainer());
454   
455 }
456 //________________________________________________________________________
457 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
458   //
459   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
460   // This is to called in Notify and should provide the path to the AOD/ESD file
461   // Copied from AliAnalysisTaskJetSpectrum2
462   //
463
464   TString file(currFile);  
465   fXsec = 0;
466   fTrials = 1;
467
468   if(file.Contains("root_archive.zip#")){
469     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
470     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
471     file.Replace(pos+1,20,"");
472   }
473   else {
474     // not an archive take the basename....
475     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
476   }
477   //  Printf("%s",file.Data());
478    
479   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
480   if(!fxsec){
481     // next trial fetch the histgram file
482     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
483     if(!fxsec){
484         // not a severe condition but inciate that we have no information
485       return kFALSE;
486     }
487     else{
488       // find the tlist we want to be independtent of the name so use the Tkey
489       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
490       if(!key){
491         fxsec->Close();
492         return kFALSE;
493       }
494       TList *list = dynamic_cast<TList*>(key->ReadObj());
495       if(!list){
496         fxsec->Close();
497         return kFALSE;
498       }
499       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
500       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
501       fxsec->Close();
502     }
503   } // no tree pyxsec.root
504   else {
505     TTree *xtree = (TTree*)fxsec->Get("Xsection");
506     if(!xtree){
507       fxsec->Close();
508       return kFALSE;
509     }
510     UInt_t   ntrials  = 0;
511     Double_t  xsection  = 0;
512     xtree->SetBranchAddress("xsection",&xsection);
513     xtree->SetBranchAddress("ntrials",&ntrials);
514     xtree->GetEntry(0);
515     fTrials = ntrials;
516     fXsec = xsection;
517     fxsec->Close();
518   }
519   return kTRUE;
520 }
521 //________________________________________________________________________
522 Bool_t AliPWG4HighPtSpectra::Notify()
523 {
524   //
525   // Implemented Notify() to read the cross sections
526   // and number of trials from pyxsec.root
527   // Copied from AliAnalysisTaskJetSpectrum2
528   // 
529
530   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
531   Float_t xsection = 0;
532   Float_t ftrials  = 1;
533
534   fAvgTrials = 1;
535   if(tree){
536     TFile *curfile = tree->GetCurrentFile();
537     if (!curfile) {
538       Error("Notify","No current file");
539       return kFALSE;
540     }
541     if(!fh1Xsec||!fh1Trials){
542       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
543       return kFALSE;
544     }
545      PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
546     fh1Xsec->Fill("<#sigma>",xsection);
547     // construct a poor man average trials 
548     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
549     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
550   }  
551   return kTRUE;
552 }
553
554 //________________________________________________________________________
555 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
556   
557   if(!mcEvent)return 0;
558   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
559   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
560   if(!pythiaGenHeader){
561     // cocktail ??
562     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
563     
564     if (!genCocktailHeader) {
565       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
566       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
567       return 0;
568     }
569     TList* headerList = genCocktailHeader->GetHeaders();
570     for (Int_t i=0; i<headerList->GetEntries(); i++) {
571       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
572       if (pythiaGenHeader)
573         break;
574     }
575     if(!pythiaGenHeader){
576       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
577       return 0;
578     }
579   }
580   return pythiaGenHeader;
581
582 }
583
584
585 //___________________________________________________________________________
586 void AliPWG4HighPtSpectra::Terminate(Option_t*)
587 {
588   // The Terminate() function is the last function to be called during
589   // a query. It always runs on the client, it can be used to present
590   // the results graphically or save the results to file.
591
592 }
593
594 //___________________________________________________________________________
595 void AliPWG4HighPtSpectra::CreateOutputObjects() {
596   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
597   //TO BE SET BEFORE THE EXECUTION OF THE TASK
598   //
599   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
600
601   Bool_t oldStatus = TH1::AddDirectoryStatus();
602   TH1::AddDirectory(kFALSE); 
603
604   //slot #1
605   OpenFile(0);
606   fHistList = new TList();
607   fHistList->SetOwner(kTRUE);
608   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
609   fHistList->Add(fNEventAll);
610   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
611   fHistList->Add(fNEventSel);
612
613   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
614   fHistList->Add(fNEventReject);
615
616   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
617   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
618   fHistList->Add(fh1Xsec);
619
620   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
621   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
622   fHistList->Add(fh1Trials);
623
624   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
625   fHistList->Add(fh1PtHard);
626   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
627   fHistList->Add(fh1PtHardTrials);
628
629   TH1::AddDirectory(oldStatus);   
630
631   PostData(0,fHistList);
632   PostData(1,fCFManagerPos->GetParticleContainer());
633   PostData(2,fCFManagerNeg->GetParticleContainer());
634
635 }
636
637 #endif