Set jet area to jets from lightweight jet finder
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskAj.cxx
1
2 // ******************************************
3 // This task searches for events with large dijet imbalance
4 // and then looks to the jet structure of the b-t-b jets.    
5 // *******************************************
6
7
8 /**************************************************************************
9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10  *                                                                        *
11  * Author: The ALICE Off-line Project.                                    *
12  * Contributors are mentioned in the code where appropriate.              *
13  *                                                                        *
14  * Permission to use, copy, modify and distribute this software and its   *
15  * documentation strictly for non-commercial purposes is hereby granted   *
16  * without fee, provided that the above copyright notice appears in all   *
17  * copies and that both the copyright notice and this permission notice   *
18  * appear in the supporting documentation. The authors make no claims     *
19  * about the suitability of this software for any purpose. It is          *
20  * provided "as is" without express or implied warranty.                  *
21  **************************************************************************/
22
23
24 #include "TChain.h"
25 #include "TTree.h"
26 #include "TMath.h"
27 #include "TH1F.h"
28 #include "TH2F.h"
29 #include "TH3F.h"
30 #include "THnSparse.h"
31 #include "TMatrixD.h"
32 #include "TMatrixDSym.h"
33 #include "TMatrixDSymEigen.h"
34
35 #include "TCanvas.h"
36
37 #include "AliLog.h"
38
39 #include "AliAnalysisTask.h"
40 #include "AliAnalysisManager.h"
41
42 #include "AliVEvent.h"
43 #include "AliESDEvent.h"
44 #include "AliESDInputHandler.h"
45 #include "AliCentrality.h"
46 #include "AliAnalysisHelperJetTasks.h"
47 #include "AliInputEventHandler.h"
48 #include "AliAODJetEventBackground.h"
49 #include "AliAnalysisTaskFastEmbedding.h"
50 #include "AliAODEvent.h"
51 #include "AliAODHandler.h"
52 #include "AliAODJet.h"
53
54 #include "AliAnalysisTaskAj.h"
55
56 using std::cout;
57 using std::endl;
58
59 ClassImp(AliAnalysisTaskAj)
60
61 AliAnalysisTaskAj::AliAnalysisTaskAj() :
62 AliAnalysisTaskSE(),
63 fESD(0x0),
64 fAODIn(0x0),
65 fAODOut(0x0),
66 fAODExtension(0x0),
67 fBackgroundBranch(""),
68 fNonStdFile(""),
69 fIsPbPb(kTRUE),
70 fOfflineTrgMask(AliVEvent::kAny),
71 fMinContribVtx(1),
72 fVtxZMin(-10.),
73 fVtxZMax(10.),
74 fEvtClassMin(0),
75 fEvtClassMax(4),
76 fFilterMask(0),
77 fRadioFrac(0.2),
78 fMinDist(0.1),
79 fCentMin(0.),
80 fCentMax(100.),
81 fNInputTracksMin(0),
82 fNInputTracksMax(-1),
83 fAngStructCloseTracks(0),
84 fCheckMethods(0),
85 fDoEventMixing(0), 
86 fJetEtaMin(-.5),
87 fJetEtaMax(.5),
88 fNevents(0x0),
89 fTindex(0x0),
90 fJetPtMin(20.),
91 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
92 fJetPtFractionMin(0.5),
93 fNMatchJets(4),
94 fMatchMaxDist(0.8),
95 fKeepJets(kFALSE),
96 fkNbranches(2),
97 fkEvtClasses(12),
98 fOutputList(0x0),
99 fbEvent(kTRUE),
100 fHistEvtSelection(0x0),
101 fh2Pt1Pt2C10(0x0),
102 fh2Pt1Pt2C40(0x0), 
103 fh3LocalCoordinates(0x0), 
104 fh2Sum2pt20(0x0),
105 fh2Sum4pt20(0x0),
106 fh2Sum6pt20(0x0),
107 fh2Sum8pt20(0x0),
108 fh2Sum12pt20(0x0),
109 fh2Sum2lpt20(0x0),
110 fh2Sum4lpt20(0x0),
111 fh2Sum6lpt20(0x0),
112 fh2Sum8lpt20(0x0),
113 fh2Sum12lpt20(0x0),
114 fh2Sum2pt40(0x0),
115 fh2Sum4pt40(0x0),
116 fh2Sum6pt40(0x0),
117 fh2Sum8pt40(0x0),
118 fh2Sum12pt40(0x0),
119 fh2Sum2lpt40(0x0),
120 fh2Sum4lpt40(0x0),
121 fh2Sum6lpt40(0x0),
122 fh2Sum8lpt40(0x0),
123 fh2Sum12lpt40(0x0),
124 fhnDeltaR(0x0)
125  {
126    // default Constructor
127
128
129    fJetBranchName[0] = "";
130    fJetBranchName[1] = "";
131
132    fListJets[0] = new TList;
133    fListJets[1] = new TList;
134 }
135
136 AliAnalysisTaskAj::AliAnalysisTaskAj(const char *name) :
137 AliAnalysisTaskSE(name),
138 fESD(0x0),
139 fAODIn(0x0),
140 fAODOut(0x0),
141 fAODExtension(0x0),
142 fBackgroundBranch(""),
143 fNonStdFile(""),
144 fIsPbPb(kTRUE),
145 fOfflineTrgMask(AliVEvent::kAny),
146 fMinContribVtx(1),
147 fVtxZMin(-10.),
148 fVtxZMax(10.),
149 fEvtClassMin(0),
150 fEvtClassMax(4),
151 fFilterMask(0),
152 fRadioFrac(0.2),
153 fMinDist(0.1),
154 fCentMin(0.),
155 fCentMax(100.),
156 fNInputTracksMin(0),
157 fNInputTracksMax(-1),
158 fAngStructCloseTracks(0),
159 fCheckMethods(0),
160 fDoEventMixing(0),
161 fJetEtaMin(-.5),
162 fJetEtaMax(.5),
163 fNevents(0x0),
164 fTindex(0x0),
165 fJetPtMin(20.),
166 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
167 fJetPtFractionMin(0.5),
168 fNMatchJets(4),
169 fMatchMaxDist(0.8),
170 fKeepJets(kFALSE),
171 fkNbranches(2),
172 fkEvtClasses(12),
173 fOutputList(0x0),
174 fbEvent(kTRUE),
175 fHistEvtSelection(0x0),
176 fh2Pt1Pt2C10(0x0),
177 fh2Pt1Pt2C40(0x0), 
178 fh3LocalCoordinates(0x0), 
179 fh2Sum2pt20(0x0),
180 fh2Sum4pt20(0x0),
181 fh2Sum6pt20(0x0),
182 fh2Sum8pt20(0x0),
183 fh2Sum12pt20(0x0),
184 fh2Sum2lpt20(0x0),
185 fh2Sum4lpt20(0x0),
186 fh2Sum6lpt20(0x0),
187 fh2Sum8lpt20(0x0),
188 fh2Sum12lpt20(0x0),
189 fh2Sum2pt40(0x0),
190 fh2Sum4pt40(0x0),
191 fh2Sum6pt40(0x0),
192 fh2Sum8pt40(0x0),
193 fh2Sum12pt40(0x0),
194 fh2Sum2lpt40(0x0),
195 fh2Sum4lpt40(0x0),
196 fh2Sum6lpt40(0x0),
197 fh2Sum8lpt40(0x0),
198 fh2Sum12lpt40(0x0),
199 fhnDeltaR(0x0)
200  {
201    // Constructor
202
203
204    
205
206
207
208
209    fJetBranchName[0] = "";
210    fJetBranchName[1] = "";
211
212    fListJets[0] = new TList;
213    fListJets[1] = new TList;
214
215    DefineOutput(1, TList::Class());
216 }
217
218 AliAnalysisTaskAj::~AliAnalysisTaskAj()
219 {
220    delete fListJets[0];
221    delete fListJets[1];
222 }
223
224 void AliAnalysisTaskAj::SetBranchNames(const TString &branch1, const TString &branch2)
225 {
226    fJetBranchName[0] = branch1;
227    fJetBranchName[1] = branch2;
228 }
229
230 void AliAnalysisTaskAj::Init()
231 {
232
233    // check for jet branches
234    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
235       AliError("Jet branch name not set.");
236    }
237
238 }
239
240 void AliAnalysisTaskAj::UserCreateOutputObjects()
241 {
242    // Create histograms
243    // Called once
244    OpenFile(1);
245    if(!fOutputList) fOutputList = new TList;
246    fOutputList->SetOwner(kTRUE);
247
248    Bool_t oldStatus = TH1::AddDirectoryStatus();
249    TH1::AddDirectory(kFALSE);
250
251
252    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
253    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
254    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
255    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
256    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
257    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
258    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
259
260      UInt_t entries = 0; // bit coded, see GetDimParams() below 
261      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7|1<<8; 
262      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
263
264     
265      //change binning in centrality
266      Double_t *xPt0 = new Double_t[4];
267      xPt0[0] = 0.;
268      for(int i = 1; i<=3;i++){
269       if(xPt0[i-1]<10)xPt0[i] = xPt0[i-1] + 10; // 1 - 5
270       else if(xPt0[i-1]<40)xPt0[i] = xPt0[i-1] + 30; // 5 - 12
271       else xPt0[i] = xPt0[i-1] + 60.; // 18 
272      }
273      fhnDeltaR->SetBinEdges(0,xPt0);
274      delete [] xPt0;
275
276
277
278      //change binning in jTtrack
279      Double_t *xPt3 = new Double_t[3];
280      xPt3[0] = 0.;
281      for(int i = 1; i<=2;i++){
282       if(xPt3[i-1]<1)xPt3[i] = xPt3[i-1] + 1.; // 1 - 5
283       else xPt3[i] = xPt3[i-1] + 149.; // 18 
284      }
285      fhnDeltaR->SetBinEdges(3,xPt3);
286      delete [] xPt3;
287
288
289
290      //change binning in pTtrack
291      Double_t *xPt4 = new Double_t[5];
292      xPt4[0] = 0.;
293      for(int i = 1; i<=4;i++){
294       if(xPt4[i-1]<0.4)xPt4[i] = xPt4[i-1] + 0.4; // 1 - 5
295       else if(xPt4[i-1]<3)xPt4[i] = xPt4[i-1] + 2.6; // 5 - 12
296       else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 7.4; // 5 - 12
297       else xPt4[i] = xPt4[i-1] + 150.; // 18 
298      }
299     fhnDeltaR->SetBinEdges(4,xPt4);
300     delete [] xPt4;
301
302
303
304
305
306
307     //change binning in HTI
308      Double_t *xPt5 = new Double_t[4];
309      xPt5[0] = 0.;
310      for(int i = 1; i<=3;i++){
311       if(xPt5[i-1]<10)xPt5[i] = xPt5[i-1] + 5; // 10 - 12
312       else xPt5[i] = xPt5[i-1] + 40.; // 13
313      }
314     fhnDeltaR->SetBinEdges(8,xPt5);
315     delete [] xPt5;
316
317     
318    
319     fh2Pt1Pt2C10 = new TH2F("Dijet spectra central","",20,0.,200.,20,0.,200.);
320     fh2Pt1Pt2C40 = new TH2F("Dijet spectra peripheral","",20,0.,200.,20,0.,200.);
321     fh3LocalCoordinates = new TH3F("Local coordinates","",10,-2,2,10,-2,2,10,0,100);
322     fh2Sum2pt20 = new TH2F("pL R<0.2 pt20","",10,0.,1.,100,0.,100.);
323     fh2Sum4pt20 = new TH2F("pL R<0.4 pt20","",10,0.,1.,100,0.,100.);
324     fh2Sum6pt20 = new TH2F("pL R<0.6 pt20","",10,0.,1.,100,0.,100.);
325     fh2Sum8pt20 = new TH2F("pL R<0.8 pt20","",10,0.,1.,100,0.,100.);
326     fh2Sum12pt20 = new TH2F("pL R<1.2 pt20","",10,0.,1.,100,0.,100.);
327     fh2Sum2lpt20 = new TH2F("pL R<0.2 low pt pt20","",10,0.,1.,100,0,100);
328     fh2Sum4lpt20 = new TH2F("pL R<0.4 low pt pt20","",10,0.,1.,100,0,100);
329     fh2Sum6lpt20 = new TH2F("pL R<0.6 low pt pt20","",10,0.,1.,100,0,100);
330     fh2Sum8lpt20 = new TH2F("pL R<0.8 low pt pt20","",10,0.,1.,100,0,100);
331     fh2Sum12lpt20 = new TH2F("pL R<1.2 low pt pt20","",10,0.,1.,100,0,100);
332     fh2Sum2pt40 = new TH2F("pL R<0.2 pt40","",10,0.,1.,100,0.,100.);
333     fh2Sum4pt40 = new TH2F("pL R<0.4 pt40","",10,0.,1.,100,0.,100.);
334     fh2Sum6pt40 = new TH2F("pL R<0.6 pt40","",10,0.,1.,100,0.,100.);
335     fh2Sum8pt40 = new TH2F("pL R<0.8 pt40","",10,0.,1.,100,0.,100.);
336     fh2Sum12pt40 = new TH2F("pL R<1.2 pt40","",10,0.,1.,100,0.,100.);
337     fh2Sum2lpt40 = new TH2F("pL R<0.2 low pt pt40","",10,0.,1.,100,0,100);
338     fh2Sum4lpt40 = new TH2F("pL R<0.4 low pt pt40","",10,0.,1.,100,0,100);
339     fh2Sum6lpt40 = new TH2F("pL R<0.6 low pt pt40","",10,0.,1.,100,0,100);
340     fh2Sum8lpt40 = new TH2F("pL R<0.8 low pt pt40","",10,0.,1.,100,0,100);
341     fh2Sum12lpt40 = new TH2F("pL R<1.2 low pt pt40","",10,0.,1.,100,0,100);
342
343
344    fOutputList->Add(fHistEvtSelection);
345    fOutputList->Add(fh2Pt1Pt2C10);
346    fOutputList->Add(fh2Pt1Pt2C40);
347    fOutputList->Add(fh3LocalCoordinates);
348
349    fOutputList->Add(fh2Sum2pt20);
350    fOutputList->Add(fh2Sum4pt20);
351    fOutputList->Add(fh2Sum6pt20);      
352    fOutputList->Add(fh2Sum8pt20);
353    fOutputList->Add(fh2Sum12pt20);  
354    fOutputList->Add(fh2Sum2lpt20);
355    fOutputList->Add(fh2Sum4lpt20);
356    fOutputList->Add(fh2Sum6lpt20);      
357    fOutputList->Add(fh2Sum8lpt20);
358    fOutputList->Add(fh2Sum12lpt20);  
359
360    fOutputList->Add(fh2Sum2pt40);
361    fOutputList->Add(fh2Sum4pt40);
362    fOutputList->Add(fh2Sum6pt40);      
363    fOutputList->Add(fh2Sum8pt40);
364    fOutputList->Add(fh2Sum12pt40);  
365    fOutputList->Add(fh2Sum2lpt40);
366    fOutputList->Add(fh2Sum4lpt40);
367    fOutputList->Add(fh2Sum6lpt40);      
368    fOutputList->Add(fh2Sum8lpt40);
369    fOutputList->Add(fh2Sum12lpt40);  
370
371    fOutputList->Add(fhnDeltaR);
372         
373       
374      
375    // fOutputList->Add(fh3specbiased);
376
377    // =========== Switch on Sumw2 for all histos ===========
378    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
379       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
380       if (h1){
381          h1->Sumw2();
382          continue;
383       }
384     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
385       if (hn){
386          hn->Sumw2();
387       }   
388    }
389    TH1::AddDirectory(oldStatus);
390
391    PostData(1, fOutputList);
392 }
393
394 void AliAnalysisTaskAj::UserExec(Option_t *)
395 {
396    
397
398    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
399       AliError("Jet branch name not set.");
400       return;
401    }
402
403    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
404    if (!fESD) {
405       AliError("ESD not available");
406       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
407       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
408
409       static AliAODEvent* aod = 0;
410        // take all other information from the aod we take the tracks from
411        if(!aod){
412        if(!fESD)aod = fAODIn;
413        else aod = fAODOut;}
414   
415
416  
417     if(fNonStdFile.Length()!=0){
418     // case that we have an AOD extension we need can fetch the jets from the extended output
419     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
420     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
421     if(!fAODExtension){
422       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
423     }}
424     
425
426
427
428
429    // -- event selection --
430    fHistEvtSelection->Fill(1); // number of events before event selection
431
432    // physics selection
433    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
434    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
435    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
436    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
437       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
438       fHistEvtSelection->Fill(2);
439       PostData(1, fOutputList);
440       return;
441    }
442
443    // vertex selection
444    if(!aod){
445      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
446      fHistEvtSelection->Fill(3);
447       PostData(1, fOutputList);
448    }
449    AliAODVertex* primVtx = aod->GetPrimaryVertex();
450
451    if(!primVtx){
452      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
453      fHistEvtSelection->Fill(3);
454      PostData(1, fOutputList);
455      return;
456    }
457
458    Int_t nTracksPrim = primVtx->GetNContributors();
459    if ((nTracksPrim < fMinContribVtx) ||
460          (primVtx->GetZ() < fVtxZMin) ||
461          (primVtx->GetZ() > fVtxZMax) ){
462       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
463       fHistEvtSelection->Fill(3);
464       PostData(1, fOutputList);
465       return;
466    }
467
468    // event class selection (from jet helper task)
469    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
470    if(fDebug) Printf("Event class %d", eventClass);
471    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
472       fHistEvtSelection->Fill(4);
473       PostData(1, fOutputList);
474       return;
475    }
476
477    // centrality selection
478    AliCentrality *cent = 0x0;
479    Double_t centValue = 0.; 
480    if(fESD) {cent = fESD->GetCentrality();
481      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
482    else     centValue=((AliVAODHeader*)aod->GetHeader())->GetCentrality();
483    
484    if(fDebug) printf("centrality: %f\n", centValue);
485    //  if (centValue < fCentMin || centValue > fCentMax){
486    //  fHistEvtSelection->Fill(4);
487    //  PostData(1, fOutputList);
488    //  return;
489    // }
490   
491
492    fHistEvtSelection->Fill(0); 
493    // accepted events  
494    // -- end event selection --
495   
496    // get background
497    AliAODJetEventBackground* externalBackground = 0;
498    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
499       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
500       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
501    }
502    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
503      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
504       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
505    }
506    
507     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
508       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
509       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
510     } 
511   Float_t rho = 0;
512    if(externalBackground)rho = externalBackground->GetBackground(0);
513
514
515    // fetch jets
516    TClonesArray *aodJets[2];
517    aodJets[0]=0;
518    if(fAODOut&&!aodJets[0]){
519    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); 
520    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
521    if(fAODExtension && !aodJets[0]){ 
522    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
523    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
524    if(fAODIn&&!aodJets[0]){
525    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data())); 
526    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data())); }
527
528    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
529    //Int_t inord[aodJets[0]->GetEntriesFast()];
530    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
531    //  ptsub[n]=0;
532    //  inord[n]=0;}   
533
534    TList ParticleList;
535    Int_t nT = GetListOfTracks(&ParticleList);
536      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
537       fListJets[iJetType]->Clear();
538       if (!aodJets[iJetType]) continue;
539
540       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
541       
542    
543       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
544          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
545          if (jet) fListJets[iJetType]->Add(jet);
546          // if(iJetType==0){
547          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
548       }}
549    
550
551    Double_t eta1=0;
552    Int_t selec=-1;
553    Double_t ptmax=-10; 
554    Double_t areaj=0;
555    //   Double_t phij=0;
556    Double_t etaj=0;
557    Double_t ptj=0;
558    Double_t ptcorrj=0;  
559    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
560            AliAODJet* jetj = (AliAODJet*)(fListJets[0]->At(i));
561            etaj  = jetj->Eta();
562            //           phij  = jetj->Phi();
563            ptj = jetj->Pt();
564            if(ptj==0) continue; 
565            areaj = jetj->EffectiveAreaCharged();
566            ptcorrj=ptj-rho*areaj;
567            if(ptcorrj<=0) continue;
568            if((etaj<fJetEtaMin)||(eta1>fJetEtaMax)) continue;
569            if(ptcorrj>ptmax){ptmax=ptcorrj;
570                              selec=i;}}
571    ///hardest jet selected
572     if(selec<0){PostData(1, fOutputList);
573                return;} 
574     AliAODJet* jet1 = (AliAODJet*)(fListJets[0]->At(selec));
575     //What is the hardest constituent track?
576                        AliAODTrack* leadtrack1; 
577                        Int_t ippt=0;
578                        Double_t ppt=-10;   
579                        TRefArray *genTrackList = jet1->GetRefTracks();
580                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
581                        AliAODTrack* genTrack;
582                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
583                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
584                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
585                          ippt=ir;}}
586                         leadtrack1=(AliAODTrack*)(genTrackList->At(ippt));
587                         //If it doesn't exist or if it is greater that 100 GeV, discard.
588                         if(!leadtrack1)  {PostData(1, fOutputList);
589                           return;}
590                         if(leadtrack1->Pt()>=100.){ PostData(1, fOutputList);
591                           return;}             
592
593            //Look to the back-to-back jet 
594            Int_t btb=-1;
595            for(Int_t j=1;j<fListJets[0]->GetEntries();j++){
596            if(j==selec) continue;
597            AliAODJet* jetb = (AliAODJet*)(fListJets[0]->At(j));
598            etaj  = jetb->Eta();
599            //           phij  = jetb->Phi();
600            ptj   = jetb->Pt();
601            if(ptj<=0) continue; 
602            areaj = jetb->EffectiveAreaCharged();
603            ptcorrj=ptj-rho*areaj;
604            if(ptcorrj<=0) continue;
605            if((etaj<fJetEtaMin)||(etaj>fJetEtaMax)) continue;
606            Double_t dphij=RelativePhi(jetb->Phi(),jet1->Phi());  
607            if(TMath::Abs(dphij)>TMath::Pi()-0.2) { btb=j;
608              break;}}
609
610             AliAODJet* jet2 = (AliAODJet*)(fListJets[0]->At(btb));
611            //the back-to-back jet is also identified
612            
613            if(btb<0){PostData(1, fOutputList);
614                           return;} 
615
616             Double_t ptcorr1=jet1->Pt()-rho*jet1->EffectiveAreaCharged();
617             Double_t ptcorr2=jet2->Pt()-rho*jet2->EffectiveAreaCharged();
618           
619             if(centValue<10.) fh2Pt1Pt2C10->Fill(ptcorr1,ptcorr2);
620             if(centValue>40.) fh2Pt1Pt2C40->Fill(ptcorr1,ptcorr2);
621  
622               
623            Double_t px2=jet2->Px();
624            Double_t py2=jet2->Py();
625            Double_t pz2=jet2->Pz();  
626            Double_t phi2=jet2->Phi();
627            Double_t eta2=jet2->Eta();
628            //Once we have have a dijet event,look to the structure of the back-to-back jet:
629
630       TVector3  ppJ1(px2, py2, pz2);
631       TVector3  ppJ3(- px2 * pz2, - py2 * pz2, px2 * px2 + py2 * py2);
632       ppJ3.SetMag(1.);
633       TVector3  ppJ2(-py2, px2, 0);
634       ppJ2.SetMag(1.);
635       Float_t mxx    = 0.;
636       Float_t myy    = 0.;
637       Float_t mxy    = 0.;
638       Int_t   nc     = 0;
639       Float_t sump2  = 0.;
640       Float_t ptMax  = 0.;
641       Float_t etaMax = 0.;
642       Float_t phiMax = 0.;
643       //      Int_t   iMax   = -1;
644       //1st loop over all tracks      
645          for(int it = 0;it<nT;++it){
646          AliVParticle *track = (AliVParticle*)ParticleList.At(it);
647          TVector3 pp(track->Px(), track->Py(), track->Pz());
648           Float_t phi = track->Phi();
649           Float_t eta = track->Eta();
650           Float_t pt  = track->Pt();
651           Float_t jT  = pp.Perp(ppJ1);
652           Float_t deta = eta - eta2;
653           Float_t dphi = phi - phi2;
654           if (dphi >   TMath::Pi()) dphi =   2. * TMath::Pi() - dphi;
655           if (dphi < - TMath::Pi()) dphi = - 2. * TMath::Pi() - dphi;
656           Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
657           /////////To compute the TM axis, we use particles with large jT
658          
659               //cout<<"Jt spectrum............"<<ptbig<<" "<<jT<<endl;
660           /////////We compute the TM with large jT particles only
661              if((jT >1.)&&(r<1.2)){
662             //longitudinal and perpendicular component of the track pT in the
663             //local frame
664               TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
665               TVector3 pPerp = pp - pLong;
666               //projection onto the two perpendicular vectors defined above
667               Float_t ppjX = pPerp.Dot(ppJ2);
668               Float_t ppjY = pPerp.Dot(ppJ3);
669               Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
670               //components of the sphericity matrix
671               mxx += (ppjX * ppjX / ppjT);
672               myy += (ppjY * ppjY / ppjT);
673               mxy += (ppjX * ppjY / ppjT);
674               nc++;
675               sump2 += ppjT;
676               // max pt
677               if (pt > ptMax) {
678                   ptMax  = pt;
679                   //              iMax   = it;
680                   etaMax = deta;
681                   phiMax = dphi;
682               }
683                 } // R < 0.4
684          } // 1st Track Loop
685             
686   // At this point we have mxx, myy, mxy
687       if (nc == 0) return;      
688 // Shericity Matrix     
689       const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};     
690       TMatrixDSym m0(2,ele);
691 // Find eigenvectors
692       TMatrixDSymEigen m(m0);
693       TVectorD eval(2);
694       TMatrixD evecm = m.GetEigenVectors();
695       eval  = m.GetEigenValues();
696 // Largest eigenvector
697       Int_t jev = 0;
698       if (eval[0] < eval[1]) jev = 1;
699       TVectorD evec0(2);
700 // Principle axis
701       evec0 = TMatrixDColumn(evecm, jev);
702       TVector2 evec(evec0[0], evec0[1]); 
703 // Principle axis from leading partice
704       Float_t phiM = TMath::ATan2(phiMax, etaMax);
705       TVector2 evecM(TMath::Cos(phiM), TMath::Sin(phiM)); 
706       Float_t phistM = evecM.DeltaPhi(evec);
707       if (TMath::Abs(phistM) > TMath::Pi()/2.) evec*=(-1.);
708            
709       //////we have now the direction 
710       /////along which the sum of the projections of the particle
711       ///momentum is higher.
712       Double_t sum2lpt20=0;
713       Double_t sum4lpt20=0;
714       Double_t sum6lpt20=0;
715       Double_t sum8lpt20=0;
716       Double_t sum12lpt20=0;
717       Double_t sum2pt20=0;
718       Double_t sum4pt20=0;
719       Double_t sum6pt20=0;
720       Double_t sum8pt20=0;
721       Double_t sum12pt20=0;
722
723        Double_t sum2lpt40=0;
724       Double_t sum4lpt40=0;
725       Double_t sum6lpt40=0;
726       Double_t sum8lpt40=0;
727       Double_t sum12lpt40=0;
728       Double_t sum2pt40=0;
729       Double_t sum4pt40=0;
730       Double_t sum6pt40=0;
731       Double_t sum8pt40=0;
732       Double_t sum12pt40=0;
733
734           for (Int_t ip = 0; ip < nT; ip++) {
735           AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
736           TVector3 pp(track->Px(), track->Py(), track->Pz());
737           Float_t phi = track->Phi();
738           Float_t eta = track->Eta();
739           Float_t pt  = track->Pt();
740           Float_t jT  = pp.Perp(ppJ1);
741           Double_t deta = eta - eta2;
742           Double_t deltaPhi = phi - phi2;
743           Double_t r = TMath::Sqrt(deltaPhi * deltaPhi + deta * deta);
744           
745             if(ptcorr2>20. && ptcorr2<40.){ 
746             if(pt<0.4){
747                  if(r<0.2) sum2lpt20=sum2lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
748                  if(r<0.4) sum4lpt20=sum4lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());   
749                  if(r<0.6) sum6lpt20=sum6lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
750                  if(r<0.8) sum8lpt20=sum8lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());
751                  if(r<1.2) sum12lpt20=sum12lpt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
752             
753                  if(r<0.2) sum2pt20=sum2pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
754                  if(r<0.4) sum4pt20=sum4pt20-1.*pt*TMath::Cos(phi-jet1->Phi());   
755                  if(r<0.6) sum6pt20=sum6pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
756                  if(r<0.8) sum8pt20=sum8pt20-1.*pt*TMath::Cos(phi-jet1->Phi());
757                  if(r<1.2) sum12pt20=sum12pt20-1.*pt*TMath::Cos(phi-jet1->Phi());}
758
759
760                 if(ptcorr2>40. && ptcorr2<60.){ 
761             if(pt<0.4){
762                  if(r<0.2) sum2lpt40=sum2lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
763                  if(r<0.4) sum4lpt40=sum4lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());   
764                  if(r<0.6) sum6lpt40=sum6lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
765                  if(r<0.8) sum8lpt40=sum8lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());
766                  if(r<1.2) sum12lpt40=sum12lpt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
767             
768                  if(r<0.2) sum2pt40=sum2pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
769                  if(r<0.4) sum4pt40=sum4pt40-1.*pt*TMath::Cos(phi-jet1->Phi());   
770                  if(r<0.6) sum6pt40=sum6pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
771                  if(r<0.8) sum8pt40=sum8pt40-1.*pt*TMath::Cos(phi-jet1->Phi());
772                  if(r<1.2) sum12pt40=sum12pt40-1.*pt*TMath::Cos(phi-jet1->Phi());}
773          
774
775
776           
777           TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
778           TVector3 pPerp = pp - pLong;
779           Float_t ppjX = pPerp.Dot(ppJ2);
780           Float_t ppjY = pPerp.Dot(ppJ3);
781           TVector2 vr(ppjX, ppjY) ;
782           //and this is the angle between the particle and the TM axis. 
783           Float_t phistr = evec.DeltaPhi(vr);
784
785           //Double_t phistr=vr.Phi()-evec.Phi();
786
787           if(centValue<10.) fh3LocalCoordinates->Fill(ppjX,ppjY,ptcorr2); 
788           Double_t deltaEta = eta2-track->Eta();
789          
790
791           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
792           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
793           Double_t jetEntries[9] = {centValue,ptcorr1,ptcorr2,jT,pt,deltaEta,deltaPhi,phistr,ptMax}; 
794           fhnDeltaR->Fill(jetEntries);
795           }
796           if(centValue<10.){
797             if(ptcorr2>20.){
798               if(ptcorr1>80.){
799           Double_t aj=(ptcorr1-ptcorr2)/(ptcorr1+ptcorr2); 
800           
801           fh2Sum2pt20->Fill(aj,sum2pt20);
802           fh2Sum4pt20->Fill(aj,sum4pt20);
803           fh2Sum6pt20->Fill(aj,sum6pt20);
804           fh2Sum8pt20->Fill(aj,sum8pt20);
805           fh2Sum12pt20->Fill(aj,sum12pt20);  
806           fh2Sum2lpt20->Fill(aj,sum2lpt20);
807           fh2Sum4lpt20->Fill(aj,sum4lpt20);
808           fh2Sum6lpt20->Fill(aj,sum6lpt20);
809           fh2Sum8lpt20->Fill(aj,sum8lpt20);
810           fh2Sum12lpt20->Fill(aj,sum12lpt20);
811           
812           fh2Sum2pt40->Fill(aj,sum2pt40);
813           fh2Sum4pt40->Fill(aj,sum4pt40);
814           fh2Sum6pt40->Fill(aj,sum6pt40);
815           fh2Sum8pt40->Fill(aj,sum8pt40);
816           fh2Sum12pt40->Fill(aj,sum12pt40);  
817           fh2Sum2lpt40->Fill(aj,sum2lpt40);
818           fh2Sum4lpt40->Fill(aj,sum4lpt40);
819           fh2Sum6lpt40->Fill(aj,sum6lpt40);
820           fh2Sum8lpt40->Fill(aj,sum8lpt40);
821           fh2Sum12lpt40->Fill(aj,sum12lpt40);
822               }}}
823
824
825
826
827
828
829
830    PostData(1, fOutputList);
831 }
832
833 void AliAnalysisTaskAj::Terminate(const Option_t *)
834 {
835    // Draw result to the screen
836    // Called once at the end of the query
837
838    if (!GetOutputData(1))
839    return;
840 }
841
842
843
844
845
846
847
848
849
850
851
852 Int_t  AliAnalysisTaskAj::GetListOfTracks(TList *list){
853
854     Int_t iCount = 0;
855    
856      AliAODEvent *aod = 0;
857      if(!fESD)aod = fAODIn;
858      else aod = fAODOut;   
859      if(!aod) return iCount;
860     
861     for(int it = 0;it < aod->GetNumberOfTracks();++it){
862       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
863       if(!tr) AliFatal("Not a standard AOD");
864       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
865       if(TMath::Abs(tr->Eta())>0.9)continue;
866       if(tr->Pt()<0.15)continue;
867       list->Add(tr);
868       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
869       iCount++;
870     }
871   
872    
873   return iCount;
874  
875 }
876
877    Int_t  AliAnalysisTaskAj::GetHardestTrackBackToJet(AliAODJet *jetbig){
878        
879       Int_t index=-1;
880      AliAODEvent *aod = 0;
881      if(!fESD)aod = fAODIn;
882      else aod = fAODOut;   
883      if(!aod) return index;
884    
885    
886     Double_t ptmax=-10;
887     Double_t dphi=0;
888     //    Double_t dif=0;
889     Int_t iCount=0;
890     for(int it = 0;it < aod->GetNumberOfTracks();++it){
891       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
892       if(!tr) AliFatal("Not a standard AOD");
893       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
894       if(TMath::Abs(tr->Eta())>0.9)continue;
895       if(tr->Pt()<0.15)continue;
896       iCount=iCount+1;
897       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
898       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
899       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
900       index=iCount-1;
901       //  dif=dphi; 
902       }}
903   
904       return index;
905
906    }
907
908
909
910
911
912
913
914
915
916  Int_t  AliAnalysisTaskAj::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
917
918     Int_t iCount = 0;
919      AliAODEvent *aod = 0;
920      if(!fESD)aod = fAODIn;
921      else aod = fAODOut;   
922      if(!aod) return iCount;
923         
924  
925   
926     for(int it = 0;it < aod->GetNumberOfTracks();++it){
927       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
928       if(!tr) AliFatal("Not a standard AOD");
929       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
930       if(TMath::Abs(tr->Eta())>0.9)continue;
931       if(tr->Pt()<0.15)continue;
932       Double_t disR=jetbig->DeltaR(tr);
933       if(disR>0.8)  continue;
934       list->Add(tr);
935       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
936       iCount++;
937     }
938   
939    list->Sort();
940    return iCount;
941
942 }
943
944
945
946
947
948
949
950
951
952
953
954 Int_t AliAnalysisTaskAj::GetNInputTracks()
955 {
956
957    Int_t nInputTracks = 0;
958       AliAODEvent *aod = 0;
959      if(!fESD)aod = fAODIn;
960      else aod = fAODOut;   
961      if(!aod) return nInputTracks;
962    
963
964    TString jbname(fJetBranchName[1]);
965    //needs complete event, use jets without background subtraction
966    for(Int_t i=1; i<=3; ++i){
967       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
968    }
969    // use only HI event
970    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
971    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
972
973    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
974    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
975    if(!tmpAODjets){
976       Printf("Jet branch %s not found", jbname.Data());
977       Printf("AliAnalysisTaskAj::GetNInputTracks FAILED");
978       return -1;
979    }
980    
981    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
982       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
983       if(!jet) continue;
984       TRefArray *trackList = jet->GetRefTracks();
985       Int_t nTracks = trackList->GetEntriesFast();
986       nInputTracks += nTracks;
987       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
988    }
989    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
990
991    return nInputTracks;  
992 }
993
994
995
996 Double_t AliAnalysisTaskAj::RelativePhi(Double_t mphi,Double_t vphi){
997
998   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
999   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1000   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1001   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1002   double dphi = mphi-vphi;
1003   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1004   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1005   return dphi;//dphi in [-Pi, Pi]
1006 }
1007
1008
1009
1010 THnSparse* AliAnalysisTaskAj::NewTHnSparseF(const char* name, UInt_t entries)
1011 {
1012    // generate new THnSparseF, axes are defined in GetDimParams()
1013
1014    Int_t count = 0;
1015    UInt_t tmp = entries;
1016    while(tmp!=0){
1017       count++;
1018       tmp = tmp &~ -tmp;  // clear lowest bit
1019    }
1020
1021    TString hnTitle(name);
1022    const Int_t dim = count;
1023    Int_t nbins[dim];
1024    Double_t xmin[dim];
1025    Double_t xmax[dim];
1026
1027    Int_t i=0;
1028    Int_t c=0;
1029    while(c<dim && i<32){
1030       if(entries&(1<<i)){
1031       
1032          TString label("");
1033          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1034          hnTitle += Form(";%s",label.Data());
1035          c++;
1036       }
1037       
1038       i++;
1039    }
1040    hnTitle += ";";
1041
1042    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1043 }
1044
1045 void AliAnalysisTaskAj::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1046 {
1047    // stores label and binning of axis for THnSparse
1048
1049    const Double_t pi = TMath::Pi();
1050    
1051    switch(iEntry){
1052       
1053    case 0:
1054       label = "V0 centrality (%)";
1055      
1056          nbins = 3;
1057          xmin = 0.;
1058          xmax = 100.;
1059          break;
1060       
1061       
1062    case 1:
1063       label = "corrected jet pt1";
1064          nbins = 10;
1065          xmin = 0.;
1066          xmax = 200.;
1067           break;
1068       
1069       case 2:
1070       label = "corrected jet pt2";
1071          nbins = 10;
1072          xmin = 0.;
1073          xmax = 200.;
1074           break;
1075       
1076
1077       
1078    case 3:
1079       label = "track jT";
1080      
1081          nbins = 2;
1082          xmin = 0.;
1083          xmax = 150.;
1084          break;
1085
1086    case 4:
1087       label = "track pT";
1088      
1089          nbins = 4;
1090          xmin = 0.;
1091          xmax = 150.;
1092          break;
1093       
1094       
1095       case 5:
1096       label = "deltaEta";
1097       nbins = 8;
1098       xmin = -1.6;
1099       xmax = 1.6;
1100       break;
1101
1102
1103       case 6:
1104       label = "deltaPhi";
1105       nbins = 60;
1106       xmin = -0.5*pi;
1107       xmax = 1.5*pi;
1108       break;   
1109    
1110       case 7:
1111       label = "deltaPhiTM";
1112       nbins = 60;
1113       xmin = 0.;
1114       xmax = 1.3*pi;
1115       break;   
1116
1117      
1118         
1119       case 8:
1120       label = "leading track";
1121       nbins = 3;
1122       xmin = 0;
1123       xmax = 50;
1124       break;
1125            
1126
1127
1128
1129    
1130   
1131
1132
1133
1134
1135    }
1136
1137 }
1138