Switched on vertex constrained again for Jet SPectrum, limit some printout in three...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.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 "AliAnalysisTaskJetSpectrum2.h"
40 #include "AliAnalysisManager.h"
41 #include "AliJetFinder.h"
42 #include "AliJetHeader.h"
43 #include "AliJetReader.h"
44 #include "AliJetReaderHeader.h"
45 #include "AliUA1JetHeaderV1.h"
46 #include "AliESDEvent.h"
47 #include "AliAODEvent.h"
48 #include "AliAODHandler.h"
49 #include "AliAODTrack.h"
50 #include "AliAODJet.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
54 #include "AliStack.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
59
60
61 #include "AliAnalysisHelperJetTasks.h"
62
63 ClassImp(AliAnalysisTaskJetSpectrum2)
64
65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
66   fJetHeaderRec(0x0),
67   fJetHeaderGen(0x0),
68   fAOD(0x0),
69   fhnCorrelation(0x0),
70   fBranchRec("jets"),
71   fBranchGen(""),
72   fUseAODInput(kFALSE),
73   fUseGlobalSelection(kFALSE),
74   fUseExternalWeightOnly(kFALSE),
75   fLimitGenJetEta(kFALSE),
76   fFilterMask(0),
77   fAnalysisType(0),
78   fTrackTypeRec(kTrackUndef),
79   fTrackTypeGen(kTrackUndef),
80   fAvgTrials(1),
81   fExternalWeight(1),    
82   fRecEtaWindow(0.5),
83   fh1Xsec(0x0),
84   fh1Trials(0x0),
85   fh1PtHard(0x0),
86   fh1PtHardNoW(0x0),  
87   fh1PtHardTrials(0x0),
88   fh1NGenJets(0x0),
89   fh1NRecJets(0x0),
90   fh1PtTrackRec(0x0),   
91   fh1SumPtTrackRec(0x0),  
92   fh1SumPtTrackAreaRec(0x0),  
93   fHistList(0x0)  
94 {
95   for(int i = 0;i < kMaxStep*2;++i){
96     fhnJetContainer[i] = 0;
97   }
98   for(int i = 0;i < kMaxJets;++i){
99     fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
100     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
101   }  
102
103 }
104
105 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
106   AliAnalysisTaskSE(name),
107   fJetHeaderRec(0x0),
108   fJetHeaderGen(0x0),
109   fAOD(0x0),
110   fhnCorrelation(0x0),
111   fBranchRec("jets"),
112   fBranchGen(""),
113   fUseAODInput(kFALSE),
114   fUseGlobalSelection(kFALSE),
115   fUseExternalWeightOnly(kFALSE),
116   fLimitGenJetEta(kFALSE),
117   fFilterMask(0),
118   fAnalysisType(0),
119   fTrackTypeRec(kTrackUndef),
120   fTrackTypeGen(kTrackUndef),
121   fAvgTrials(1),
122   fExternalWeight(1),    
123   fRecEtaWindow(0.5),
124   fh1Xsec(0x0),
125   fh1Trials(0x0),
126   fh1PtHard(0x0),
127   fh1PtHardNoW(0x0),  
128   fh1PtHardTrials(0x0),
129   fh1NGenJets(0x0),
130   fh1NRecJets(0x0),
131   fh1PtTrackRec(0x0),   
132   fh1SumPtTrackRec(0x0),  
133   fh1SumPtTrackAreaRec(0x0),  
134   fHistList(0x0)
135 {
136
137   for(int i = 0;i < kMaxStep*2;++i){
138     fhnJetContainer[i] = 0;
139   }  
140   for(int i = 0;i < kMaxJets;++i){
141     fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
142     fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
143   }
144   DefineOutput(1, TList::Class());  
145 }
146
147
148
149 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
150 {
151   //
152   // Implemented Notify() to read the cross sections
153   // and number of trials from pyxsec.root
154   // 
155
156   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
157   Float_t xsection = 0;
158   Float_t ftrials  = 1;
159
160   fAvgTrials = 1;
161   if(tree){
162     TFile *curfile = tree->GetCurrentFile();
163     if (!curfile) {
164       Error("Notify","No current file");
165       return kFALSE;
166     }
167     if(!fh1Xsec||!fh1Trials){
168       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
169       return kFALSE;
170     }
171     AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
172     fh1Xsec->Fill("<#sigma>",xsection);
173     // construct a poor man average trials 
174     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
175     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
176   }  
177   return kTRUE;
178 }
179
180 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
181 {
182
183   //
184   // Create the output container
185   //
186
187
188   // Connect the AOD
189
190
191   MakeJetContainer();
192
193
194   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
195
196   OpenFile(1);
197   if(!fHistList)fHistList = new TList();
198
199   Bool_t oldStatus = TH1::AddDirectoryStatus();
200   TH1::AddDirectory(kFALSE);
201
202   //
203   //  Histogram
204     
205   const Int_t nBinPt = 100;
206   Double_t binLimitsPt[nBinPt+1];
207   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
208     if(iPt == 0){
209       binLimitsPt[iPt] = 0.0;
210     }
211     else {// 1.0
212       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
213     }
214   }
215   
216   const Int_t nBinPhi = 90;
217   Double_t binLimitsPhi[nBinPhi+1];
218   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
219     if(iPhi==0){
220       binLimitsPhi[iPhi] = -0.5*TMath::Pi();
221     }
222     else{
223       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
224     }
225   }
226
227
228
229   const Int_t nBinEta = 40;
230   Double_t binLimitsEta[nBinEta+1];
231   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
232     if(iEta==0){
233       binLimitsEta[iEta] = -2.0;
234     }
235     else{
236       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 1/(Float_t)nBinEta + 0.1;
237     }
238   }
239
240   const Int_t nBinFrag = 25;
241
242
243   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
244   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
245
246   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
247   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
248
249   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
250   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
251   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
252
253   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
254   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
255
256   fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
257   fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
258   fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
259
260   // 
261   for(int ij = 0;ij < kMaxJets;++ij){
262     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
263     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
264
265     fh2PhiPt[ij] =  new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
266                              nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
267
268     fh2PhiEta[ij] =  new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;p_{T,jet}",
269                               nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
270
271
272     fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
273                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
274     fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
275                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
276
277     fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
278                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
279     fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
280                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
281   }
282
283
284   const Int_t saveLevel = 3; // large save level more histos
285   if(saveLevel>0){
286     fHistList->Add(fh1Xsec);
287     fHistList->Add(fh1Trials);
288     fHistList->Add(fh1PtHard);
289     fHistList->Add(fh1PtHardNoW);
290     fHistList->Add(fh1PtHardTrials);
291     fHistList->Add(fh1NGenJets);
292     fHistList->Add(fh1NRecJets);
293     fHistList->Add(fh1PtTrackRec);
294     fHistList->Add(fh1SumPtTrackRec);
295     fHistList->Add(fh1SumPtTrackAreaRec);
296     for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
297     for(int ij = 0;ij<kMaxJets;++ij){
298       fHistList->Add( fh1PtRecIn[ij]);
299       fHistList->Add( fh1PtGenIn[ij]);
300       fHistList->Add( fh2PhiPt[ij]);
301       fHistList->Add( fh2PhiEta[ij]);
302       fHistList->Add( fh2FragRec[ij]);
303       fHistList->Add( fh2FragLnRec[ij]);
304       fHistList->Add( fh2FragGen[ij]);
305       fHistList->Add( fh2FragLnGen[ij]);
306     }
307     fHistList->Add(fhnCorrelation);
308   }
309
310   // =========== Switch on Sumw2 for all histos ===========
311   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
312     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
313     if (h1){
314       h1->Sumw2();
315       continue;
316     }
317     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
318     if(hn)hn->Sumw2();
319   }
320   TH1::AddDirectory(oldStatus);
321 }
322
323 void AliAnalysisTaskJetSpectrum2::Init()
324 {
325   //
326   // Initialization
327   //
328
329   Printf(">>> AnalysisTaskJetSpectrum2::Init() debug level %d\n",fDebug);
330   if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
331
332 }
333
334 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
335 {
336
337   if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
338     // no selection by the service task, we continue
339     if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
340     PostData(1, fHistList);
341     return;
342   }
343
344   //
345   // Execute analysis for current event
346   //
347   AliESDEvent *fESD = 0;
348   if(fUseAODInput){    
349     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
350     if(!fAOD){
351       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
352       return;
353     }
354     // fethc the header
355   }
356   else{
357     //  assume that the AOD is in the general output...
358     fAOD  = AODEvent();
359     if(!fAOD){
360       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
361       return;
362     }
363     if(fDebug>0){
364       fESD = dynamic_cast<AliESDEvent*> (InputEvent());
365     }
366   }
367   
368
369
370
371   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
372
373   
374   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
375
376   if(!aodH){
377     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
378     return;
379   }
380
381   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
382   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
383   if(!aodRecJets){
384     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
385     return;
386   }
387
388   // ==== General variables needed
389
390
391   // We use statice array, not to fragment the memory
392   AliAODJet genJets[kMaxJets];
393   Int_t nGenJets = 0;
394   AliAODJet recJets[kMaxJets];
395   Int_t nRecJets = 0;
396   ///////////////////////////
397
398
399   Double_t eventW = 1;
400   Double_t ptHard = 0; 
401   Double_t nTrials = 1; // Trials for MC trigger 
402
403   if(fUseExternalWeightOnly){
404     eventW = fExternalWeight;
405   }
406
407   fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
408
409   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
410   if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
411     // this is the part we only use when we have MC information
412     AliMCEvent* mcEvent = MCEvent();
413     //    AliStack *pStack = 0; 
414     if(!mcEvent){
415       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
416       return;
417     }
418     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
419     Int_t iCount = 0;  
420     if(pythiaGenHeader){
421       nTrials = pythiaGenHeader->Trials();
422       ptHard  = pythiaGenHeader->GetPtHard();
423       int iProcessType = pythiaGenHeader->ProcessType();
424       // 11 f+f -> f+f
425       // 12 f+barf -> f+barf
426       // 13 f+barf -> g+g
427       // 28 f+g -> f+g
428       // 53 g+g -> f+barf
429       // 68 g+g -> g+g
430       if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
431       if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
432       
433       // fetch the pythia generated jets only to be used here
434       Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
435       AliAODJet pythiaGenJets[kMaxJets];
436       for(int ip = 0;ip < nPythiaGenJets;++ip){
437         if(iCount>=kMaxJets)continue;
438         Float_t p[4];
439         pythiaGenHeader->TriggerJet(ip,p);
440         pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
441         
442         if(fLimitGenJetEta){
443           if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
444              pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
445         }
446         
447
448         if(fBranchGen.Length()==0){
449           // if we have MC particles and we do not read from the aod branch
450           // use the pythia jets
451           genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
452         }
453         iCount++;
454       }
455     }
456     if(fBranchGen.Length()==0)nGenJets = iCount;    
457   }// (fAnalysisType&kMCESD)==kMCESD)
458
459
460   // we simply fetch the tracks/mc particles as a list of AliVParticles
461
462   TList recParticles;
463   TList genParticles;
464
465   GetListOfTracks(&recParticles,fTrackTypeRec);
466   GetListOfTracks(&genParticles,fTrackTypeGen);
467
468   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
469   fh1PtHard->Fill(ptHard,eventW);
470   fh1PtHardNoW->Fill(ptHard,1);
471   fh1PtHardTrials->Fill(ptHard,nTrials);
472
473   // If we set a second branch for the input jets fetch this 
474   if(fBranchGen.Length()>0){
475     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
476     if(aodGenJets){
477       Int_t iCount = 0;
478       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
479         if(iCount>=kMaxJets)continue;
480         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
481         if(!tmp)continue;
482         if(fLimitGenJetEta){
483           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
484              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
485         }
486         genJets[iCount] = *tmp;
487         iCount++;
488       }
489       nGenJets = iCount;
490     }
491     else{
492       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
493       fAOD->Print();
494     }
495   }
496
497   fh1NGenJets->Fill(nGenJets);
498   // We do not want to exceed the maximum jet number
499   nGenJets = TMath::Min(nGenJets,kMaxJets);
500
501   // Fetch the reconstructed jets...
502
503   nRecJets = aodRecJets->GetEntries();
504
505   
506
507
508   fh1NRecJets->Fill(nRecJets);
509   nRecJets = TMath::Min(nRecJets,kMaxJets);
510
511   for(int ir = 0;ir < nRecJets;++ir){
512     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
513     if(!tmp)continue;
514     recJets[ir] = *tmp;
515   }
516
517
518   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
519   // Relate the jets
520   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
521   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
522   
523   for(int i = 0;i<kMaxJets;++i){
524     iGenIndex[i] = iRecIndex[i] = -1;
525   }
526
527
528   AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
529                                             iGenIndex,iRecIndex,fDebug);
530   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
531
532   if(fDebug){
533     for(int i = 0;i<kMaxJets;++i){
534       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
535       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
536     }
537   }
538
539
540
541
542   Double_t container[6];
543
544   // loop over generated jets
545
546   // radius; tmp, get from the jet header itself
547   // at some point, how todeal woht FastJet on MC?
548   Float_t radiusGen =0.4;
549   Float_t radiusRec =0.4;
550
551   for(int ig = 0;ig < nGenJets;++ig){
552     Double_t ptGen = genJets[ig].Pt();
553     Double_t phiGen = genJets[ig].Phi();
554     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
555     Double_t etaGen = genJets[ig].Eta();
556     
557     container[3] = ptGen;
558     container[4] = etaGen;
559     container[5] = phiGen;
560     fhnJetContainer[kStep0]->Fill(&container[3],eventW);
561     Int_t ir = iRecIndex[ig];
562
563     if(TMath::Abs(etaGen)<fRecEtaWindow){
564       fhnJetContainer[kStep1]->Fill(&container[3],eventW);
565       fh1PtGenIn[ig]->Fill(ptGen,eventW);
566       // fill the fragmentation function
567       for(int it = 0;it<genParticles.GetEntries();++it){
568         AliVParticle *part = (AliVParticle*)genParticles.At(it);
569         if(genJets[ig].DeltaR(part)<radiusGen){
570           Float_t z = part->Pt()/ptGen;
571           Float_t lnz =  -1.*TMath::Log(z);
572           fh2FragGen[ig]->Fill(z,ptGen,eventW);
573           fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
574         }
575       }
576       if(ir>=0&&ir<nRecJets){
577         fhnJetContainer[kStep3]->Fill(&container[3],eventW);
578       }
579     }    
580
581     if(ir>=0&&ir<nRecJets){   
582       fhnJetContainer[kStep2]->Fill(&container[3],eventW);
583       Double_t etaRec = recJets[ir].Eta();
584       if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
585     }
586   }// loop over generated jets
587
588   
589   Float_t sumPt = 0;
590   for(int it = 0;it<recParticles.GetEntries();++it){
591     AliVParticle *part = (AliVParticle*)recParticles.At(it);
592     // fill sum pt and P_t of all paritles
593     if(TMath::Abs(part->Eta())<0.9){
594       Float_t pt = part->Pt();
595       fh1PtTrackRec->Fill(pt,eventW);
596       sumPt += pt;
597     }
598   }
599   fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
600   fh1SumPtTrackRec->Fill(sumPt,eventW);
601
602
603   // loop over reconstructed jets
604   for(int ir = 0;ir < nRecJets;++ir){
605     Double_t etaRec = recJets[ir].Eta();
606     Double_t ptRec = recJets[ir].Pt();
607     Double_t phiRec = recJets[ir].Phi();
608     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
609     container[0] = ptRec;
610     container[1] = etaRec;
611     container[2] = phiRec;
612
613     if(ptRec>10.&&fDebug>0){
614       // need to cast to int, otherwise the printf overwrites
615       Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
616       Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
617       fAOD->GetHeader()->Print();
618       Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
619       for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
620         AliAODTrack *tr = fAOD->GetTrack(it);
621         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
622         tr->Print();
623         tr->Dump();
624         if(fESD){
625           AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
626           esdTr->Print("");
627           esdTr->Dump();
628         }
629       }
630     }
631   
632
633     fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
634     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
635     if(TMath::Abs(etaRec)<fRecEtaWindow){
636       fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
637       fh1PtRecIn[ir]->Fill(ptRec,eventW);
638       // fill the fragmentation function
639       for(int it = 0;it<recParticles.GetEntries();++it){
640         AliVParticle *part = (AliVParticle*)recParticles.At(it);
641         Float_t eta = part->Eta();
642         if(TMath::Abs(eta)<0.9){
643           Float_t phi = part->Phi();
644           if(phi<0)phi+=TMath::Pi()*2.;    
645           Float_t dPhi = phi - phiRec;
646           Float_t dEta = eta - etaRec;
647           if(dPhi<(-1.*TMath::Pi()/2))phiRec+=TMath::Pi()*2.;    
648           fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
649           fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
650         }
651         if(recJets[ir].DeltaR(part)<radiusRec){
652           Float_t z = part->Pt()/ptRec;
653           Float_t lnz =  -1.*TMath::Log(z);
654           fh2FragRec[ir]->Fill(z,ptRec,eventW);
655           fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
656         }
657       }
658     }
659
660
661     // Fill Correlation
662     Int_t ig = iGenIndex[ir];
663     if(ig>=0 && ig<nGenJets){
664       fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
665       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
666       Double_t ptGen  = genJets[ig].Pt();
667       Double_t phiGen = genJets[ig].Phi();
668       if(phiGen<0)phiGen+=TMath::Pi()*2.; 
669       Double_t etaGen = genJets[ig].Eta();
670
671       container[3] = ptGen;
672       container[4] = etaGen;
673       container[5] = phiGen;
674
675       // 
676       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
677       // 
678
679       if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
680       if(TMath::Abs(etaRec)<fRecEtaWindow){
681         fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
682         fhnCorrelation->Fill(container);
683       }// if etarec in window
684
685     } 
686     ////////////////////////////////////////////////////  
687     else{
688
689     }
690   }// loop over reconstructed jets
691
692
693   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
694   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
695   PostData(1, fHistList);
696 }
697
698
699 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
700   //
701   // Create the particle container for the correction framework manager and 
702   // link it
703   //
704   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
705   const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
706   const Double_t kEtamin = -3.0, kEtamax = 3.0;
707   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
708
709   // can we neglect migration in eta and phi?
710   // phi should be no problem since we cover full phi and are phi symmetric
711   // eta migration is more difficult  due to needed acceptance correction
712   // in limited eta range
713
714
715   //arrays for the number of bins in each dimension
716   Int_t iBin[kNvar];
717   iBin[0] = 100; //bins in pt
718   iBin[1] =  1; //bins in eta 
719   iBin[2] = 1; // bins in phi
720
721   //arrays for lower bounds :
722   Double_t* binEdges[kNvar];
723   for(Int_t ivar = 0; ivar < kNvar; ivar++)
724     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
725
726   //values for bin lower bounds
727   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);  
728   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
729   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
730   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
731
732
733   for(int i = 0;i < kMaxStep*2;++i){
734     fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
735     for (int k=0; k<kNvar; k++) {
736       fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
737     }
738   }
739   //create correlation matrix for unfolding
740   Int_t thnDim[2*kNvar];
741   for (int k=0; k<kNvar; k++) {
742     //first half  : reconstructed 
743     //second half : MC
744     thnDim[k]      = iBin[k];
745     thnDim[k+kNvar] = iBin[k];
746   }
747
748   fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
749   for (int k=0; k<kNvar; k++) {
750     fhnCorrelation->SetBinEdges(k,binEdges[k]);
751     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
752   }
753   fhnCorrelation->Sumw2();
754
755   // Add a histogram for Fake jets
756   //  thnDim[3] = AliPID::kSPECIES;
757   //  fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
758   // for(Int_t idim = 0; idim < kNvar; idim++)
759   //  fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
760 }
761
762 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
763 {
764 // Terminate analysis
765 //
766     if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
767 }
768
769
770 Int_t  AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
771
772
773   Int_t iCount = 0;
774   if(type==kTrackAODIn||type==kTrackAODOut){
775     AliAODEvent *aod = 0;
776     if(type==kTrackAODIn)aod = dynamic_cast<AliAODEvent*>(InputEvent());
777     else if (type==kTrackAODOut) aod = AODEvent();
778     if(!aod){
779       return iCount;
780     }
781     for(int it = 0;it < aod->GetNumberOfTracks();++it){
782       AliAODTrack *tr = aod->GetTrack(it);
783       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
784       list->Add(tr);
785       iCount++;
786     }
787   }
788   else if (type ==  kTrackKineAll||type == kTrackKineCharged){
789     AliMCEvent* mcEvent = MCEvent();
790     if(!mcEvent)return iCount;
791     // we want to have alivpartilces so use get track
792     for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
793       if(!mcEvent->IsPhysicalPrimary(it))continue;
794       AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
795       if(type == kTrackKineAll){
796         list->Add(part);
797         iCount++;
798       }
799       else if(type == kTrackKineCharged){
800         if(part->Particle()->GetPDG()->Charge()==0)continue;
801         list->Add(part);
802         iCount++;
803       }
804     }
805   }
806   else if (type == kTrackAODMCCharged || type == kTrackAODMCAll ) {
807     AliAODEvent *aod = dynamic_cast<AliAODEvent*>(InputEvent());
808     if(!aod) aod = AODEvent();
809     if(!aod)return iCount;
810     TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
811     if(!tca)return iCount;
812     for(int it = 0;it < tca->GetEntriesFast();++it){
813       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
814       if(!part->IsPhysicalPrimary())continue;
815       if(type == kTrackAODMCAll){
816         list->Add(part);
817         iCount++;
818       }
819       else if (type == kTrackAODMCCharged){
820         if(part->Charge()==0)continue;
821         list->Add(part);
822         iCount++;
823       }
824     }
825   }// AODMCparticle
826   return iCount;
827
828 }