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