PWGJE/UserTasks/AliAnalysisTaskJetAntenna.cxx

   1 // ******************************************
   2 // This task computes several jet observables like
   3 // the fraction of energy in inner and outer coronnas,
   4 // jet-track correlations,triggered jet shapes and
   5 // correlation strength distribution of particles inside jets.
   6 // Author: lcunquei@cern.ch
   7 // *******************************************
   8
   9
  10 /**************************************************************************
  11  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  12  *                                                                        *
  13  * Author: The ALICE Off-line Project.                                    *
  14  * Contributors are mentioned in the code where appropriate.              *
  15  *                                                                        *
  16  * Permission to use, copy, modify and distribute this software and its   *
  17  * documentation strictly for non-commercial purposes is hereby granted   *
  18  * without fee, provided that the above copyright notice appears in all   *
  19  * copies and that both the copyright notice and this permission notice   *
  20  * appear in the supporting documentation. The authors make no claims     *
  21  * about the suitability of this software for any purpose. It is          *
  22  * provided "as is" without express or implied warranty.                  *
  23  **************************************************************************/
  24
  25
  26 #include "TChain.h"
  27 #include "TTree.h"
  28 #include "TMath.h"
  29 #include "TH1F.h"
  30 #include "TH2F.h"
  31 #include "TH3F.h"
  32 #include "THnSparse.h"
  33 #include "TMatrixD.h"
  34 #include "TMatrixDSym.h"
  35 #include "TMatrixDSymEigen.h"
  36 #include "TCanvas.h"
  37 #include "TRandom3.h"
  38 #include "AliLog.h"
  39
  40 #include "AliAnalysisTask.h"
  41 #include "AliAnalysisManager.h"
  42
  43 #include "AliVEvent.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliCentrality.h"
  47 #include "AliAnalysisHelperJetTasks.h"
  48 #include "AliInputEventHandler.h"
  49 #include "AliAODJetEventBackground.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliAnalysisTaskFastEmbedding.h"
  52 #include "AliAODEvent.h"
  53 #include "AliAODHandler.h"
  54 #include "AliAODJet.h"
  55
  56 #include "AliAnalysisTaskJetAntenna.h"
  57
  58 using std::cout;
  59 using std::endl;
  60
  61 ClassImp(AliAnalysisTaskJetAntenna)
  62
  63 AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna() :
  64 AliAnalysisTaskSE(),
  65 fESD(0x0),
  66 fAODIn(0x0),
  67 fAODOut(0x0),
  68 fAODExtension(0x0),
  69 fBackgroundBranch(""),
  70 fNonStdFile(""),
  71 fIsPbPb(kTRUE),
  72 fOfflineTrgMask(AliVEvent::kAny),
  73 fMinContribVtx(1),
  74 fVtxZMin(-10.),
  75 fVtxZMax(10.),
  76 fEvtClassMin(0),
  77 fEvtClassMax(4),
  78 fFilterMask(0),
  79 fFilterMaskBestPt(0),
  80 fFilterType(0),
  81 fCentMin(0.),
  82 fCentMax(100.),
  83 fNInputTracksMin(0),
  84 fNInputTracksMax(-1),
  85 fRequireITSRefit(0),
  86 fApplySharedClusterCut(0),
  87 fTrackTypeRec(kTrackUndef),
  88 fRPAngle(0),
  89 fNRPBins(50),
  90 fSemigoodCorrect(0),
  91 fDoMatching(kFALSE),
  92 fHolePos(4.71),
  93 fHoleWidth(0.2),
  94 fCutTM(0.15),
  95 fJetEtaMin(-.5),
  96 fJetEtaMax(.5),
  97 fNevents(0),
  98 fTindex(0),
  99 fJetPtMin(20.),
 100 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 101 fJetPtFractionMin(0.5),
 102 fNMatchJets(4),
 103 fMatchMaxDist(0.8),
 104 fKeepJets(kFALSE),
 105 fkNbranches(2),
 106 fkEvtClasses(12),
 107 fOutputList(0x0),
 108 fHistEvtSelection(0x0),
 109 fh2JetEntries(0x0),
 110 fh2Circularity(0x0),
 111 fh2JetAxisPhi(0x0),
 112 fhnJetTM(0x0)
 113 {
 114    // default Constructor
 115
 116    fJetBranchName[0] = "";
 117    fJetBranchName[1] = "";
 118
 119    fListJets[0] = new TList;
 120    fListJets[1] = new TList;
 121 }
 122
 123 AliAnalysisTaskJetAntenna::AliAnalysisTaskJetAntenna(const char *name) :
 124 AliAnalysisTaskSE(name),
 125 fESD(0x0),
 126 fAODIn(0x0),
 127 fAODOut(0x0),
 128 fAODExtension(0x0),
 129 fBackgroundBranch(""),
 130 fNonStdFile(""),
 131 fIsPbPb(kTRUE),
 132 fOfflineTrgMask(AliVEvent::kAny),
 133 fMinContribVtx(1),
 134 fVtxZMin(-10.),
 135 fVtxZMax(10.),
 136 fEvtClassMin(0),
 137 fEvtClassMax(4),
 138 fFilterMask(0),
 139 fFilterMaskBestPt(0),
 140 fFilterType(0),
 141 fCentMin(0.),
 142 fCentMax(100.),
 143 fNInputTracksMin(0),
 144 fNInputTracksMax(-1),
 145 fRequireITSRefit(0),
 146 fApplySharedClusterCut(0),
 147 fTrackTypeRec(kTrackUndef),
 148 fRPAngle(0),
 149 fNRPBins(50),
 150 fSemigoodCorrect(0),
 151 fDoMatching(kFALSE),
 152 fHolePos(4.71),
 153 fHoleWidth(0.2),
 154 fCutTM(0.15),
 155 fJetEtaMin(-.5),
 156 fJetEtaMax(.5),
 157 fNevents(0),
 158 fTindex(0),
 159 fJetPtMin(20.),
 160 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 161 fJetPtFractionMin(0.5),
 162 fNMatchJets(4),
 163 fMatchMaxDist(0.8),
 164 fKeepJets(kFALSE),
 165 fkNbranches(2),
 166 fkEvtClasses(12),
 167 fOutputList(0x0),
 168 fHistEvtSelection(0x0),
 169 fh2JetEntries(0x0),
 170 fh2Circularity(0x0),
 171 fh2JetAxisPhi(0x0),
 172 fhnJetTM(0x0)
 173  {
 174    // Constructor
 175
 176
 177    fJetBranchName[0] = "";
 178    fJetBranchName[1] = "";
 179
 180    fListJets[0] = new TList;
 181    fListJets[1] = new TList;
 182
 183    DefineOutput(1, TList::Class());
 184 }
 185
 186 AliAnalysisTaskJetAntenna::~AliAnalysisTaskJetAntenna()
 187 {
 188    delete fListJets[0];
 189    delete fListJets[1];
 190 }
 191
 192 void AliAnalysisTaskJetAntenna::SetBranchNames(const TString &branch1, const TString &branch2)
 193 {
 194    fJetBranchName[0] = branch1;
 195    fJetBranchName[1] = branch2;
 196 }
 197
 198 void AliAnalysisTaskJetAntenna::Init()
 199 {
 200
 201    // check for jet branches
 202    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 203       AliError("Jet branch name not set.");
 204    }
 205
 206 }
 207
 208 void AliAnalysisTaskJetAntenna::UserCreateOutputObjects()
 209 {
 210   // Create histograms
 211   // Called once
 212   OpenFile(1);
 213   if(!fOutputList) fOutputList = new TList;
 214   fOutputList->SetOwner(kTRUE);
 215
 216   Bool_t oldStatus = TH1::AddDirectoryStatus();
 217   TH1::AddDirectory(kFALSE);
 218
 219   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 220   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 221   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 222   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 223   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 224   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 225   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 226   fOutputList->Add(fHistEvtSelection);
 227   fh2JetEntries=new TH2F("JetEntries","",150,0,150,10,-0.5,9.5);
 228   fOutputList->Add(fh2JetEntries);
 229   fh2Circularity=new TH2F("Circcularity","",10,0,1,150,0,150);
 230   fOutputList->Add(fh2Circularity);
 231   fh2JetAxisPhi=new TH2F("JetAxisSmearPhi","",9,0,TMath::Pi(),10,-0.5,9.5);
 232   fOutputList->Add(fh2JetAxisPhi);
 233
 234
 235   Int_t nbinsJet[10]={3,9,75,9,36,5,7,10,50,2};
 236   Double_t binlowJet[10]= {0,0, 0, 0,-0.5*TMath::Pi(),0,0,-0.5,0,0};
 237   Double_t binupJet[10]= {100,0.9, 150,150,1.5*TMath::Pi(),1,150,9.5,200,2};
 238   fhnJetTM = new THnSparseF("fhnJetTM", "fhnJetTM; cent;dr;pt_jet;pt_track;phi;circ;ptrue;pthard;ptlead;isemebed",10,nbinsJet,binlowJet,binupJet);
 239
 240  Double_t *xPt3=new Double_t[10];
 241   xPt3[0] = 0.;
 242   for(Int_t i = 1;i<=9;i++){
 243     if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
 244     else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
 245     else xPt3[i] = xPt3[i-1] + 150.; // 18
 246   }
 247   fhnJetTM->SetBinEdges(3,xPt3);
 248
 249   // Double_t *xPt2=new Double_t[10];
 250   // xPt2[0] = 0.;
 251   // xPt2[1]=20;
 252   // xPt2[2]=40;
 253   // xPt2[3]=60;
 254   // xPt2[4]=80;
 255   // xPt2[5]=100;
 256   // xPt2[6]=120;
 257   // xPt2[7]=150;
 258
 259   //fhnJetTM->SetBinEdges(2,xPt2);
 260
 261
 262   Double_t *xPt4=new Double_t[4];
 263   xPt4[0] = 0.;
 264   xPt4[1]=10;
 265   xPt4[2]=30;
 266   xPt4[3]=50;
 267
 268
 269   fhnJetTM->SetBinEdges(0,xPt4);
 270
 271   Double_t *xPt5=new Double_t[10];
 272   xPt5[0] = 0.;
 273   xPt5[1]=20;
 274   xPt5[2]=40;
 275   xPt5[3]=60;
 276   xPt5[4]=80;
 277   xPt5[5]=100;
 278   xPt5[6]=120;
 279   xPt5[7]=150;
 280
 281   fhnJetTM->SetBinEdges(6,xPt5);
 282
 283
 284    fOutputList->Add(fhnJetTM);
 285    delete [] xPt3;
 286    // delete [] xPt2;
 287    delete [] xPt4;
 288    delete [] xPt5;
 289   // =========== Switch on Sumw2 for all histos ===========
 290   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 291     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 292     if (h1){
 293       h1->Sumw2();
 294       continue;
 295     }
 296     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 297     if (hn){
 298       hn->Sumw2();
 299     }
 300   }
 301   TH1::AddDirectory(oldStatus);
 302
 303   PostData(1, fOutputList);
 304 }
 305
 306 void AliAnalysisTaskJetAntenna::UserExec(Option_t *)
 307 {
 308
 309
 310   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 311     AliError("Jet branch name not set.");
 312     return;
 313   }
 314
 315   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 316   if (!fESD) {
 317     AliError("ESD not available");
 318     fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 319   }
 320   fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 321
 322   static AliAODEvent* aod = 0;
 323   // take all other information from the aod we take the tracks from
 324   if(!aod){
 325     if(!fESD)aod = fAODIn;
 326     else aod = fAODOut;}
 327
 328
 329
 330   if(fNonStdFile.Length()!=0){
 331     // case that we have an AOD extension we need can fetch the jets from the extended output
 332     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 333     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 334     if(!fAODExtension){
 335       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 336     }
 337   }
 338
 339
 340   // -- event selection --
 341   fHistEvtSelection->Fill(1); // number of events before event selection
 342
 343
 344   Bool_t selected=kTRUE;
 345   selected = AliAnalysisHelperJetTasks::Selected();
 346   if(!selected){
 347     // no selection by the service task, we continue
 348     PostData(1,fOutputList);
 349     return;}
 350
 351
 352
 353   // physics selection: this is now redundant, all should appear as accepted after service task selection
 354   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 355     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 356   std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
 357   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 358     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 359     fHistEvtSelection->Fill(2);
 360     PostData(1, fOutputList);
 361     return;
 362   }
 363
 364
 365
 366   // vertex selection
 367   if(!aod){
 368     if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 369     fHistEvtSelection->Fill(3);
 370     PostData(1, fOutputList);
 371   }
 372   AliAODVertex* primVtx = aod->GetPrimaryVertex();
 373
 374   if(!primVtx){
 375     if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 376     fHistEvtSelection->Fill(3);
 377     PostData(1, fOutputList);
 378     return;
 379   }
 380
 381   Int_t nTracksPrim = primVtx->GetNContributors();
 382   if ((nTracksPrim < fMinContribVtx) ||
 383       (primVtx->GetZ() < fVtxZMin) ||
 384       (primVtx->GetZ() > fVtxZMax) ){
 385     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 386     fHistEvtSelection->Fill(3);
 387     PostData(1, fOutputList);
 388     return;
 389   }
 390
 391
 392
 393   // centrality selection
 394   AliCentrality *cent = 0x0;
 395   Double_t centValue = 0.;
 396   if(fIsPbPb){
 397     if(fESD) {cent = fESD->GetCentrality();
 398       if(cent) centValue = cent->GetCentralityPercentile("V0M");}
 399     else     centValue=((AliVAODHeader*)aod->GetHeader())->GetCentrality();
 400
 401     if(fDebug) printf("centrality: %f\n", centValue);
 402     if (centValue < fCentMin || centValue > fCentMax){
 403       fHistEvtSelection->Fill(4);
 404       PostData(1, fOutputList);
 405       return;
 406     }}
 407
 408
 409   fHistEvtSelection->Fill(0);
 410   // accepted events
 411   // -- end event selection --
 412   // pt hard
 413   Double_t pthard=0;
 414   Double_t pthardbin=0;
 415   if(fDoMatching){
 416   pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
 417   pthardbin = GetPtHardBin(pthard);}
 418   // get background
 419   AliAODJetEventBackground* externalBackground = 0;
 420   if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
 421     externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 422     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 423   }
 424   if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 425     externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 426     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 427   }
 428
 429   if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
 430     externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 431     if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 432   }
 433
 434   Float_t rho = 0;
 435
 436
 437   if(externalBackground)rho = externalBackground->GetBackground(0);
 438
 439   // fetch jets
 440   TClonesArray *aodJets[2];
 441   aodJets[0]=0;
 442   if(fAODOut&&!aodJets[0]){
 443     aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
 444     aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));
 445   }
 446   if(fAODExtension && !aodJets[0]){
 447     aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 448     aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
 449   }
 450   if(fAODIn&&!aodJets[0]){
 451     aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
 452     aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));
 453   }
 454
 455
 456
 457   Int_t nT=0;
 458   TList ParticleList;
 459   if(!fDoMatching) nT = GetListOfTracks(&ParticleList);
 460   if(fDoMatching)nT=GetListOfTracksExtra(&ParticleList);
 461   if(nT<0){
 462     PostData(1, fOutputList);
 463     return;
 464   }
 465
 466   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 467     fListJets[iJetType]->Clear();
 468     if (!aodJets[iJetType]) continue;
 469     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 470     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 471       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 472       if (jet) fListJets[iJetType]->Add(jet);
 473     }
 474   }
 475
 476 // jet matching
 477   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 478   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 479   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 480   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 481
 482   if(fDoMatching){
 483
 484
 485   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 486   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 487                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 488                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);}
 489
 490
 491
 492
 493
 494
 495  // loop over matched jets
 496   Int_t ir = -1;  // index of matched reconstruced jet
 497   Float_t fraction = -1.;
 498    AliAODJet *jetmatched =0x0;
 499
 500
 501
 502
 503
 504     for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 505
 506     Double_t etabig=0;
 507     Double_t ptbig=0;
 508     Double_t areabig=0;
 509     Double_t phibig=0.;
 510     Double_t pxbig,pybig,pzbig;
 511     Double_t phitrue=0;
 512     Double_t etatrue=0;
 513     Double_t smearphi=0;
 514     Double_t ptrue=0;
 515     AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 516     etabig  = jetbig->Eta();
 517     phibig  = jetbig->Phi();
 518     phitrue=phibig;
 519     etatrue=etabig;
 520     ptbig   = jetbig->Pt();
 521     if(ptbig==0) continue;
 522     areabig = jetbig->EffectiveAreaCharged();
 523     if(!fDoMatching) ptbig=ptbig-rho*areabig;
 524     if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 525
 526     if(fSemigoodCorrect){
 527       if((phibig>fHolePos-fHoleWidth) && (phibig<fHolePos+fHoleWidth)) continue;
 528     }
 529     pxbig=jetbig->Px();
 530     pybig=jetbig->Py();
 531     pzbig=jetbig->Pz();
 532
 533    if(fDoMatching){
 534      Bool_t jetAccepted=kTRUE;
 535     ir = aMatchIndex[i];
 536     if(ir>=0) jetmatched = (AliAODJet*)(fListJets[1]->At(ir));
 537     else    continue;
 538     fraction  = aPtFraction[i];
 539     // minimum fraction required
 540
 541     if(fraction<fJetPtFractionMin) jetAccepted = kFALSE;
 542     if(jetAccepted){
 543     // jet acceptance + minimum pT check
 544     if(jetmatched->Eta()>fJetEtaMax || jetmatched->Eta()<fJetEtaMin)  jetAccepted = kFALSE;
 545
 546     }
 547     ptrue=ptbig;
 548     if(!jetAccepted) continue;
 549
 550     etabig=jetmatched->Eta();
 551     phibig=jetmatched->Phi();
 552     pxbig=jetmatched->Px();
 553     pybig=jetmatched->Py();
 554     pzbig=jetmatched->Pz();
 555     ptbig=jetmatched->Pt()-rho*jetmatched->EffectiveAreaCharged();
 556
 557
 558
 559
 560
 561     smearphi=RelativePhi(phitrue,jetmatched->Phi());
 562     smearphi=TMath::Abs(smearphi);
 563
 564
 565 }
 566
 567
 568
 569
 570     TVector3  ppJ1(pxbig, pybig, pzbig);
 571     TVector3  ppJ3(- pxbig * pzbig, - pybig * pzbig, pxbig * pxbig + pybig * pybig);
 572     ppJ3.SetMag(1.);
 573     TVector3  ppJ2(-pybig, pxbig, 0);
 574     ppJ2.SetMag(1.);
 575
 576
 577
 578
 579     Float_t mxx    = 0.;
 580     Float_t myy    = 0.;
 581     Float_t mxy    = 0.;
 582     Int_t   nc     = 0;
 583     Float_t sump2  = 0.;
 584     Float_t ptmax=-10;
 585     for(int it = 0;it<ParticleList.GetEntries();++it){
 586       AliVParticle *track = (AliVParticle*)ParticleList.At(it);
 587       TVector3 pp(track->Px(), track->Py(), track->Pz());
 588       Float_t phi = track->Phi();
 589       Float_t eta = track->Eta();
 590       Float_t pt  = track->Pt();
 591       Float_t deta = eta - etabig;
 592       Float_t dphi = RelativePhi(phi,phibig);
 593       if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
 594       Float_t r = TMath::Sqrt(dphi * dphi + deta * deta);
 595
 596        if (r < 0.4 && pt>fCutTM) {
 597         //longitudinal and perpendicular component of the track pT in the
 598         //local frame
 599         if(pt>ptmax) ptmax=pt;
 600         TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 601         TVector3 pPerp = pp - pLong;
 602         //projection onto the two perpendicular vectors defined above
 603         Float_t ppjX = pPerp.Dot(ppJ2);
 604         Float_t ppjY = pPerp.Dot(ppJ3);
 605         Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
 606         //components of the 2D symmetrical sphericity matrix
 607         mxx += (ppjX * ppjX / ppjT);
 608         myy += (ppjY * ppjY / ppjT);
 609         mxy += (ppjX * ppjY / ppjT);
 610         nc++;
 611         sump2 += ppjT;}
 612
 613       if(nc<2) continue;
 614
 615     } // 1st Track Loop
 616
 617
 618     // Sphericity Matrix
 619     const Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
 620     TMatrixDSym m0(2,ele);
 621     // Find eigenvectors
 622     TMatrixDSymEigen m(m0);
 623     TVectorD eval(2);
 624     TMatrixD evecm = m.GetEigenVectors();
 625     eval  = m.GetEigenValues();
 626     // Largest eigenvector
 627     Int_t jev = 0;
 628     if (eval[0] < eval[1]) jev = 1;
 629     TVectorD evec0(2);
 630     // Principle axis
 631     evec0 = TMatrixDColumn(evecm, jev);
 632     TVector2 evec(evec0[0], evec0[1]);
 633     Float_t circ=0;
 634     if(jev==1) circ=2*eval[0];
 635     if(jev==0) circ=2*eval[1];
 636     fh2Circularity->Fill(circ,ptbig);
 637     fh2JetEntries->Fill(ptbig,pthardbin);
 638
 639
 640     if(fDoMatching) fh2JetAxisPhi->Fill(smearphi,pthardbin);
 641
 642
 643
 644     for (Int_t ip = 0; ip < ParticleList.GetEntries(); ip++) {
 645       AliVParticle *track = (AliVParticle*)ParticleList.At(ip);
 646       Float_t isembed=0.5;
 647       if(fDoMatching) if(ip<nT) isembed=1.5;
 648       TVector3 pp(track->Px(), track->Py(), track->Pz());
 649       Float_t phi = track->Phi();
 650       Float_t eta = track->Eta();
 651       Float_t pt  = track->Pt();
 652
 653       Float_t deta = eta - etabig;
 654       Float_t dphi = RelativePhi(phi,phibig);
 655       if(TMath::Abs(dphi)>=0.5*TMath::Pi()) continue;
 656
 657       Float_t dRR = TMath::Sqrt(dphi * dphi + deta * deta);
 658       TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
 659       TVector3 pPerp = pp - pLong;
 660       Float_t ppjX = pPerp.Dot(ppJ2);
 661       Float_t ppjY = pPerp.Dot(ppJ3);
 662       TVector2 vr(ppjX, ppjY) ;
 663       //and this is the angle between the particle and the TM axis.
 664       float phistr=evec.Phi()-vr.Phi();
 665       if(phistr>2*TMath::Pi()) phistr -= 2*TMath::Pi();
 666       if(phistr<-2*TMath::Pi()) phistr += 2*TMath::Pi();
 667       if(phistr<-0.5*TMath::Pi()) phistr += 2*TMath::Pi();
 668       if(phistr>1.5*TMath::Pi()) phistr -= 2*TMath::Pi();
 669
 670       double jetEntries[10] = {centValue,dRR,ptbig,pt,phistr,circ,ptrue,pthardbin,ptmax,isembed};
 671       fhnJetTM->Fill(jetEntries);
 672
 673     } // 2nd Track loop
 674   }//jet loop
 675
 676   PostData(1, fOutputList);
 677 }
 678
 679 void AliAnalysisTaskJetAntenna::Terminate(const Option_t *)
 680 {
 681    // Draw result to the screen
 682    // Called once at the end of the query
 683
 684    if (!GetOutputData(1))
 685    return;
 686 }
 687
 688
 689 Int_t  AliAnalysisTaskJetAntenna::GetListOfTracks(TList *list){
 690
 691   Int_t iCount = 0;
 692   AliAODEvent *aod = 0;
 693
 694   if(!fESD)aod = fAODIn;
 695   else aod = fAODOut;
 696   if(!aod)return 0;
 697
 698
 699     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 700     AliAODTrack *tr = aod->GetTrack(it);
 701     Bool_t bGood = false;
 702     if(fFilterType == 0)bGood = true;
 703     else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
 704     else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
 705     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 706     if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
 707     if(bGood==false) continue;
 708     if (fApplySharedClusterCut) {
 709       Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
 710       if (frac > 0.4) continue;
 711     }
 712     if(TMath::Abs(tr->Eta())>0.9)continue;
 713     if(tr->Pt()<0.15)continue;
 714     list->Add(tr);
 715     iCount++;
 716     }
 717   return iCount;
 718 }
 719
 720
 721
 722
 723 Int_t  AliAnalysisTaskJetAntenna::GetListOfTracksExtra(TList *list){
 724
 725   Int_t iCount = 0;
 726   Int_t nEmbed=0;
 727   AliAODEvent *aod = 0;
 728
 729   if(!fESD)aod = fAODIn;
 730   else aod = fAODOut;
 731   if(!aod)return 0;
 732
 733
 734       TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
 735       if(!aodExtraTracks)return iCount;
 736       for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
 737         AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
 738         if (!track) continue;
 739         AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
 740         if(!trackAOD)continue;
 741         Bool_t bGood = false;
 742         if(fFilterType == 0)bGood = true;
 743         else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
 744         else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
 745         if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
 746         if(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
 747          if (fApplySharedClusterCut) {
 748            Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls());
 749            if (frac > 0.4) continue;
 750          }
 751
 752         if(TMath::Abs(trackAOD->Eta())>0.9) continue;
 753         if(trackAOD->Pt()<0.15) continue;
 754         if(fDebug) printf("pt extra track %.2f \n", trackAOD->Pt());
 755         list->Add(trackAOD);
 756         iCount++;
 757         }
 758
 759       nEmbed=iCount-1;
 760
 761
 762
 763
 764     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 765     AliAODTrack *tr = aod->GetTrack(it);
 766     Bool_t bGood = false;
 767     if(fFilterType == 0)bGood = true;
 768     else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
 769     else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
 770     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 771     if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
 772     if(bGood==false) continue;
 773     if (fApplySharedClusterCut) {
 774       Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
 775       if (frac > 0.4) continue;
 776     }
 777     if(TMath::Abs(tr->Eta())>0.9)continue;
 778     if(tr->Pt()<0.15)continue;
 779     list->Add(tr);
 780     iCount++;
 781     }
 782
 783      return nEmbed;
 784 }
 785
 786 Double_t AliAnalysisTaskJetAntenna::RelativePhi(Double_t mphi,Double_t vphi){
 787
 788   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
 789   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
 790   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
 791   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
 792   double dphi = mphi-vphi;
 793   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
 794   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
 795   return dphi;//dphi in [-Pi, Pi]
 796 }
 797
 798 Int_t AliAnalysisTaskJetAntenna::GetPhiBin(Double_t phi)
 799 {
 800     Int_t phibin=-1;
 801     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
 802     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
 803     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
 804     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
 805     return phibin;
 806 }
 807
 808 Int_t AliAnalysisTaskJetAntenna::GetPtHardBin(Double_t ptHard){
 809
 810    const Int_t nBins = 10;
 811   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
 812
 813   Int_t bin = -1;
 814   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
 815     bin++;
 816   }
 817
 818   return bin;
 819 }