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