Add sigma2 jet shape and fill constituent info. for subtracted jets
[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=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 = aod->GetTrack(it);
863       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
864       if(TMath::Abs(tr->Eta())>0.9)continue;
865       if(tr->Pt()<0.15)continue;
866       list->Add(tr);
867       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
868       iCount++;
869     }
870   
871    
872   return iCount;
873  
874 }
875
876    Int_t  AliAnalysisTaskAj::GetHardestTrackBackToJet(AliAODJet *jetbig){
877        
878       Int_t index=-1;
879      AliAODEvent *aod = 0;
880      if(!fESD)aod = fAODIn;
881      else aod = fAODOut;   
882      if(!aod) return index;
883    
884    
885     Double_t ptmax=-10;
886     Double_t dphi=0;
887     //    Double_t dif=0;
888     Int_t iCount=0;
889     for(int it = 0;it < aod->GetNumberOfTracks();++it){
890       AliAODTrack *tr = aod->GetTrack(it);
891       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
892       if(TMath::Abs(tr->Eta())>0.9)continue;
893       if(tr->Pt()<0.15)continue;
894       iCount=iCount+1;
895       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
896       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
897       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
898       index=iCount-1;
899       //  dif=dphi; 
900       }}
901   
902       return index;
903
904    }
905
906
907
908
909
910
911
912
913
914  Int_t  AliAnalysisTaskAj::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
915
916     Int_t iCount = 0;
917      AliAODEvent *aod = 0;
918      if(!fESD)aod = fAODIn;
919      else aod = fAODOut;   
920      if(!aod) return iCount;
921         
922  
923   
924     for(int it = 0;it < aod->GetNumberOfTracks();++it){
925       AliAODTrack *tr = aod->GetTrack(it);
926       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
927       if(TMath::Abs(tr->Eta())>0.9)continue;
928       if(tr->Pt()<0.15)continue;
929       Double_t disR=jetbig->DeltaR(tr);
930       if(disR>0.8)  continue;
931       list->Add(tr);
932       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
933       iCount++;
934     }
935   
936    list->Sort();
937    return iCount;
938
939 }
940
941
942
943
944
945
946
947
948
949
950
951 Int_t AliAnalysisTaskAj::GetNInputTracks()
952 {
953
954    Int_t nInputTracks = 0;
955       AliAODEvent *aod = 0;
956      if(!fESD)aod = fAODIn;
957      else aod = fAODOut;   
958      if(!aod) return nInputTracks;
959    
960
961    TString jbname(fJetBranchName[1]);
962    //needs complete event, use jets without background subtraction
963    for(Int_t i=1; i<=3; ++i){
964       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
965    }
966    // use only HI event
967    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
968    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
969
970    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
971    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
972    if(!tmpAODjets){
973       Printf("Jet branch %s not found", jbname.Data());
974       Printf("AliAnalysisTaskAj::GetNInputTracks FAILED");
975       return -1;
976    }
977    
978    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
979       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
980       if(!jet) continue;
981       TRefArray *trackList = jet->GetRefTracks();
982       Int_t nTracks = trackList->GetEntriesFast();
983       nInputTracks += nTracks;
984       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
985    }
986    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
987
988    return nInputTracks;  
989 }
990
991
992
993 Double_t AliAnalysisTaskAj::RelativePhi(Double_t mphi,Double_t vphi){
994
995   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
996   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
997   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
998   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
999   double dphi = mphi-vphi;
1000   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1001   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1002   return dphi;//dphi in [-Pi, Pi]
1003 }
1004
1005
1006
1007 THnSparse* AliAnalysisTaskAj::NewTHnSparseF(const char* name, UInt_t entries)
1008 {
1009    // generate new THnSparseF, axes are defined in GetDimParams()
1010
1011    Int_t count = 0;
1012    UInt_t tmp = entries;
1013    while(tmp!=0){
1014       count++;
1015       tmp = tmp &~ -tmp;  // clear lowest bit
1016    }
1017
1018    TString hnTitle(name);
1019    const Int_t dim = count;
1020    Int_t nbins[dim];
1021    Double_t xmin[dim];
1022    Double_t xmax[dim];
1023
1024    Int_t i=0;
1025    Int_t c=0;
1026    while(c<dim && i<32){
1027       if(entries&(1<<i)){
1028       
1029          TString label("");
1030          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1031          hnTitle += Form(";%s",label.Data());
1032          c++;
1033       }
1034       
1035       i++;
1036    }
1037    hnTitle += ";";
1038
1039    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1040 }
1041
1042 void AliAnalysisTaskAj::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1043 {
1044    // stores label and binning of axis for THnSparse
1045
1046    const Double_t pi = TMath::Pi();
1047    
1048    switch(iEntry){
1049       
1050    case 0:
1051       label = "V0 centrality (%)";
1052      
1053          nbins = 3;
1054          xmin = 0.;
1055          xmax = 100.;
1056          break;
1057       
1058       
1059    case 1:
1060       label = "corrected jet pt1";
1061          nbins = 10;
1062          xmin = 0.;
1063          xmax = 200.;
1064           break;
1065       
1066       case 2:
1067       label = "corrected jet pt2";
1068          nbins = 10;
1069          xmin = 0.;
1070          xmax = 200.;
1071           break;
1072       
1073
1074       
1075    case 3:
1076       label = "track jT";
1077      
1078          nbins = 2;
1079          xmin = 0.;
1080          xmax = 150.;
1081          break;
1082
1083    case 4:
1084       label = "track pT";
1085      
1086          nbins = 4;
1087          xmin = 0.;
1088          xmax = 150.;
1089          break;
1090       
1091       
1092       case 5:
1093       label = "deltaEta";
1094       nbins = 8;
1095       xmin = -1.6;
1096       xmax = 1.6;
1097       break;
1098
1099
1100       case 6:
1101       label = "deltaPhi";
1102       nbins = 60;
1103       xmin = -0.5*pi;
1104       xmax = 1.5*pi;
1105       break;   
1106    
1107       case 7:
1108       label = "deltaPhiTM";
1109       nbins = 60;
1110       xmin = 0.;
1111       xmax = 1.3*pi;
1112       break;   
1113
1114      
1115         
1116       case 8:
1117       label = "leading track";
1118       nbins = 3;
1119       xmin = 0;
1120       xmax = 50;
1121       break;
1122            
1123
1124
1125
1126    
1127   
1128
1129
1130
1131
1132    }
1133
1134 }
1135