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