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