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 // *******************************************
8 /**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
30 #include "THnSparse.h"
32 #include "TMatrixDSym.h"
33 #include "TMatrixDSymEigen.h"
39 #include "AliAnalysisTask.h"
40 #include "AliAnalysisManager.h"
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"
54 #include "AliAnalysisTaskAj.h"
59 ClassImp(AliAnalysisTaskAj)
61 AliAnalysisTaskAj::AliAnalysisTaskAj() :
67 fBackgroundBranch(""),
70 fOfflineTrgMask(AliVEvent::kAny),
83 fAngStructCloseTracks(0),
91 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
92 fJetPtFractionMin(0.5),
100 fHistEvtSelection(0x0),
103 fh3LocalCoordinates(0x0),
126 // default Constructor
129 fJetBranchName[0] = "";
130 fJetBranchName[1] = "";
132 fListJets[0] = new TList;
133 fListJets[1] = new TList;
136 AliAnalysisTaskAj::AliAnalysisTaskAj(const char *name) :
137 AliAnalysisTaskSE(name),
142 fBackgroundBranch(""),
145 fOfflineTrgMask(AliVEvent::kAny),
157 fNInputTracksMax(-1),
158 fAngStructCloseTracks(0),
166 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
167 fJetPtFractionMin(0.5),
175 fHistEvtSelection(0x0),
178 fh3LocalCoordinates(0x0),
209 fJetBranchName[0] = "";
210 fJetBranchName[1] = "";
212 fListJets[0] = new TList;
213 fListJets[1] = new TList;
215 DefineOutput(1, TList::Class());
218 AliAnalysisTaskAj::~AliAnalysisTaskAj()
224 void AliAnalysisTaskAj::SetBranchNames(const TString &branch1, const TString &branch2)
226 fJetBranchName[0] = branch1;
227 fJetBranchName[1] = branch2;
230 void AliAnalysisTaskAj::Init()
233 // check for jet branches
234 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
235 AliError("Jet branch name not set.");
240 void AliAnalysisTaskAj::UserCreateOutputObjects()
245 if(!fOutputList) fOutputList = new TList;
246 fOutputList->SetOwner(kTRUE);
248 Bool_t oldStatus = TH1::AddDirectoryStatus();
249 TH1::AddDirectory(kFALSE);
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)");
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);
265 //change binning in centrality
266 Double_t *xPt0 = new Double_t[4];
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
273 fhnDeltaR->SetBinEdges(0,xPt0);
278 //change binning in jTtrack
279 Double_t *xPt3 = new Double_t[3];
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
285 fhnDeltaR->SetBinEdges(3,xPt3);
290 //change binning in pTtrack
291 Double_t *xPt4 = new Double_t[5];
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
299 fhnDeltaR->SetBinEdges(4,xPt4);
307 //change binning in HTI
308 Double_t *xPt5 = new Double_t[4];
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
314 fhnDeltaR->SetBinEdges(8,xPt5);
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);
344 fOutputList->Add(fHistEvtSelection);
345 fOutputList->Add(fh2Pt1Pt2C10);
346 fOutputList->Add(fh2Pt1Pt2C40);
347 fOutputList->Add(fh3LocalCoordinates);
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);
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);
371 fOutputList->Add(fhnDeltaR);
375 // fOutputList->Add(fh3specbiased);
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));
384 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
389 TH1::AddDirectory(oldStatus);
391 PostData(1, fOutputList);
394 void AliAnalysisTaskAj::UserExec(Option_t *)
398 if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
399 AliError("Jet branch name not set.");
403 fESD=dynamic_cast<AliESDEvent*>(InputEvent());
405 AliError("ESD not available");
406 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());}
407 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
409 static AliAODEvent* aod = 0;
410 // take all other information from the aod we take the tracks from
412 if(!fESD)aod = fAODIn;
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);
422 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
429 // -- event selection --
430 fHistEvtSelection->Fill(1); // number of events before event 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);
445 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
446 fHistEvtSelection->Fill(3);
447 PostData(1, fOutputList);
449 AliAODVertex* primVtx = aod->GetPrimaryVertex();
452 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
453 fHistEvtSelection->Fill(3);
454 PostData(1, fOutputList);
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);
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);
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();
484 if(fDebug) printf("centrality: %f\n", centValue);
485 // if (centValue < fCentMin || centValue > fCentMax){
486 // fHistEvtSelection->Fill(4);
487 // PostData(1, fOutputList);
492 fHistEvtSelection->Fill(0);
494 // -- end event selection --
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());;
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());;
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());;
512 if(externalBackground)rho = externalBackground->GetBackground(0);
516 TClonesArray *aodJets[2];
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())); }
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++){
535 Int_t nT = GetListOfTracks(&ParticleList);
536 for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
537 fListJets[iJetType]->Clear();
538 if (!aodJets[iJetType]) continue;
540 if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
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);
547 // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
559 for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
560 AliAODJet* jetj = (AliAODJet*)(fListJets[0]->At(i));
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;
571 ///hardest jet selected
572 if(selec<0){PostData(1, fOutputList);
574 AliAODJet* jet1 = (AliAODJet*)(fListJets[0]->At(selec));
575 //What is the hardest constituent track?
576 AliAODTrack* leadtrack1;
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();
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);
590 if(leadtrack1->Pt()>=100.){ PostData(1, fOutputList);
593 //Look to the back-to-back jet
595 for(Int_t j=1;j<fListJets[0]->GetEntries();j++){
596 if(j==selec) continue;
597 AliAODJet* jetb = (AliAODJet*)(fListJets[0]->At(j));
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;
610 AliAODJet* jet2 = (AliAODJet*)(fListJets[0]->At(btb));
611 //the back-to-back jet is also identified
613 if(btb<0){PostData(1, fOutputList);
616 Double_t ptcorr1=jet1->Pt()-rho*jet1->EffectiveAreaCharged();
617 Double_t ptcorr2=jet2->Pt()-rho*jet2->EffectiveAreaCharged();
619 if(centValue<10.) fh2Pt1Pt2C10->Fill(ptcorr1,ptcorr2);
620 if(centValue>40.) fh2Pt1Pt2C40->Fill(ptcorr1,ptcorr2);
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:
630 TVector3 ppJ1(px2, py2, pz2);
631 TVector3 ppJ3(- px2 * pz2, - py2 * pz2, px2 * px2 + py2 * py2);
633 TVector3 ppJ2(-py2, px2, 0);
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
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
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);
686 // At this point we have mxx, myy, mxy
689 const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
690 TMatrixDSym m0(2,ele);
692 TMatrixDSymEigen m(m0);
694 TMatrixD evecm = m.GetEigenVectors();
695 eval = m.GetEigenValues();
696 // Largest eigenvector
698 if (eval[0] < eval[1]) jev = 1;
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.);
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;
721 Double_t sum12pt20=0;
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;
732 Double_t sum12pt40=0;
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);
745 if(ptcorr2>20. && ptcorr2<40.){
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());}
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());}
760 if(ptcorr2>40. && ptcorr2<60.){
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());}
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());}
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);
785 //Double_t phistr=vr.Phi()-evec.Phi();
787 if(centValue<10.) fh3LocalCoordinates->Fill(ppjX,ppjY,ptcorr2);
788 Double_t deltaEta = eta2-track->Eta();
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);
799 Double_t aj=(ptcorr1-ptcorr2)/(ptcorr1+ptcorr2);
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);
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);
830 PostData(1, fOutputList);
833 void AliAnalysisTaskAj::Terminate(const Option_t *)
835 // Draw result to the screen
836 // Called once at the end of the query
838 if (!GetOutputData(1))
852 Int_t AliAnalysisTaskAj::GetListOfTracks(TList *list){
856 AliAODEvent *aod = 0;
857 if(!fESD)aod = fAODIn;
859 if(!aod) return iCount;
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;
867 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
876 Int_t AliAnalysisTaskAj::GetHardestTrackBackToJet(AliAODJet *jetbig){
879 AliAODEvent *aod = 0;
880 if(!fESD)aod = fAODIn;
882 if(!aod) return index;
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;
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();
913 Int_t AliAnalysisTaskAj::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
916 AliAODEvent *aod = 0;
917 if(!fESD)aod = fAODIn;
919 if(!aod) return iCount;
923 for(int it = 0;it < aod->GetNumberOfTracks();++it){
924 AliAODTrack *tr = aod->GetTrack(it);
925 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
926 if(TMath::Abs(tr->Eta())>0.9)continue;
927 if(tr->Pt()<0.15)continue;
928 Double_t disR=jetbig->DeltaR(tr);
929 if(disR>0.8) continue;
931 //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
950 Int_t AliAnalysisTaskAj::GetNInputTracks()
953 Int_t nInputTracks = 0;
954 AliAODEvent *aod = 0;
955 if(!fESD)aod = fAODIn;
957 if(!aod) return nInputTracks;
960 TString jbname(fJetBranchName[1]);
961 //needs complete event, use jets without background subtraction
962 for(Int_t i=1; i<=3; ++i){
963 if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
966 if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
967 if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
969 if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
970 TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
972 Printf("Jet branch %s not found", jbname.Data());
973 Printf("AliAnalysisTaskAj::GetNInputTracks FAILED");
977 for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
978 AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
980 TRefArray *trackList = jet->GetRefTracks();
981 Int_t nTracks = trackList->GetEntriesFast();
982 nInputTracks += nTracks;
983 if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
985 if(fDebug) Printf("---> input tracks: %d", nInputTracks);
992 Double_t AliAnalysisTaskAj::RelativePhi(Double_t mphi,Double_t vphi){
994 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
995 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
996 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
997 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
998 double dphi = mphi-vphi;
999 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1000 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1001 return dphi;//dphi in [-Pi, Pi]
1006 THnSparse* AliAnalysisTaskAj::NewTHnSparseF(const char* name, UInt_t entries)
1008 // generate new THnSparseF, axes are defined in GetDimParams()
1011 UInt_t tmp = entries;
1014 tmp = tmp &~ -tmp; // clear lowest bit
1017 TString hnTitle(name);
1018 const Int_t dim = count;
1025 while(c<dim && i<32){
1029 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1030 hnTitle += Form(";%s",label.Data());
1038 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1041 void AliAnalysisTaskAj::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1043 // stores label and binning of axis for THnSparse
1045 const Double_t pi = TMath::Pi();
1050 label = "V0 centrality (%)";
1059 label = "corrected jet pt1";
1066 label = "corrected jet pt2";
1107 label = "deltaPhiTM";
1116 label = "leading track";