Coverity fix (forward Null)
[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                                                                 if(!jets)continue;
345                                                                 Int_t nj = jets->GetEntriesFast();
346                                                                 if (fDebug) printf("There are %5d jets in the event \n", nj);
347
348
349                                                                 Jet_n[algorithm] = nj;
350                                                                 for(int njet =0;njet<nj;njet++){
351                                                                                                 jetsAOD = (AliAODJet*) (jets->At(njet));
352                                                                                                 Jet_pt   [algorithm][njet] = jetsAOD->Pt();
353                                                                                                 Jet_phi  [algorithm][njet] = jetsAOD->Phi();  
354                                                                                                 Jet_eta  [algorithm][njet] = jetsAOD->Eta();
355                                                                                                 double eta_cut_Jet=0.5;
356                                                                                                 TRefArray *reftracks = jetsAOD->GetRefTracks();
357                                                                                                 int ntracks=reftracks->GetEntriesFast();
358                                                                                                 if(TMath::Abs(Jet_eta[algorithm][njet])<eta_cut_Jet){
359                                                                                                                                 //------------calc max pt in Jet----------------
360                                                                                                                                 double trackpt;
361                                                                                                                                 double sumtrackpt=0;//test
362                                                                                                                                 for(int ntr=0;ntr<ntracks;ntr++){// calc. max pt of track which is in Jet
363                                                                                                                                                                 AliAODTrack *AODtrack = dynamic_cast<AliAODTrack*>(reftracks->At(ntr));
364                                                                                                                                                                 if(AODtrack){
365                                                                                                                                                                                                 if(AODtrack->TestFilterMask(256)){
366                                                                                                                                                                                                                                 trackpt = AODtrack->Pt();
367                                                                                                                                                                                                                                 sumtrackpt += trackpt;
368                                                                                                                                                                                                                                 if(trackpt>maxpt[algorithm][njet]){
369                                                                                                                                                                                                                                                                 maxpt[algorithm][njet] = trackpt;
370                                                                                                                                                                                                                                 }
371                                                                                                                                                                                                 }
372                                                                                                                                                                 }
373                                                                                                                                 }// track Loop
374
375                                                                                                                                 if(algorithm==0){
376                                                                                                                                                                 fH1jetMCAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
377                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetMCAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
378                                                                                                                                                                 if(ntracks>1)     fH1jetMCAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
379                                                                                                                                                                 if(ntracks>2)     fH1jetMCAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
380                                                                                                                                                                 if(ntracks>3)     fH1jetMCAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
381                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetMCAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
382                                                                                                                                                                 fH2jetMCAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
383                                                                                                                                 }
384                                                                                                                                 if(algorithm==1){
385                                                                                                                                                                 fH1jetAKT04_pt[0]->Fill(Jet_pt[algorithm][njet]);
386                                                                                                                                                                 if(maxpt[algorithm][njet]<40)fH1jetAKT04_pt[1]->Fill(Jet_pt[algorithm][njet]);
387                                                                                                                                                                 if(ntracks>1)     fH1jetAKT04_pt[2]->Fill(Jet_pt[algorithm][njet]);
388                                                                                                                                                                 if(ntracks>2)     fH1jetAKT04_pt[3]->Fill(Jet_pt[algorithm][njet]);
389                                                                                                                                                                 if(ntracks>3)     fH1jetAKT04_pt[4]->Fill(Jet_pt[algorithm][njet]);
390                                                                                                                                                                 if((maxpt[algorithm][njet]<40)&&(ntracks>1))fH1jetAKT04_pt[5]->Fill(Jet_pt[algorithm][njet]);
391                                                                                                                                                                 fH2jetAKT04_Jetpt_maxpt->Fill(maxpt[algorithm][njet],Jet_pt[algorithm][njet]);
392                                                                                                                                 }
393                                                                                                 }
394                                                                 }
395                                                                 if(!(IsMC&&algorithm==1))continue;
396                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
397                                                                                                 double eta_cut_Jet=0.5;
398                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
399                                                                                                                                 //Find muched jet pare=====================================
400                                                                                                                                 for(int cut=0;cut<6;cut++){
401                                                                                                                                                                 double min_R=10.;
402                                                                                                                                                                 for(int njetAOD=0;njetAOD<Jet_n[1];njetAOD++){
403                                                                                                                                                                                                 fJetBranch   = "clustersAOD_ANTIKT04_B0_Filter00256_Cut00150_Skip00";
404                                                                                                                                                                                                 jets = dynamic_cast <TClonesArray*> (fAODIn->FindListObject(fJetBranch.Data()));
405                                                                                                                                                                                                 if(jets)continue;
406                                                                                                                                                                                                 jetsAOD = (AliAODJet*) (jets->At(njetAOD));
407                                                                                                                                                                                                 TRefArray *reftracks = jetsAOD->GetRefTracks();
408                                                                                                                                                                                                 int ntracks=reftracks->GetEntriesFast();
409                                                                                                                                                                                                 if(cut==1){if(maxpt[1][njetAOD]>=40.)continue;}
410                                                                                                                                                                                                 if(cut==2){if(ntracks==1)continue;}
411                                                                                                                                                                                                 if(cut==3){if(ntracks<=2)continue;}
412                                                                                                                                                                                                 if(cut==4){if(ntracks<=3)continue;}
413                                                                                                                                                                                                 if(cut==5){if(maxpt[1][njetAOD]>=40.)continue;if(ntracks==1)continue;}
414                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njetMC],Jet_phi[1][njetAOD],Jet_eta[0][njetMC],Jet_eta[1][njetAOD]);
415                                                                                                                                                                                                 if(DelR<min_R){
416                                                                                                                                                                                                                                 nearrecID[njetMC]=njetAOD;
417                                                                                                                                                                                                                                 min_R=DelR;
418                                                                                                                                                                                                 }
419                                                                                                                                                                 }
420                                                                                                                                                                 if(min_R<0.4){
421                                                                                                                                                                                                 min_R=10.;
422                                                                                                                                                                                                 int neargenID=99999;
423                                                                                                                                                                                                 for(int njet =0;njet<Jet_n[0];njet++){
424                                                                                                                                                                                                                                 double DelR = DeltaR(Jet_phi[0][njet],Jet_phi[1][nearrecID[njetMC]],Jet_eta[0][njet],Jet_eta[1][nearrecID[njetMC]]);
425                                                                                                                                                                                                                                 if(DelR<min_R){
426                                                                                                                                                                                                                                                                 neargenID=njet;
427                                                                                                                                                                                                                                                                 min_R=DelR;
428                                                                                                                                                                                                                                 }
429                                                                                                                                                                                                 }
430                                                                                                                                                                                                 if((min_R<0.4)&&(neargenID==njetMC))fFIND[cut][njetMC]=true;
431                                                                                                                                                                 }
432                                                                                                                                 }
433                                                                                                                                 //======================================================
434                                                                                                 }
435                                                                 }
436                                 }//algorithm
437                                 if(IsMC){
438                                                                 for(int njetMC =0;njetMC<Jet_n[0];njetMC++){
439                                                                                                 double eta_cut_Jet=0.5;
440                                                                                                 if(TMath::Abs(Jet_eta[0][njetMC])<eta_cut_Jet){
441                                                                                                                                 for(int cut=0;cut<6;cut++){
442                                                                                                                                                                 if(fFIND[cut][njetMC]==true){
443                                                                                                                                                                                                 fH1jetMCAKT04_match[cut]->Fill(Jet_pt[0][njetMC]);
444                                                                                                                                                                                                 fH2jetMCAKT04_Eratio[cut]->Fill(Jet_pt[0][njetMC],Jet_pt[1][nearrecID[njetMC]]/Jet_pt[0][njetMC]);
445                                                                                                                                                                 }
446                                                                                                                                 }
447                                                                                                 }//eta window
448                                                                 }//njet loop
449                                 }
450
451                                 // End of jet analysis ---------
452                                 // Post output data.
453                                 PostData(1, fHistList);
454                                 return;
455 }      
456
457 //________________________________________________________________________
458 void AliAnalysisTaskCheckSingleTrackJetRejection::Terminate(Option_t *) 
459 {
460                                 // Terminate analysis
461                                 //
462                                 if (fDebug) printf("AnalysisTaskPt: Terminate() \n");
463
464 }
465
466
467 Bool_t  AliAnalysisTaskCheckSingleTrackJetRejection::JetSelected(AliAODJet *jet){
468                                 Bool_t selected = false;
469
470                                 if(!jet)return selected;
471
472                                 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
473                                                                 selected = kTRUE;
474                                 }
475                                 return selected;
476
477 }
478
479
480 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaPhi(Double_t phi1,Double_t phi2){
481                                 Float_t pi=TMath::Pi();
482                                 Double_t dphi = phi1-phi2;
483                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
484                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
485                                 return dphi;
486 }
487 Double_t AliAnalysisTaskCheckSingleTrackJetRejection::DeltaR(Double_t phi1,Double_t phi2,Double_t eta1,Double_t eta2){
488                                 Float_t pi=TMath::Pi();
489                                 Double_t dphi = DeltaPhi(phi1,phi2);
490                                 if     (dphi<(-1./2*pi))dphi = dphi +2*pi;
491                                 else if(dphi>( 3./2*pi))dphi = dphi -2*pi;
492                                 Double_t Deta = eta1 - eta2;
493                                 Double_t deltaR = sqrt(pow(dphi,2)+pow(Deta,2));
494                                 return deltaR;
495 }
496
497 Bool_t AliAnalysisTaskCheckSingleTrackJetRejection::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
498                                 //
499                                 // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
500                                 // This is to called in Notify and should provide the path to the AOD/ESD file
501                                 // Copied from AliAnalysisTaskJetSpectrum2
502                                 //
503
504                                 TString file(currFile);
505                                 fXsec = 0;
506                                 fTrials = 1;
507
508                                 if(file.Contains("root_archive.zip#")){
509                                                                 Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
510                                                                 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
511                                                                 file.Replace(pos+1,20,"");
512                                 }
513                                 else {
514                                                                 // not an archive take the basename....
515                                                                 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
516                                 }
517                                 //  Printf("%s",file.Data());
518
519
520                                 TFile *filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really 
521                                 //test the existance of a file in a archive so we have to lvie with open error message from root
522                                                                 if(!filexsec){
523                                                                                                 // next trial fetch the histgram file
524                                                                                                 filexsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
525                                                                                                 if(!filexsec){
526                                                                                                                                 // not a severe condition but inciate that we have no information
527                                                                                                                                 return kFALSE;
528                                                                                                 }
529                                                                                                 else{
530                                                                                                                                 // find the tlist we want to be independtent of the name so use the Tkey
531                                                                                                                                 TKey* key = (TKey*)filexsec->GetListOfKeys()->At(0);
532                                                                                                                                 if(!key){
533                                                                                                                                                                 filexsec->Close();
534                                                                                                                                                                 return kFALSE;
535                                                                                                                                 }
536                                                                                                                                 TList *list = dynamic_cast<TList*>(key->ReadObj());
537                                                                                                                                 if(!list){
538                                                                                                                                                                 filexsec->Close();
539                                                                                                                                                                 return kFALSE;
540                                                                                                                                 }
541                                                                                                                                 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
542                                                                                                                                 fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
543                                                                                                                                 filexsec->Close();
544                                                                                                 }
545                                                                 } // no tree pyxsec.root
546                                                                 else {
547                                                                                                 TTree *xtree = (TTree*)filexsec->Get("Xsection");
548                                                                                                 if(!xtree){
549                                                                                                                                 filexsec->Close();
550                                                                                                                                 return kFALSE;
551                                                                                                 }
552                                                                                                 UInt_t   ntrials  = 0;
553                                                                                                 Double_t  xsection  = 0;
554                                                                                                 xtree->SetBranchAddress("xsection",&xsection);
555                                                                                                 xtree->SetBranchAddress("ntrials",&ntrials);
556                                                                                                 xtree->GetEntry(0);
557                                                                                                 fTrials = ntrials;
558                                                                                                 fXsec = xsection;
559                                                                                                 filexsec->Close();
560                                                                 }
561                                 return kTRUE;
562 }
563
564
565 Float_t AliAnalysisTaskCheckSingleTrackJetRejection::GetTotalEvents(const char* currFile){
566                                 Float_t totalevent;
567                                 TString file_es(currFile);
568                                 if(file_es.Contains("root_archive.zip#")){
569                                                                 Ssiz_t pos1 = file_es.Index("root_archive",12,TString::kExact);
570                                                                 Ssiz_t pos = file_es.Index("#",1,pos1,TString::kExact);
571                                                                 file_es.Replace(pos+1,20,"");
572                                 }
573                                 else {
574                                                                 // not an archive take the basename....
575                                                                 file_es.ReplaceAll(gSystem->BaseName(file_es.Data()),"");
576                                 }
577
578                                 TString cAdd = "";
579                                 cAdd += Form("%02d_",(int)((Radius+0.01)*10.));
580                                 cAdd += Form("B%d",(int)BackM);
581                                 cAdd += Form("_Filter%05d",Filtermask);
582                                 cAdd += Form("_Cut%05d",(int)(1000.*TrackPtcut));
583                                 cAdd += Form("_Skip%02d",SkipCone);
584                                 TString Dirname,Listname;
585                                 Dirname  = Form("PWG4_cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
586                                 Listname = Form("pwg4cluster_AOD__%s%s",JFAlg.Data(),cAdd.Data());
587
588                                 TFile *feventstat = TFile::Open(Form("%s%s",file_es.Data(),"JetTasksOutput.root"));
589                                 //gROOT->Cd("PWG4_cluster_AOD__KT06_B0_Filter00256_Cut00150_Skip00");
590                                 //TList *templist     = (TList*)gROOT->FindObject("pwg4cluster_AOD__KT06_B0_Filter00256_Cut00150_Skip00");
591                                 gROOT->Cd(Dirname.Data());
592                                 TList *templist     = (TList*)gROOT->FindObject(Listname.Data());
593                                 TH1F* temphist = (TH1F*)templist->FindObject("fh1Trials");
594                                 totalevent = temphist->Integral();
595                                 //cout<<temphist->Integral()<<endl;
596                                 delete feventstat;
597                                 return totalevent;
598
599 }
600
601