peak finder code from Per Thomas/Yale
[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   fh1Xsec(0x0),
74   fh1Trials(0x0),
75   fh1PtHard(0x0),
76   fh1PtHardTrials(0x0),
77   fh2TriggerCount(0x0),
78   fh2ESDTriggerCount(0x0),
79   fh2TriggerVtx(0x0),
80   fh2ESDTriggerVtx(0x0),
81   fh2ESDTriggerRun(0x0),
82   fh2VtxXY(0x0),
83   fHistList(0x0)  
84 {
85   fRunRange[0] = fRunRange[1] = 0; 
86 }
87
88 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
89   AliAnalysisTaskSE(name),
90   fUseAODInput(kFALSE),
91   fAvgTrials(1),
92   fZVtxCut(8.),
93   fRealData(kFALSE),
94   fh1Xsec(0x0),
95   fh1Trials(0x0),
96   fh1PtHard(0x0),
97   fh1PtHardTrials(0x0),
98   fh2TriggerCount(0x0),
99   fh2ESDTriggerCount(0x0),
100   fh2TriggerVtx(0x0),
101   fh2ESDTriggerVtx(0x0),
102   fh2ESDTriggerRun(0x0),
103   fh2VtxXY(0x0),
104   fHistList(0x0)  
105 {
106   fRunRange[0] = fRunRange[1] = 0; 
107   DefineOutput(1,TList::Class());
108 }
109
110
111
112 Bool_t AliAnalysisTaskJetServices::Notify()
113 {
114   //
115   // Implemented Notify() to read the cross sections
116   // and number of trials from pyxsec.root
117   // 
118
119   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
120   Float_t xsection = 0;
121   Float_t ftrials  = 1;
122
123   fAvgTrials = 1;
124   if(tree){
125     TFile *curfile = tree->GetCurrentFile();
126     if (!curfile) {
127       Error("Notify","No current file");
128       return kFALSE;
129     }
130     if(!fh1Xsec||!fh1Trials){
131       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
132       return kFALSE;
133     }
134     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
135     fh1Xsec->Fill("<#sigma>",xsection);
136     // construct a poor man average trials 
137     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
138     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
139   }  
140   return kTRUE;
141 }
142
143 void AliAnalysisTaskJetServices::UserCreateOutputObjects()
144 {
145
146   //
147   // Create the output container
148   //
149
150
151   if (fDebug > 1) printf("AnalysisTaskJetServices::UserCreateOutputObjects() \n");
152
153   OpenFile(1);
154   if(!fHistList)fHistList = new TList();
155
156   Bool_t oldStatus = TH1::AddDirectoryStatus();
157   TH1::AddDirectory(kFALSE);
158   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
159   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
160   fHistList->Add(fh1Xsec);
161
162   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
163   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
164   fHistList->Add(fh1Trials);
165
166   const Int_t nBinPt = 100;
167   Double_t binLimitsPt[nBinPt+1];
168   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
169     if(iPt == 0){
170       binLimitsPt[iPt] = 0.0;
171     }
172     else {// 1.0
173       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
174     }
175   }
176   
177   fh2TriggerCount = new TH2F("fh2TriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
178   fHistList->Add(fh2TriggerCount);
179
180   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
181   fHistList->Add(fh2ESDTriggerCount);
182
183   fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
184   fHistList->Add(fh2TriggerVtx);
185
186   fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
187   fHistList->Add(fh2ESDTriggerVtx);
188   
189
190   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
191   fHistList->Add(fh1PtHard);
192   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
193   fHistList->Add(fh1PtHardTrials);
194
195   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
196   // 3 triggers BB BE/EB EE
197
198   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);
199   fHistList->Add(fh2ESDTriggerRun);
200
201   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
202   fHistList->Add(fh2VtxXY);
203   // =========== Switch on Sumw2 for all histos ===========
204   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
205     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
206     if (h1){
207       h1->Sumw2();
208       continue;
209     }
210     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
211     if(hn)hn->Sumw2();
212   }
213
214
215   TH1::AddDirectory(oldStatus);
216 }
217
218 void AliAnalysisTaskJetServices::Init()
219 {
220   //
221   // Initialization
222   //
223
224   Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
225   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
226
227 }
228
229 void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
230 {
231
232   //
233   // Execute analysis for current event
234   //
235  
236   AliAODEvent *aod = 0;
237   AliESDEvent *esd = 0;
238
239   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
240
241   if(fUseAODInput){    
242     aod = dynamic_cast<AliAODEvent*>(InputEvent());
243     if(!aod){
244       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
245       return;
246     }
247     // fethc the header
248   }
249   else{
250     //  assume that the AOD is in the general output...
251     aod  = AODEvent();
252     if(!aod){
253       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
254       return;
255     }
256     esd = dynamic_cast<AliESDEvent*>(InputEvent());
257   }
258   
259   if(esd){
260     Float_t run = (Float_t)esd->GetRunNumber();
261     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
262     Float_t zvtx = -999;
263     Float_t xvtx = -999;
264     Float_t yvtx = -999;
265
266     if(vtxESD->GetNContributors()>0){
267       zvtx = vtxESD->GetZ();
268       yvtx = vtxESD->GetY();
269       xvtx = vtxESD->GetX();
270     }
271
272     Int_t iTrig = -1;
273     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
274        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
275        ||esd->GetFiredTriggerClasses().Contains("MB1")
276        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
277       iTrig = 0;
278     }
279     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
280             ||esd->GetFiredTriggerClasses().Contains("CSMBA")
281             ||esd->GetFiredTriggerClasses().Contains("CINT6A")
282             ||esd->GetFiredTriggerClasses().Contains("CINT1C")
283             ||esd->GetFiredTriggerClasses().Contains("CSMBC")
284             ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
285       // empty bunch or bunch empty
286       iTrig = 1;
287     }
288     if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
289        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
290       iTrig = 2;
291     }
292
293     
294     if(iTrig>=0){
295       iTrig *= 3;
296       fh2ESDTriggerRun->Fill(run,iTrig+1);
297       if(vtxESD->GetNContributors()>0){
298         fh2ESDTriggerRun->Fill(run,iTrig+2);
299         fh2VtxXY->Fill(xvtx,yvtx);
300       }
301       if(TMath::Abs(zvtx)<fZVtxCut&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5)fh2ESDTriggerRun->Fill(run,iTrig+3);
302     }
303     else{
304       fh2ESDTriggerRun->Fill(run,0);
305     }
306
307
308   }
309   
310   // loop over all possible trigger and 
311   for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
312     Bool_t esdTrig = kFALSE;
313     if(esd){
314       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
315       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
316     }
317     Bool_t aodTrig = kFALSE;
318     if(aod){
319       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
320       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
321     }
322
323     // Check wether we have also an SPD vertex
324     
325     if(aod){
326       const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
327       //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
328       
329       if(vtxAOD->GetNContributors()>0){
330         if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
331         Float_t zvtx = vtxAOD->GetZ();
332         Float_t yvtx = vtxAOD->GetY();
333         Float_t xvtx = vtxAOD->GetX();
334         fh2TriggerVtx->Fill(it,zvtx);
335         if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
336           fh2TriggerCount->Fill(it,kTriggeredVertexIn);
337         }
338       }
339     }
340     if(esd){
341       const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
342       //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
343       Bool_t cand = fInputHandler->IsEventSelected();
344       if(cand){
345         fh2ESDTriggerCount->Fill(it,kSelectedALICE); 
346         fh2ESDTriggerCount->Fill(it,kSelected);
347         AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
348       }
349       if(vtxESD->GetNContributors()>0){
350         if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
351         Float_t zvtx = vtxESD->GetZ();
352         Float_t yvtx = vtxESD->GetY();
353         Float_t xvtx = vtxESD->GetX();
354         Float_t r2   = yvtx*yvtx+xvtx*xvtx;  
355         fh2ESDTriggerVtx->Fill(it,zvtx);
356         if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&r2<1){
357           fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
358           if(cand){
359             fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
360             fh2ESDTriggerCount->Fill(it,kSelected);
361             AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
362           }
363           // here we select based on ESD info...
364         }
365       }
366     }
367   }
368
369
370
371   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
372
373   
374   Double_t ptHard = 0; 
375   Double_t nTrials = 1; // Trials for MC trigger 
376
377   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
378   AliMCEvent* mcEvent = MCEvent();
379   //    AliStack *pStack = 0; 
380   if(mcEvent){
381     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
382     if(pythiaGenHeader){
383       nTrials = pythiaGenHeader->Trials();
384       ptHard  = pythiaGenHeader->GetPtHard();
385       int iProcessType = pythiaGenHeader->ProcessType();
386       // 11 f+f -> f+f
387       // 12 f+barf -> f+barf
388       // 13 f+barf -> g+g
389       // 28 f+g -> f+g
390       // 53 g+g -> f+barf
391       // 68 g+g -> g+g
392       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
393       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
394       fh1PtHard->Fill(ptHard);
395       fh1PtHardTrials->Fill(ptHard,nTrials);
396
397     }// if pythia gen header
398   }
399
400   // trigger selection
401   
402
403   PostData(1, fHistList);
404 }
405
406
407 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
408 {
409   // Terminate analysis
410   //
411 }