Compatibility with the Root trunk
[u/mrichter/AliRoot.git] / PWGJE / UserTasks / AliAnalysisTaskCheckSingleTrackJetRejection.cxx
1 //#include <string.h>
2 //#include <TStyle.h>
3 #include <list>
4 #include <string>
5
6 #include "TTree.h"
7 #include "TCanvas.h"
8 #include "AliAnalysisTask.h"
9 #include "AliInputEventHandler.h"
10 #include "AliESDtrack.h"
11 #include "AliAODVertex.h"
12 #include "AliAODCluster.h"
13
14 #include <TROOT.h>
15 #include <TRandom.h>
16 #include <TSystem.h>
17 #include <TInterpreter.h>
18 #include <TChain.h>
19 #include <TFile.h>
20 #include <TKey.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TProfile.h>
25 #include <TList.h>
26 #include <TLorentzVector.h>
27 #include <TClonesArray.h>
28 #include <TRefArray.h>
29
30 #include "TDatabasePDG.h"
31 #include "AliAnalysisManager.h"
32 #include "AliJetFinder.h"
33 #include "AliJetHeader.h"
34 #include "AliJetReader.h"
35 #include "AliJetReaderHeader.h"
36 #include "AliUA1JetHeaderV1.h"
37 #include "AliSISConeJetHeader.h"
38 #include "AliESDEvent.h"
39 #include "AliAODEvent.h"
40 #include "AliAODHandler.h"
41 #include "AliAODInputHandler.h"
42 #include "AliAODTrack.h"
43 #include "AliVParticle.h"
44 #include "AliAODJet.h"
45 #include "AliAODJetEventBackground.h"
46 #include "AliMCParticle.h"
47 #include "AliAODMCParticle.h"
48 #include "AliMCEventHandler.h"
49 #include "AliMCEvent.h"
50 #include "AliStack.h"
51
52 #include "AliAODHeader.h"
53 #include "AliAODMCHeader.h"
54 #include "AliJetKineReaderHeader.h"
55 #include "AliGenCocktailEventHeader.h"
56 #include "AliInputEventHandler.h"
57 #include "AliGenEventHeader.h"
58 #include "AliGenDPMjetEventHeader.h"
59
60
61
62 #include "AliAnalysisTaskCheckSingleTrackJetRejection.h"
63 #include "AliAnalysisTaskPhiCorrelations.h"
64 #include "AliAnalysisHelperJetTasks.h"
65 #include "AliPWG4HighPtQAMC.h"
66
67
68
69
70 using std::cout;
71 using std::endl;
72
73 ClassImp(AliAnalysisTaskCheckSingleTrackJetRejection)
74
75                                 //________________________________________________________________________
76                                 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection(): 
77                                                                 AliAnalysisTaskSE(),
78                                                                 fUseAODInput(kFALSE),
79                                                                 fFillAOD(kFALSE),
80                                                                 fNonStdFile(""),
81                                                                 fJetBranch("jets"),
82                                                                 fAODIn(0x0),
83                                                                 fAODOut(0x0),
84                                                                 fAODExtension(0x0),
85                                                                 JFAlg("ANTIKT"),
86                                                                 Radius(0.4),
87                                                                 Filtermask(256),
88                                                                 BackM(0),
89                                                                 TrackPtcut(0.15),
90                                                                 SkipCone(0),
91                                                                 IsMC(kTRUE),
92                                                                 fHistList(0x0),
93                                                                 fxsec(0.),
94                                                                 ftrial(1.),
95                                                                 fJetRecEtaWindow(0.5),
96                                                                 fMinJetPt(0.15),
97                                                                 fH1Xsec(0x0),
98                                                                 fH1Trials(0x0),
99                                                                 fH1Events(0x0),
100
101                                                                 fH2jetMCAKT04_Jetpt_maxpt(0x0),
102                                                                 fH2jetAKT04_Jetpt_maxpt  (0x0)
103
104 {
105                                 for(int j=0;j<6;j++){
106                                                                 fH1jetMCAKT04_pt         [j]=0;
107                                                                 fH1jetAKT04_pt           [j]=0;
108
109                                                                 fH2jetMCAKT04_Eratio       [j]=0;
110                                                                 fH1jetMCAKT04_match        [j]=0;
111                                 }
112
113                                 // Default constructor
114 }
115
116 //________________________________________________________________________
117 AliAnalysisTaskCheckSingleTrackJetRejection::AliAnalysisTaskCheckSingleTrackJetRejection(const char *name): 
118                                 AliAnalysisTaskSE(name),
119                                 fUseAODInput(kFALSE),
120                                 fFillAOD(kFALSE),
121                                 fNonStdFile(""),
122                                 fJetBranch("jets"),
123                                 fAODIn(0x0), 
124                                 fAODOut(0x0), 
125                                 fAODExtension(0x0),
126                                 JFAlg("ANTIKT"),
127                                 Radius(0.4),
128                                 Filtermask(256),
129                                 BackM(0),
130                                 TrackPtcut(0.15),
131                                 SkipCone(0),
132                                 IsMC(kTRUE),
133                                 fHistList(0x0),
134                                 fxsec(0.),
135                                 ftrial(1.),
136                                 fJetRecEtaWindow(0.5),
137                                 fMinJetPt(0.15),
138                                 fH1Xsec(0x0),
139                                 fH1Trials(0x0),
140                                 fH1Events(0x0),
141                                 fH2jetMCAKT04_Jetpt_maxpt(0x0),
142                                 fH2jetAKT04_Jetpt_maxpt  (0x0)
143
144 {
145
146                                 for(int j=0;j<6;j++){
147                                                                 fH1jetMCAKT04_pt         [j]=0;
148                                                                 fH1jetAKT04_pt           [j]=0;
149
150                                                                 fH2jetMCAKT04_Eratio       [j]=0;
151                                                                 fH1jetMCAKT04_match        [j]=0;
152                                 }
153
154                                 // Constructor
155                                 // Define input and output slots here
156                                 // Input slot #0 works with a TChain
157                                 DefineInput(0, TChain::Class());
158                                 // Output slot #0 id reserved by the base class for AOD
159                                 // Output slot #1 writes into a TH1 container
160                                 DefineOutput(1, TList::Class());
161
162 }
163
164 //________________________________________________________________________
165 void AliAnalysisTaskCheckSingleTrackJetRejection::UserCreateOutputObjects()
166 {
167                                 // Create histograms
168                                 // Called once
169
170                                 if (!fHistList){ fHistList = new TList();fHistList->SetOwner(kTRUE); cout<<"TList is created for output "<<endl;}
171
172                                 Bool_t oldStatus = TH1::AddDirectoryStatus();
173                                 TH1::AddDirectory(kFALSE);
174
175                                 char *histname;
176                                 if(IsMC){
177                                                                 fH1Xsec           = new TProfile("Xsec","Xsec",1,0,1);
178                                                                 fH1Trials         = new TH1F    ("Trials","Trials",1,0,1);
179                                                                 histname = Form("fH2jetMCAKT04_Jetpt_maxpt");
180                                                                 fH2jetMCAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
181                                                                 histname = Form("fH2jetAKT04_Jetpt_maxpt");
182                                                                 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
183                                                                 for(int j=0;j<6;j++){
184                                                                                                 histname = Form("fH1jetMCAKT04_pt%d",j);
185                                                                                                 fH1jetMCAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
186                                                                                                 histname = Form("fH2jetAKT04_pt%d",j);
187                                                                                                 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
188
189                                                                                                 histname = Form("fH2jetMCAKT04_Eratio%d",j);
190                                                                                                 fH2jetMCAKT04_Eratio[j] = new TH2F(histname,histname,200,0,200,100,0,30);
191                                                                                                 histname = Form("fH1jetMCAKT04_match%d",j);
192                                                                                                 fH1jetMCAKT04_match[j] = new TH1F(histname,histname,200,0,200);
193                                                                 }
194                                                                 fHistList->Add(fH1Xsec);
195                                                                 fHistList->Add(fH1Trials);
196                                                                 fHistList->Add(fH2jetMCAKT04_Jetpt_maxpt);
197                                                                 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
198                                                                 for(int j=0;j<6;j++){
199                                                                                                 fHistList->Add(fH1jetAKT04_pt[j]);
200                                                                                                 fHistList->Add(fH1jetMCAKT04_pt[j]);
201                                                                                                 fHistList->Add(fH2jetMCAKT04_Eratio[j]);
202                                                                                                 fHistList->Add(fH1jetMCAKT04_match[j]);
203                                                                 }
204                                 }
205                                 else    {
206                                                                 fH1Events         = new TH1F("Events","Number of Events",1,0,1);
207                                                                 histname = Form("fH2jetAKT04_Jetpt_maxpt");
208                                                                 fH2jetAKT04_Jetpt_maxpt = new TH2F(histname,histname,400,0,400,400,0,400);
209                                                                 for(int j=0;j<6;j++){
210                                                                                                 histname = Form("fH2jetAKT04_pt%d",j);
211                                                                                                 fH1jetAKT04_pt[j] = new TH1F(histname,histname,200,0,200);
212                                                                 }
213                                                                 fHistList->Add(fH1Events);
214                                                                 fHistList->Add(fH2jetAKT04_Jetpt_maxpt);
215                                                                 for(int j=0;j<6;j++){
216                                                                                                 fHistList->Add(fH1jetAKT04_pt[j]);
217                                                                 }
218                                 }
219
220
221                                 // =========== Switch on Sumw2 for all histos ===========
222                                 for (Int_t i=0; i<fHistList->GetEntries(); ++i) 
223                                 {
224                                                                 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
225                                                                 if (h1)
226                                                                 {
227                                                                                                 h1->Sumw2();
228                                                                                                 continue;
229                                                                 }
230                                                                 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
231                                                                 if(hn)hn->Sumw2();
232                                 }
233                                 TH1::AddDirectory(oldStatus);
234
235
236                                 PostData(1,fHistList);
237
238 }
239
240 //----------------------------------------------------------------------                                                 
241 void AliAnalysisTaskCheckSingleTrackJetRejection::Init()
242 {
243                                 // Initialization                                                                                                    
244                                 if (fDebug) printf("AnalysisTaskCheckSingleTrackJetRejection::Init() \n");
245
246 }
247
248 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::Notify()
249 {
250                                 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
251                                 fAODOut = AODEvent();
252                                 if(fNonStdFile.Length()!=0){
253                                                                 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
254                                                                 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
255                                                                 if(fAODExtension){
256                                                                                                 if(fDebug>1)Printf("AODExtension found for %s ",fNonStdFile.Data());
257                                                                 }
258                                 }
259
260                                 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
261                                 fxsec=0.;
262                                 ftrial=1.;
263
264                                 if(tree){
265                                                                 TFile *curfile = tree->GetCurrentFile();
266                                                                 if(!curfile){
267                                                                                                 Error("Notify","No current file");
268                                                                                                 return kFALSE;
269                                                                 }
270
271                                                                 if(IsMC){
272                                                                                                 PythiaInfoFromFile(curfile->GetName(),fxsec,ftrial);
273                                                                                                 //cout<<" Xsec "<<fxsec<<" trial "<<ftrial<<endl;
274                                                                                                 fH1Xsec  ->Fill(0.,fxsec);
275                                                                                                 fH1Trials->Fill(0.,ftrial);
276                                                                 }
277                                                                 else{
278                                                                                                 //Float_t totalEvent;
279                                                                                                 //totalEvent = GetTotalEvents(curfile->GetName());
280                                                                                                 //fH1Events->Fill(0.,totalEvent);
281                                                                 }
282
283                                 }
284
285                                 printf("Reading File %s ",fInputHandler->GetTree()->GetCurrentFile()->GetName());
286                                 return kTRUE;
287 }
288
289 void AliAnalysisTaskCheckSingleTrackJetRejection::FinishTaskOutput()
290 {
291 }
292
293
294
295 //________________________________________________________________________
296 void AliAnalysisTaskCheckSingleTrackJetRejection::UserExec(Option_t *) 
297 {
298
299
300                                 // Main loop (called each event)
301                                 // Execute analysis for current event
302
303                                 AliAODEvent *fAOD;
304                                 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
305                                 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
306                                                                 fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
307                                                                 if(fUseAODInput) fAODIn = fAOD;
308                                                                 if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
309                                 }
310                                 else {
311                                                                 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
312                                                                 if( handler && handler->InheritsFrom("AliAODHandler") ) {
313                                                                                                 fAOD = ((AliAODHandler*)handler)->GetAOD();
314                                                                                                 fAODIn = fAOD;
315                                                                                                 if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
316                                                                 }
317                                 }
318
319                                 if(!fAODIn && !fUseAODInput){ // case we have AOD in input & output and want jets from output
320                                                                 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
321                                                                 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
322                                                                                                 fAODIn = ((AliAODHandler*)outHandler)->GetAOD();
323                                                                                                 if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
324                                                                 }
325                                 }
326
327                                 //if(fNonStdFile.Length()!=0){
328                                 //                              // case we have an AOD extension - fetch the jets from the extended output
329                                 //                              AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
330                                 //                              fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
331                                 //                              if(!fAODExtension){
332                                 //                                                              if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
333                                 //                              }
334                                 //}
335                                 //fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
336                                 if (!fAODIn) {
337                                                                 Printf("ERROR %s : fAODIn not available",(char*)__FILE__);
338                                                                 return;
339                                 }
340
341                                 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
342                                                                 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
343                                 if(!(inputHandler->IsEventSelected() & AliVEvent::kMB)){
344                                                                 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
345                                                                 return;
346                                 }
347                                 if(!IsMC)fH1Events->Fill(0);
348
349
350
351                                 // start jet analysis
352
353                                 Double_t Jet_n  [20];
354                                 Double_t Jet_pt [20][10000];
355                                 Double_t Jet_eta[20][10000];
356                                 Double_t Jet_phi[20][10000];
357
358                                 for(int i=0;i<20;i++){
359                                                                 Jet_n[i]=0;
360                                                                 for(int j=0;j<10000;j++){
361                                                                                                 Jet_pt[i][j]=0.;
362                                                                                                 Jet_phi[i][j]=999.;
363                                                                                                 Jet_eta[i][j]=999.;
364                                                                 }
365                                 }
366
367                                 // get AOD event from input/ouput
368                                 TString cAdd = "";
369                                 cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
370                                 cAdd += Form("B%d",(int)BackM);
371                                 cAdd += Form("_Filter%05d",Filtermask);
372                                 cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
373                                 cAdd += Form("_Skip%02d",SkipCone);
374                                 TString Branchname_gen,Branchname_rec;
375                                 //Branchname_gen = Form("clustersMCKINE2_%s%s",JFAlg.Data(),cAdd.Data());
376                                 Branchname_gen = Form("clustersAODMC2_%s%s",JFAlg.Data(),cAdd.Data());
377                                 Branchname_rec = Form("clustersAOD_%s%s",JFAlg.Data(),cAdd.Data());
378
379
380                                 bool fFIND[6][1000];
381                                 double maxpt[2][1000];for(int i=0;i<1000;i++){maxpt[0][i]=0;maxpt[1][i]=0;}
382                                 int nearrecID[1000];  for(int i=0;i<1000;i++){nearrecID[i]=99999;}
383                                 AliAODJet* jetsAOD;
384                                 for(int algorithm=0;algorithm<2;algorithm++){
385                                                                 //for LHC11a1  LHC11a2 official
386                                                                 if((!IsMC&&algorithm==0))continue;
387
388                                                                 if(algorithm==0)fJetBranch   = Branchname_gen.Data();
389                                                                 if(algorithm==1)fJetBranch   = Branchname_rec.Data();
390
391                                                                 TClonesArray* jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
392                                                                 if(!jets){
393                                                                                                 printf(" Tere are no Branch named %s \n",fJetBranch.Data());
394                                                                                                 continue;
395                                                                 }
396                                                                 Int_t nj = jets->GetEntriesFast();
397                                                                 if (fDebug) printf("There are %5d jets in the event \n", nj);
398
399
400                                                                 Jet_n[algorithm] = nj;
401                                                                 for(int njet =0;njet<nj;njet++){
402                                                                                                 jetsAOD = (AliAODJet*) (jets->At(njet));
403                                                                                                 Jet_pt   [algorithm][njet] = jetsAOD->Pt();
404                                                                                                 Jet_phi  [algorithm][njet] = jetsAOD->Phi();  
405                                                                                                 Jet_eta  [algorithm][njet] = jetsAOD->Eta();
406                                                                                                 double eta_cut_Jet=0.5;
407                                                                                                 TRefArray *reftracks = jetsAOD->GetRefTracks();
408                                                                                                 int ntracks=reftracks->GetEntriesFast();
409                                                                                                 if(TMath::Abs(Jet_eta[algorithm][njet])<eta_cut_Jet){
410                                                                                                                                 //------------calc max pt in Jet----------------
411                                                                                                                                 double trackpt;
412                                                                                                                                 double sumtrackpt=0;//test
413                                                                                                                                 for(int ntr=0;ntr<ntracks;ntr++){// calc. max pt of track which is in Jet
414                                                                                                                                                                 AliAODTrack *AODtrack = dynamic_cast<AliAODTrack*>(reftracks->At(ntr));
415                                                                                                                                                                 if(AODtrack){
416                                                                                                                                                                                                 Bool_t bgoodT=false;
417                                                                                                                                                                                                 if(Filtermask!=768){
418                                                                                                                                                                                                                                 if(AODtrack->TestFilterMask(Filtermask))bgoodT=true;
419                                                                                                                                                                                                 }
420                                                                                                                                                                                                 else{
421                                                                                                                                                                                                                                 if(AODtrack->IsHybridGlobalConstrainedGlobal())bgoodT=true; //for hybrid Track cuts
422                                                                                                                                                                                                 }
423                                                                                                                                                                                                 if(!bgoodT)continue;
424                                                                                                                                                                                                 trackpt = AODtrack->Pt();
425                                                                                                                                                                                                 sumtrackpt += trackpt;
426                                                                                                                                                                                                 if(trackpt>maxpt[algorithm][njet]){
427                                                                                                                                                                                                                                 maxpt[algorithm][njet] = trackpt;
428                                                                                                                                                                                                 }
429                                                                                                                                                                 }
430                                                                                                                                 }// track Loop
431
432                                                                                                                                 if(algorithm==0){
433                                                                                                                                                                 fH1jetMCAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
434                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetMCAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
435                                                                                                                                                                 if(ntracks>1)     fH1jetMCAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
436                                                                                                                                                                 if(ntracks>2)     fH1jetMCAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
437                                                                                                                                                                 if(ntracks>3)     fH1jetMCAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
438                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetMCAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
439                                                                                                                                                                 fH2jetMCAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
440                                                                                                                                 }
441                                                                                                                                 if(algorithm==1){
442                                                                                                                                                                 fH1jetAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
443                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
444                                                                                                                                                                 if(ntracks>1)     fH1jetAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
445                                                                                                                                                                 if(ntracks>2)     fH1jetAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
446                                                                                                                                                                 if(ntracks>3)     fH1jetAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
447                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
448                                                                                                                                                                 fH2jetAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
449                                                                                                                                 }
450                                                                                                 }
451                                                                 }
452                                                                 if(!(IsMC&&algorithm==1))continue;
453                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
454                                                                                                 double eta_cut_Jet=0.5;
455                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
456                                                                                                                                 //Find muched jet pare=====================================
457                                                                                                                                 for(int cut=0;cut<6;cut++){
458                                                                                                                                                                 double min_R=10.;
459                                                                                                                                                                 for(int njetAOD=0;njetAOD<Jet_n[1];njetAOD++){
460                                                                                                                                                                   jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
461                                                                                                                                                                   if(!jets)continue;
462                                                                                                                                                                   jetsAOD = (AliAODJet*) (jets->At(njetAOD));
463                                                                                                                                                                   TRefArray *reftracks = jetsAOD->GetRefTracks();
464                                                                                                                                                                   int ntracks=reftracks->GetEntriesFast();
465                                                                                                                                                                                                 if(cut==1){if(maxpt[1][njetAOD]>=40.)continue;}
466                                                                                                                                                                                                 if(cut==2){if(ntracks==1)continue;}
467                                                                                                                                                                                                 if(cut==3){if(ntracks<=2)continue;}
468                                                                                                                                                                                                 if(cut==4){if(ntracks<=3)continue;}
469                                                                                                                                                                                                 if(cut==5){if(maxpt[1][njetAOD]>=40.)continue;if(ntracks==1)continue;}
470                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njetMC],Jet_phi[1][njetAOD],Jet_eta[0][njetMC],Jet_eta[1][njetAOD]);
471                                                                                                                                                                                                 if(DelR<min_R){
472                                                                                                                                                                                                                                 nearrecID[njetMC]=njetAOD;
473                                                                                                                                                                                                                                 min_R=DelR;
474                                                                                                                                                                                                 }
475                                                                                                                                                                 }
476                                                                                                                                                                 if(min_R<0.4){
477                                                                                                                                                                                                 min_R=10.;
478                                                                                                                                                                                                 int neargenID=99999;
479                                                                                                                                                                                                 for(int njet =0;njet<Jet_n[0];njet++){
480                                                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njet],Jet_phi[1][nearrecID[njetMC]],Jet_eta[0][njet],Jet_eta[1][nearrecID[njetMC]]);
481                                                                                                                                                                                                                                 if(DelR<min_R){
482                                                                                                                                                                                                                                                                 neargenID=njet;
483                                                                                                                                                                                                                                                                 min_R=DelR;
484                                                                                                                                                                                                                                 }
485                                                                                                                                                                                                 }
486                                                                                                                                                                                                 if((min_R<0.4)&&(neargenID==njetMC))fFIND[cut][njetMC]=true;
487                                                                                                                                                                 }
488                                                                                                                                 }
489                                                                                                                                 //======================================================
490                                                                                                 }
491                                                                 }
492                                 }//algorithm
493                                 if(IsMC){
494                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
495                                                                                                 double eta_cut_Jet=0.5;
496                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
497                                                                                                                                 for(int cut=0;cut<6;cut++){
498                                                                                                                                                                 if(fFIND[cut][njetMC]==true){
499                                                                                                                                                                                                 fH1jetMCAKT04_match[cut]->Fill(Jet_pt[0][njetMC]);
500                                                                                                                                                                                                 fH2jetMCAKT04_Eratio[cut]->Fill(Jet_pt[0][njetMC],Jet_pt[1][nearrecID[njetMC]]/Jet_pt[0][njetMC]);
501                                                                                                                                                                 }
502                                                                                                                                 }
503                                                                                                 }//eta window
504                                                                 }//njet loop
505                                 }
506
507                                 // End of jet analysis ---------
508                                 // Post output data.
509                                 PostData(1, fHistList);
510                                 return;
511 }      
512
513 //________________________________________________________________________
514 void AliAnalysisTaskCheckSingleTrackJetRejection::Terminate(Option_t *) 
515 {
516                                 // Terminate analysis
517                                 //
518                                 if (fDebug) printf("AnalysisTaskPt: Terminate() \n");
519
520 }
521
522
523 Bool_t  AliAnalysisTaskCheckSingleTrackJetRejection::JetSelected(AliAODJet *jet){
524                                 Bool_t selected = false;
525
526                                 if(!jet)return selected;
527
528                                 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
529                                                                 selected = kTRUE;
530                                 }
531                                 return selected;
532
533 }
534
535
536 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaPhi(Double_t phi1,Double_t phi2){
537                                 Float_t pi=TMath::Pi();
538                                 Double_t dphi = phi1-phi2;
539                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
540                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
541                                 return dphi;
542 }
543 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaR(Double_t phi1,Double_t phi2,Double_t eta1,Double_t eta2){
544                                 Float_t pi=TMath::Pi();
545                                 Double_t dphi = DeltaPhi(phi1,phi2);
546                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
547                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
548                                 Double_t Deta = eta1 - eta2;
549                                 Double_t deltaR = sqrt(pow(dphi,2)+pow(Deta,2));
550                                 return deltaR;
551 }
552
553 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
554                                 //
555                                 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
556                                 // This is to called in Notify and should provide the path to the AOD/ESD file
557                                 // Copied from AliAnalysisTaskJetSpectrum2
558                                 //
559
560                                 TString file(currFile);
561                                 fXsec = 0;
562                                 fTrials = 1;
563
564                                 if(file.Contains("root_archive.zip#")){
565                                                                 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
566                                                                 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
567                                                                 file.Replace(pos+1,20,"");
568                                 }
569                                 else {
570                                                                 // not an archive take the basename....
571                                                                 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
572                                 }
573                                 //  Printf("%s",file.Data());
574
575
576                                 TFile *filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really 
577                                 //test the existance of a file in a archive so we have to lvie with open error message from root
578                                                                 if(!filexsec){
579                                                                                                 // next trial fetch the histgram file
580                                                                                                 filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
581                                                                                                 if(!filexsec){
582                                                                                                                                 // not a severe condition but inciate that we have no information
583                                                                                                                                 return kFALSE;
584                                                                                                 }
585                                                                                                 else{
586                                                                                                                                 // find the tlist we want to be independtent of the name so use the Tkey
587                                                                                                                                 TKey* key = (TKey*)filexsec->GetListOfKeys()->At(0);
588                                                                                                                                 if(!key){
589                                                                                                                                                                 filexsec->Close();
590                                                                                                                                                                 return kFALSE;
591                                                                                                                                 }
592                                                                                                                                 TList *list = dynamic_cast<TList*>(key->ReadObj());
593                                                                                                                                 if(!list){
594                                                                                                                                                                 filexsec->Close();
595                                                                                                                                                                 return kFALSE;
596                                                                                                                                 }
597                                                                                                                                 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
598                                                                                                                                 fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
599                                                                                                                                 filexsec->Close();
600                                                                                                 }
601                                                                 } // no tree pyxsec.root
602                                                                 else {
603                                                                                                 TTree *xtree = (TTree*)filexsec->Get("Xsection");
604                                                                                                 if(!xtree){
605                                                                                                                                 filexsec->Close();
606                                                                                                                                 return kFALSE;
607                                                                                                 }
608                                                                                                 UInt_t   ntrials  = 0;
609                                                                                                 Double_t  xsection  = 0;
610                                                                                                 xtree->SetBranchAddress("xsection",&xsection);
611                                                                                                 xtree->SetBranchAddress("ntrials",&ntrials);
612                                                                                                 xtree->GetEntry(0);
613                                                                                                 fTrials = ntrials;
614                                                                                                 fXsec = xsection;
615                                                                                                 filexsec->Close();
616                                                                 }
617                                 return kTRUE;
618 }
619
620
621 //Float_t AliAnalysisTaskCheckSingleTrackJetRejection::GetTotalEvents(const char* currFile){
622 //                              Float_t totalevent;
623 //                              TString file_es(currFile);
624 //                              if(file_es.Contains("root_archive.zip#")){
625 //                                                              Ssiz_t pos1 = file_es.Index("root_archive",12,TString::kExact);
626 //                                                              Ssiz_t pos = file_es.Index("#",1,pos1,TString::kExact);
627 //                                                              file_es.Replace(pos+1,20,"");
628 //                              }
629 //                              else {
630 //                                                              // not an archive take the basename....
631 //                                                              file_es.ReplaceAll(gSystem->BaseName(file_es.Data()),"");
632 //                              }
633 //
634 //                              TString cAdd = "";
635 //                              cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
636 //                              cAdd += Form("B%d",(int)BackM);
637 //                              cAdd += Form("_Filter%05d",Filtermask);
638 //                              cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
639 //                              cAdd += Form("_Skip%02d",SkipCone);
640 //                              TString Dirname,Listname;
641 //                              Dirname  = Form("PWG4_cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
642 //                              Listname = Form("pwg4cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
643 //
644 //                              TFile *feventstat = TFile::Open(Form("%s%s",file_es.Data(),"JetTasksOutput.root"));
645 //                              gROOT->Cd(Dirname.Data());
646 //                              TList *templist     = (TList*)gROOT->FindObject(Listname.Data());
647 //                              TH1F* temphist = (TH1F*)templist->FindObject("fh1Trials");
648 //                              totalevent = temphist->Integral();
649 //                              //cout<<temphist->Integral()<<endl;
650 //                              delete feventstat;
651 //                              return totalevent;
652 //
653 //}
654
655