Added AliPhysicsSelection
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TRandom.h>
25 #include <TSystem.h>
26 #include <TInterpreter.h>
27 #include <TChain.h>
28 #include <TFile.h>
29 #include <TKey.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TH3F.h>
33 #include <TProfile.h>
34 #include <TList.h>
35 #include <TLorentzVector.h>
36 #include <TClonesArray.h>
37 #include  "TDatabasePDG.h"
38
39 #include "AliAnalysisTaskJetServices.h"
40 #include "AliAnalysisDataContainer.h"
41 #include "AliAnalysisDataSlot.h"
42 #include "AliAnalysisManager.h"
43 #include "AliJetFinder.h"
44 #include "AliJetHeader.h"
45 #include "AliJetReader.h"
46 #include "AliJetReaderHeader.h"
47 #include "AliUA1JetHeaderV1.h"
48 #include "AliESDEvent.h"
49 #include "AliAODEvent.h"
50 #include "AliAODHandler.h"
51 #include "AliAODTrack.h"
52 #include "AliAODJet.h"
53 #include "AliAODMCParticle.h"
54 #include "AliMCEventHandler.h"
55 #include "AliMCEvent.h"
56 #include "AliStack.h"
57 #include "AliGenPythiaEventHeader.h"
58 #include "AliJetKineReaderHeader.h"
59 #include "AliGenCocktailEventHeader.h"
60 #include "AliInputEventHandler.h"
61 #include "AliPhysicsSelection.h"
62
63
64 #include "AliAnalysisHelperJetTasks.h"
65
66 ClassImp(AliAnalysisTaskJetServices)
67
68 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
69   fUseAODInput(kFALSE),
70   fAvgTrials(1),
71   fZVtxCut(8.),
72   fRealData(kFALSE),
73   fPhysicsSelection(0x0),
74   fh1Xsec(0x0),
75   fh1Trials(0x0),
76   fh1PtHard(0x0),
77   fh1PtHardTrials(0x0),
78   fh2TriggerCount(0x0),
79   fh2ESDTriggerCount(0x0),
80   fh2TriggerVtx(0x0),
81   fh2ESDTriggerVtx(0x0),
82   fh2ESDTriggerRun(0x0),
83   fh2VtxXY(0x0),
84   fHistList(0x0)  
85 {
86   fRunRange[0] = fRunRange[1] = 0; 
87 }
88
89 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
90   AliAnalysisTaskSE(name),
91   fUseAODInput(kFALSE),
92   fAvgTrials(1),
93   fZVtxCut(8.),
94   fRealData(kFALSE),
95   fPhysicsSelection(0x0),
96   fh1Xsec(0x0),
97   fh1Trials(0x0),
98   fh1PtHard(0x0),
99   fh1PtHardTrials(0x0),
100   fh2TriggerCount(0x0),
101   fh2ESDTriggerCount(0x0),
102   fh2TriggerVtx(0x0),
103   fh2ESDTriggerVtx(0x0),
104   fh2ESDTriggerRun(0x0),
105   fh2VtxXY(0x0),
106   fHistList(0x0)  
107 {
108   fRunRange[0] = fRunRange[1] = 0; 
109   DefineOutput(1,TList::Class());
110 }
111
112
113
114 Bool_t AliAnalysisTaskJetServices::Notify()
115 {
116   //
117   // Implemented Notify() to read the cross sections
118   // and number of trials from pyxsec.root
119   // 
120
121   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
122   Float_t xsection = 0;
123   Float_t ftrials  = 1;
124
125   fAvgTrials = 1;
126   if(tree){
127     TFile *curfile = tree->GetCurrentFile();
128     if (!curfile) {
129       Error("Notify","No current file");
130       return kFALSE;
131     }
132     if(!fh1Xsec||!fh1Trials){
133       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
134       return kFALSE;
135     }
136     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
137     fh1Xsec->Fill("<#sigma>",xsection);
138     // construct a poor man average trials 
139     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
140     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
141   }  
142   return kTRUE;
143 }
144
145 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
146 {
147
148   //
149   // Create the output container
150   //
151
152
153   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
154
155   OpenFile(1);
156   if(!fHistList)fHistList = new TList();
157
158   Bool_t oldStatus = TH1::AddDirectoryStatus();
159   TH1::AddDirectory(kFALSE);
160
161   if(!fPhysicsSelection)
162     fPhysicsSelection = new AliPhysicsSelection();
163   fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back
164   //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
165   fHistList->Add(fPhysicsSelection);
166
167   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
168   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
169   fHistList->Add(fh1Xsec);
170
171   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
172   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
173   fHistList->Add(fh1Trials);
174
175   const Int_t nBinPt = 100;
176   Double_t binLimitsPt[nBinPt+1];
177   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
178     if(iPt == 0){
179       binLimitsPt[iPt] = 0.0;
180     }
181     else {// 1.0
182       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
183     }
184   }
185   
186   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
187   fHistList->Add(fh2TriggerCount);
188
189   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
190   fHistList->Add(fh2ESDTriggerCount);
191
192   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
193   fHistList->Add(fh2TriggerVtx);
194
195   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
196   fHistList->Add(fh2ESDTriggerVtx);
197   
198
199   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
200   fHistList->Add(fh1PtHard);
201   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
202   fHistList->Add(fh1PtHardTrials);
203
204   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
205   // 3 triggers BB BE/EB EE
206
207   fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
208   fHistList->Add(fh2ESDTriggerRun);
209
210   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
211   fHistList->Add(fh2VtxXY);
212   // =========== Switch on Sumw2 for all histos ===========
213   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
214     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
215     if (h1){
216       h1->Sumw2();
217       continue;
218     }
219     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
220     if(hn)hn->Sumw2();
221   }
222
223
224   TH1::AddDirectory(oldStatus);
225 }
226
227 void AliAnalysisTaskJetServices::Init()
228 {
229   //
230   // Initialization
231   //
232
233   Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
234   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
235
236 }
237
238 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
239 {
240
241   //
242   // Execute analysis for current event
243   //
244  
245   AliAODEvent *aod = 0;
246   AliESDEvent *esd = 0;
247
248   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
249
250   if(fUseAODInput){    
251     aod = dynamic_cast<AliAODEvent*>(InputEvent());
252     if(!aod){
253       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
254       return;
255     }
256     // fethc the header
257   }
258   else{
259     //  assume that the AOD is in the general output...
260     aod  = AODEvent();
261     if(!aod){
262       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
263       return;
264     }
265     esd = dynamic_cast<AliESDEvent*>(InputEvent());
266   }
267   
268   if(esd){
269     Float_t run = (Float_t)esd->GetRunNumber();
270     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
271     Float_t zvtx = -999;
272     Float_t xvtx = -999;
273     Float_t yvtx = -999;
274
275     if(vtxESD->GetNContributors()>0){
276       zvtx = vtxESD->GetZ();
277       yvtx = vtxESD->GetY();
278       xvtx = vtxESD->GetX();
279     }
280
281     Int_t iTrig = -1;
282     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
283        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
284        ||esd->GetFiredTriggerClasses().Contains("MB1")
285        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
286       iTrig = 0;
287     }
288     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
289             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
290             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
291             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
292             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
293             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
294       // empty bunch or bunch empty
295       iTrig = 1;
296     }
297     if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
298        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
299       iTrig = 2;
300     }
301
302     
303     if(iTrig>=0){
304       iTrig *= 3;
305       fh2ESDTriggerRun->Fill(run,iTrig+1);
306       if(vtxESD->GetNContributors()>0){
307         fh2ESDTriggerRun->Fill(run,iTrig+2);
308         fh2VtxXY->Fill(xvtx,yvtx);
309       }
310       if(TMath::Abs(zvtx)<fZVtxCut&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5)fh2ESDTriggerRun->Fill(run,iTrig+3);
311     }
312     else{
313       fh2ESDTriggerRun->Fill(run,0);
314     }
315
316
317   }
318   
319   // loop over all possible trigger and 
320   for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
321     Bool_t esdTrig = kFALSE;
322     if(esd){
323       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
324       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
325     }
326     Bool_t aodTrig = kFALSE;
327     if(aod){
328       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
329       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
330     }
331
332     // Check wether we have also an SPD vertex
333     
334     if(aod){
335       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
336       //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
337       
338       if(vtxAOD->GetNContributors()>0){
339         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
340         Float_t zvtx = vtxAOD->GetZ();
341         Float_t yvtx = vtxAOD->GetY();
342         Float_t xvtx = vtxAOD->GetX();
343         fh2TriggerVtx->Fill(it,zvtx);
344         if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
345           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
346         }
347       }
348     }
349     if(esd){
350       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
351       //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
352       Bool_t cand = fPhysicsSelection->IsCollisionCandidate(esd);
353       if(cand) fh2ESDTriggerCount->Fill(it,kSelectedALICE);
354       if(vtxESD->GetNContributors()>0){
355         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
356         Float_t zvtx = vtxESD->GetZ();
357         Float_t yvtx = vtxESD->GetY();
358         Float_t xvtx = vtxESD->GetX();
359         fh2ESDTriggerVtx->Fill(it,zvtx);
360         if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
361           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
362           // here we select based on ESD info...
363           if(fRealData){
364             if(cand){
365               fh2ESDTriggerCount->Fill(it,kSelected);
366               AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
367             }
368           }
369           else{
370             fh2ESDTriggerCount->Fill(it,kSelected);
371             AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
372           }
373         }
374
375
376       }
377
378     }
379
380   }
381
382
383
384   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
385
386   
387   Double_t ptHard = 0; 
388   Double_t nTrials = 1; // Trials for MC trigger 
389
390   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
391   AliMCEvent* mcEvent = MCEvent();
392   //    AliStack *pStack = 0; 
393   if(mcEvent){
394     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
395     if(pythiaGenHeader){
396       nTrials = pythiaGenHeader->Trials();
397       ptHard  = pythiaGenHeader->GetPtHard();
398       int iProcessType = pythiaGenHeader->ProcessType();
399       // 11 f+f -> f+f
400       // 12 f+barf -> f+barf
401       // 13 f+barf -> g+g
402       // 28 f+g -> f+g
403       // 53 g+g -> f+barf
404       // 68 g+g -> g+g
405       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
406       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
407       fh1PtHard->Fill(ptHard);
408       fh1PtHardTrials->Fill(ptHard,nTrials);
409
410     }// if pythia gen header
411   }
412
413   // trigger selection
414   
415
416   PostData(1, fHistList);
417 }
418
419
420 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
421 {
422   // Terminate analysis
423   //
424
425   TDirectory* owd = gDirectory;
426
427   if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
428
429   fHistList = dynamic_cast<TList*> (GetOutputData(1));
430   if (!fHistList)
431     Printf("ERROR: fOutput not available");
432
433
434
435   AliAnalysisDataContainer *cont = GetOutputSlot(1)->GetContainer();
436   TString filename = cont->GetFileName();
437   TFile *f = NULL;
438   // Check first if the file is already opened
439   f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
440   if (f) {
441     // Cd to file
442     f->cd();
443     // Check for a folder request
444     TString dir = cont->GetFolderName(); 
445     if (!dir.IsNull()) {
446       if (!f->GetDirectory(dir)) f->mkdir(dir);
447       f->cd(dir);
448     }
449   }
450
451   if (fHistList)
452   {
453     fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (fHistList->FindObject("AliPhysicsSelection_outputlist"));
454   }
455
456   if (fPhysicsSelection)
457     {
458       fPhysicsSelection->SaveHistograms("physics_selection");
459       fPhysicsSelection->Print();
460     }
461
462   owd->cd();
463
464 }