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