PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx

   1 #include <TROOT.h>
   2 #include <TSystem.h>
   3 #include <TInterpreter.h>
   4 #include <TChain.h>
   5 #include <TFile.h>
   6 #include <TH1F.h>
   7 #include <TH1.h>
   8 #include <TH2.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TList.h>
  12 #include <TTree.h>
  13 #include <TBranch.h>
  14 #include <TLorentzVector.h>
  15 #include <TClonesArray.h>
  16 #include <TObjArray.h>
  17 #include <TRefArray.h>
  18 #include <TArrayD.h>
  19 #include <fstream>
  20 #include <TVector3.h>
  21 #include <TVectorD.h>
  22 #include <TMatrixDSym.h>
  23 #include <TMatrixDSymEigen.h>
  24 #include <TStyle.h>
  25 #include <TProfile.h>
  26 #include  "TDatabasePDG.h"
  27
  28 #include "AliAnalysisTaskThreeJets.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliJetFinder.h"
  31 #include "AliJetReader.h"
  32 #include "AliJetHeader.h"
  33 #include "AliJetReaderHeader.h"
  34 #include "AliUA1JetHeaderV1.h"
  35 #include "AliESDEvent.h"
  36 #include "AliAODEvent.h"
  37 #include "AliAODVertex.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODTrack.h"
  40 #include "AliAODJet.h"
  41 #include "AliMCEventHandler.h"
  42 #include "AliMCEvent.h"
  43 #include "AliStack.h"
  44 #include "AliGenPythiaEventHeader.h"
  45 #include "AliJetKineReaderHeader.h"
  46 #include "AliGenCocktailEventHeader.h"
  47 #include "AliAODPid.h"
  48 #include "AliExternalTrackParam.h"
  49
  50 #include "AliAnalysisTaskJetSpectrum.h"
  51
  52 #include "AliAnalysisHelperJetTasks.h"
  53
  54
  55 ClassImp(AliAnalysisTaskThreeJets)
  56
  57   AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets() : AliAnalysisTaskSE(),
  58
  59                                                          fAOD(0x0),
  60
  61                                                          fBranchRec(""),
  62                                                          fBranchGen(""),
  63
  64                                                          fUseAODInput(kFALSE),
  65
  66                                                          fR(0x0),
  67                                                          fList(0x0),
  68
  69                                                          fGlobVar(1),
  70                                                          fXsection(1),
  71
  72                                                          fhX3X4Rec(0x0),
  73                                                          fhX3X4Gen(0x0),
  74
  75                                                          fhMu35Rec(0x0),
  76                                                          fhMu34Rec(0x0),
  77                                                          fhMu45Rec(0x0),
  78
  79                                                          fhMu35Gen(0x0),
  80                                                          fhMu34Gen(0x0),
  81                                                          fhMu45Gen(0x0),
  82
  83                                                          fhInOut(0x0),
  84
  85                                                          fhThrustRec2(0x0),
  86                                                          fhThrustRec3(0x0),
  87
  88                                                          fhThrustGen2(0x0),
  89                                                          fhThrustGen3(0x0),
  90
  91                                                          fhCGen2(0x0),
  92                                                          fhCGen3(0x0),
  93
  94                                                          fhSGen2(0x0),
  95                                                          fhSGen3(0x0),
  96
  97                                                          fhAGen2(0x0),
  98                                                          fhAGen3(0x0),
  99
 100                                                          fhCRec2(0x0),
 101                                                          fhCRec3(0x0),
 102
 103                                                          fhSRec2(0x0),
 104                                                          fhSRec3(0x0),
 105
 106                                                          fhARec2(0x0),
 107                                                          fhARec3(0x0),
 108
 109                                                          fhX3(0x0),
 110                                                          fhX4(0x0),
 111                                                          fhX5(0x0),
 112
 113                                                          fhXSec(0x0),
 114
 115                                                          fhX3X4Rec60(0x0),
 116                                                          fhX3X4Rec60100(0x0),
 117                                                          fhX3X4Rec100(0x0),
 118
 119                                                          fhX3X4Gen60(0x0),
 120                                                          fhX3X4Gen60100(0x0),
 121                                                          fhX3X4Gen100(0x0),
 122
 123                                                          fhdPhiThrustGen(0x0),
 124                                                          fhdPhiThrustGenALL(0x0),
 125
 126                                                          fhdPhiThrustRec(0x0),
 127                                                          fhdPhiThrustRecALL(0x0)
 128 {
 129
 130 }
 131
 132 AliAnalysisTaskThreeJets::AliAnalysisTaskThreeJets(const char * name):
 133   AliAnalysisTaskSE(name),
 134
 135   fAOD(0x0),
 136
 137   fBranchRec(""),
 138   fBranchGen(""),
 139
 140   fUseAODInput(kFALSE),
 141
 142   fR(0x0),
 143   fList(0x0),
 144
 145   fGlobVar(1),
 146   fXsection(1),
 147
 148   fhX3X4Rec(0x0),
 149   fhX3X4Gen(0x0),
 150
 151   fhMu35Rec(0x0),
 152   fhMu34Rec(0x0),
 153   fhMu45Rec(0x0),
 154
 155   fhMu35Gen(0x0),
 156   fhMu34Gen(0x0),
 157   fhMu45Gen(0x0),
 158
 159   fhInOut(0x0),
 160
 161   fhThrustRec2(0x0),
 162   fhThrustRec3(0x0),
 163
 164   fhThrustGen2(0x0),
 165   fhThrustGen3(0x0),
 166
 167   fhCGen2(0x0),
 168   fhCGen3(0x0),
 169
 170   fhSGen2(0x0),
 171   fhSGen3(0x0),
 172
 173   fhAGen2(0x0),
 174   fhAGen3(0x0),
 175
 176   fhCRec2(0x0),
 177   fhCRec3(0x0),
 178
 179   fhSRec2(0x0),
 180   fhSRec3(0x0),
 181
 182   fhARec2(0x0),
 183   fhARec3(0x0),
 184
 185   fhX3(0x0),
 186   fhX4(0x0),
 187   fhX5(0x0),
 188
 189   fhXSec(0x0),
 190
 191   fhX3X4Rec60(0x0),
 192   fhX3X4Rec60100(0x0),
 193   fhX3X4Rec100(0x0),
 194
 195   fhX3X4Gen60(0x0),
 196   fhX3X4Gen60100(0x0),
 197   fhX3X4Gen100(0x0),
 198
 199   fhdPhiThrustGen(0x0),
 200   fhdPhiThrustGenALL(0x0),
 201
 202   fhdPhiThrustRec(0x0),
 203   fhdPhiThrustRecALL(0x0)
 204 {
 205   DefineOutput(1, TList::Class());
 206 }
 207
 208
 209
 210 Bool_t AliAnalysisTaskThreeJets::Notify()
 211 {
 212   //
 213   // Implemented Notify() to read the cross sections
 214   // and number of trials from pyxsec.root
 215   //
 216   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 217   UInt_t   ntrials  = 0;
 218   if(tree){
 219     TFile *curfile = tree->GetCurrentFile();
 220     if (!curfile) {
 221       Error("Notify","No current file");
 222       return kFALSE;
 223     }
 224
 225     TString fileName(curfile->GetName());
 226     if(fileName.Contains("AliESDs.root")){
 227         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
 228     }
 229     else if(fileName.Contains("AliAOD.root")){
 230         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
 231     }
 232     else if(fileName.Contains("galice.root")){
 233         // for running with galice and kinematics alone...
 234         fileName.ReplaceAll("galice.root", "pyxsec.root");
 235     }
 236     TFile *fxsec = TFile::Open(fileName.Data());
 237     if(!fxsec){
 238       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
 239       // no a severe condition
 240       return kTRUE;
 241     }
 242     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 243     if(!xtree){
 244       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
 245     }
 246     xtree->SetBranchAddress("xsection",&fXsection);
 247     xtree->SetBranchAddress("ntrials",&ntrials);
 248     xtree->GetEntry(0);
 249   }
 250
 251   return kTRUE;
 252 }
 253
 254
 255 //___________________________________________________________________________________________________________________________________
 256 void AliAnalysisTaskThreeJets::UserCreateOutputObjects()
 257 {
 258   //
 259   // Create the output container
 260   //
 261   Printf("Analysing event  %s :: # %5d\n", gSystem->pwd(), (Int_t) fEntry);
 262
 263   if(fUseAODInput){
 264     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 265     if(!fAOD){
 266       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
 267       return;
 268     }
 269   }
 270   else{
 271     //  assume that the AOD is in the general output...
 272     fAOD  = AODEvent();
 273     if(!fAOD){
 274       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
 275       return;
 276     }
 277   }
 278
 279   printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
 280
 281   fList = new TList();
 282
 283   fhX3X4Gen = new TH2F("X3vsX4Gen", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 284   fhX3X4Gen->SetXTitle("X_{3}");
 285   fhX3X4Gen->SetYTitle("X_{4}");
 286   fhX3X4Gen->Sumw2();
 287   fList->Add(fhX3X4Gen);
 288
 289   fhX3X4Rec = new TH2F("X3vsX4Rec", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 290   fhX3X4Rec->SetXTitle("X_{3}");
 291   fhX3X4Rec->SetYTitle("X_{4}");
 292   fhX3X4Rec->Sumw2();
 293   fList->Add(fhX3X4Rec);
 294
 295   fhMu35Rec = new TH1F("Mu35Rec", "", 20,0.1,0.8);
 296   fhMu35Rec->Sumw2();
 297   fhMu35Rec->SetXTitle("#mu_{35}");
 298   fhMu35Rec->SetYTitle("#frac{dN}{d#mu_{35Rec}}");
 299   fList->Add(fhMu35Rec);
 300
 301   fhMu34Rec = new TH1F("Mu34Rec", "", 20,0.5,1);
 302   fhMu34Rec->Sumw2();
 303   fhMu34Rec->SetXTitle("#mu_{34}");
 304   fhMu34Rec->SetYTitle("#frac{dN}{d#mu_{34}}");
 305   fList->Add(fhMu34Rec);
 306
 307   fhMu45Rec = new TH1F("Mu45Rec", "", 20,0,0.65);
 308   fhMu45Rec->Sumw2();
 309   fhMu45Rec->SetXTitle("#mu_{45}");
 310   fhMu45Rec->SetYTitle("#frac{dN}{d#mu_{45}}");
 311   fList->Add(fhMu45Rec);
 312
 313   fhMu35Gen = new TH1F("Mu35Gen", "", 20,0.1,0.8);
 314   fhMu35Gen->Sumw2();
 315   fhMu35Gen->SetXTitle("#mu_{35Gen}");
 316   fhMu35Gen->SetYTitle("#frac{dN}{d#mu_{35Gen}}");
 317   fList->Add(fhMu35Gen);
 318
 319   fhMu34Gen = new TH1F("Mu34Gen", "", 20,0.5,1);
 320   fhMu34Gen->Sumw2();
 321   fhMu34Gen->SetXTitle("#mu_{34Gen}");
 322   fhMu34Gen->SetYTitle("#frac{dN}{d#mu_{34Gen}}");
 323   fList->Add(fhMu34Gen);
 324
 325   fhMu45Gen = new TH1F("Mu45Gen", "", 20,0,0.65);
 326   fhMu45Gen->Sumw2();
 327   fhMu45Gen->SetXTitle("#mu_{45Gen}");
 328   fhMu45Gen->SetYTitle("#frac{dN}{d#mu_{45}}");
 329   fList->Add(fhMu45Gen);
 330
 331   fhInOut = new TH1I("InOut", "", 6, 0, 6);
 332   fhInOut->SetXTitle("#RecJets_{GenJets=3}");
 333   fhInOut->SetYTitle("#Entries");
 334   fList->Add(fhInOut);
 335
 336   fhThrustGen2 = new TH1F("ThrustGen2", "", 50, 0.5, 1);
 337   fhThrustGen2->Sumw2();
 338   fList->Add(fhThrustGen2);
 339
 340   fhThrustGen3 = new TH1F("ThrustGen3", "", 50, 0.5, 1);
 341   fhThrustGen3->Sumw2();
 342   fList->Add(fhThrustGen3);
 343
 344   fhThrustRec2 = new TH1F("ThrustRec2", "", 50, 0.5, 1);
 345   fhThrustRec2->Sumw2();
 346   fList->Add(fhThrustRec2);
 347
 348   fhThrustRec3 = new TH1F("ThrustRec3", "", 50, 0.5, 1);
 349   fhThrustRec3->Sumw2();
 350   fList->Add(fhThrustRec3);
 351
 352   fhCGen2 = new TH1F("CGen2", "", 100, 0, 1);
 353   fhCGen2->Sumw2();
 354   fList->Add(fhCGen2);
 355
 356   fhCGen3 = new TH1F("CGen3", "", 100, 0, 1);
 357   fhCGen3->Sumw2();
 358   fList->Add(fhCGen3);
 359
 360   fhCRec2 = new TH1F("CRec2", "", 100, 0, 1);
 361   fhCRec2->Sumw2();
 362   fList->Add(fhCRec2);
 363
 364   fhCRec3 = new TH1F("CRec3", "", 100, 0, 1);
 365   fhCRec3->Sumw2();
 366   fList->Add(fhCRec3);
 367
 368   fhSGen2 = new TH1F("SGen2", "", 100, 0, 1);
 369   fList->Add(fhSGen2);
 370
 371   fhSGen3 = new TH1F("SGen3", "", 100, 0, 1);
 372   fList->Add(fhSGen3);
 373
 374   fhSRec2 = new TH1F("SRec2", "", 100, 0, 1);
 375   fList->Add(fhSRec2);
 376
 377   fhSRec3 = new TH1F("SRec3", "", 100, 0, 1);
 378   fList->Add(fhSRec3);
 379
 380   fhAGen2 = new TH1F("AGen2", "", 50, 0, 0.5);
 381   fList->Add(fhAGen2);
 382
 383   fhAGen3 = new TH1F("AGen3", "", 50, 0, 0.5);
 384   fList->Add(fhAGen3);
 385
 386   fhARec2 = new TH1F("ARec2", "", 50, 0, 0.5);
 387   fList->Add(fhARec2);
 388
 389   fhARec3 = new TH1F("ARec3", "", 50, 0, 0.5);
 390   fList->Add(fhARec3);
 391
 392   fhX3 = new TH2F("X3", "", 22, 0.6, 1.02, 100, 0, 1);
 393   fhX3->SetYTitle("|X_{3}^{MC} - X_{3}^{AOD}|/X_{3}^{MC}");
 394   fhX3->SetXTitle("X_{3}");
 395   fhX3->Sumw2();
 396   fList->Add(fhX3);
 397
 398   fhX4 = new TH2F("X4", "",33, 0.4, 1.02, 100, 0, 1);
 399   fhX4->SetYTitle("|X_{4}^{MC} - X_{4}^{AOD}|/X_{4}^{MC}");
 400   fhX4->SetXTitle("X_{4}");
 401   fhX4->Sumw2();
 402   fList->Add(fhX4);
 403
 404   fhX5 = new TH2F("X5", "",100, 0., 1., 100, 0, 1);
 405   fhX5->SetYTitle("|X_{5}^{MC} - X_{5}^{AOD}|/X_{5}^{MC}");
 406   fhX5->SetXTitle("X_{5}");
 407   fhX5->Sumw2();
 408   fList->Add(fhX5);
 409
 410   fhXSec = new TProfile("XSec", "", 200, 0, 200, 0, 1);
 411   fList->Add(fhXSec);
 412
 413   fhX3X4Rec60 = new TH2F("X3vsX4Rec60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 414   fhX3X4Rec60->SetXTitle("X_{3}");
 415   fhX3X4Rec60->SetYTitle("X_{4}");
 416   fhX3X4Rec60->Sumw2();
 417   fList->Add(fhX3X4Rec60);
 418
 419   fhX3X4Rec60100 = new TH2F("X3vsX4Rec60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 420   fhX3X4Rec60100->SetXTitle("X_{3}");
 421   fhX3X4Rec60100->SetYTitle("X_{4}");
 422   fhX3X4Rec60100->Sumw2();
 423   fList->Add(fhX3X4Rec60100);
 424
 425   fhX3X4Rec100 = new TH2F("X3vsX4Rec100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 426   fhX3X4Rec100->SetXTitle("X_{3}");
 427   fhX3X4Rec100->SetYTitle("X_{4}");
 428   fhX3X4Rec100->Sumw2();
 429   fList->Add(fhX3X4Rec100);
 430
 431   fhX3X4Gen60 = new TH2F("X3vsX4Gen60", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 432   fhX3X4Gen60->SetXTitle("X_{3}");
 433   fhX3X4Gen60->SetYTitle("X_{4}");
 434   fhX3X4Gen60->Sumw2();
 435   fList->Add(fhX3X4Gen60);
 436
 437   fhX3X4Gen60100 = new TH2F("X3vsX4Gen60100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 438   fhX3X4Gen60100->SetXTitle("X_{3}");
 439   fhX3X4Gen60100->SetYTitle("X_{4}");
 440   fhX3X4Gen60100->Sumw2();
 441   fList->Add(fhX3X4Gen60100);
 442
 443   fhX3X4Gen100 = new TH2F("X3vsX4Gen100", "", 22, 0.6, 1.02, 33, 0.4, 1.02);
 444   fhX3X4Gen100->SetXTitle("X_{3}");
 445   fhX3X4Gen100->SetYTitle("X_{4}");
 446   fhX3X4Gen100->Sumw2();
 447   fList->Add(fhX3X4Gen100);
 448
 449   fhdPhiThrustGen = new TH2F("dPhiThrustGen", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
 450   fhdPhiThrustGen->Sumw2();
 451   fList->Add(fhdPhiThrustGen);
 452
 453   fhdPhiThrustGenALL = new TH2F("dPhiThrustGenALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0,  150);
 454   fhdPhiThrustGenALL->Sumw2();
 455   fList->Add(fhdPhiThrustGenALL);
 456
 457   fhdPhiThrustRec = new TH2F("dPhiThrustRec", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
 458   fhdPhiThrustRec->Sumw2();
 459   fList->Add(fhdPhiThrustRec);
 460
 461   fhdPhiThrustRecALL = new TH2F("dPhiThrustRecALL", "", 32, -1*TMath::Pi(), TMath::Pi(), 25, 0, 150);
 462   fhdPhiThrustRecALL->Sumw2();
 463   fList->Add(fhdPhiThrustRecALL);
 464
 465   Printf("UserCreateOutputObjects finished\n");
 466 }
 467
 468 //__________________________________________________________________________________________________________________________________________
 469 void AliAnalysisTaskThreeJets::Init()
 470 {
 471   printf("AliAnalysisJetCut::Init() \n");
 472 }
 473
 474 //____________________________________________________________________________________________________________________________________________
 475 void AliAnalysisTaskThreeJets::UserExec(Option_t * )
 476 {
 477   //  if (fDebug > 1) printf("Analysing event # %5d\n", (Int_t) fEntry);
 478
 479
 480   //create an AOD handler
 481   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 482
 483   if(!aodH)
 484     {
 485       Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
 486       return;
 487     }
 488
 489   AliMCEvent* mcEvent =MCEvent();
 490   if(!mcEvent){
 491     Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
 492     return;
 493   }
 494
 495   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
 496
 497   //primary vertex
 498   AliAODVertex * pvtx = dynamic_cast<AliAODVertex*>(fAOD->GetPrimaryVertex());
 499   if(!pvtx) return;
 500
 501   AliAODJet genJetsPythia[kMaxJets];
 502   Int_t nPythiaGenJets = 0;
 503
 504   AliAODJet genJets[kMaxJets];
 505   Int_t nGenJets = 0;
 506
 507   AliAODJet recJets[kMaxJets];
 508   Int_t nRecJets = 0;
 509
 510   //array of reconstructed jets from the AOD input
 511   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
 512   if(!aodRecJets){
 513     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
 514     return;
 515   }
 516
 517   // reconstructed jets
 518   nRecJets = aodRecJets->GetEntries();
 519   nRecJets = TMath::Min(nRecJets, kMaxJets);
 520
 521   for(int ir = 0;ir < nRecJets;++ir)
 522     {
 523      AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
 524      if(!tmp)continue;
 525      recJets[ir] = *tmp;
 526     }
 527
 528   // If we set a second branch for the input jets fetch this
 529   TClonesArray * aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
 530
 531   if(!aodGenJets)
 532     {
 533       printf("NO MC jets Found\n");
 534       return;
 535     }
 536
 537   //   //Generated jets
 538   nGenJets = aodGenJets->GetEntries();
 539   nGenJets = TMath::Min(nGenJets, kMaxJets);
 540
 541   for(Int_t ig =0 ; ig < nGenJets; ++ig)
 542     {
 543       AliAODJet * tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
 544       if(!tmp)continue;
 545       genJets[ig] = * tmp;
 546     }
 547
 548   AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 549   if(!pythiaGenHeader){
 550     Printf("!!!NO GEN HEADER AVALABLE!!!");
 551     return;
 552   }
 553
 554   // Int_t ProcessType = pythiaGenHeader->ProcessType();
 555   // if(ProcessType != 28) return;
 556   fhXSec->Fill(pythiaGenHeader->GetPtHard(), fXsection);
 557   nPythiaGenJets = pythiaGenHeader->NTriggerJets();
 558   nPythiaGenJets = TMath::Min(nPythiaGenJets, kMaxJets);
 559
 560   Double_t eRec[kMaxJets];
 561   Double_t eGen[kMaxJets];
 562
 563   Double_t eJetRec[kMaxJets];
 564   //  Double_t EJetGen[kMaxJets];
 565
 566   AliAODJet jetRec[kMaxJets];
 567   AliAODJet jetGen[kMaxJets];
 568
 569   Int_t idxRec[kMaxJets];
 570   Int_t idxGen[kMaxJets];
 571
 572   Double_t xRec[kMaxJets];
 573   Double_t xGen[kMaxJets];
 574
 575   Double_t eSumRec = 0;
 576   Double_t eSumGen = 0;
 577
 578   TLorentzVector vRec[kMaxJets];
 579   TLorentzVector vRestRec[kMaxJets];
 580
 581   TLorentzVector vGen[kMaxJets];
 582   TLorentzVector vRestGen[kMaxJets];
 583
 584   TLorentzVector vsumRec;
 585   TLorentzVector vsumGen;
 586
 587   TVector3 pRec[kMaxJets];
 588   TVector3 pGen[kMaxJets];
 589
 590   TVector3 pTrack[kTracks];
 591
 592   TVector3 pRestRec[kMaxJets];
 593   TVector3 pRestGen[kMaxJets];
 594
 595   Double_t psumRestRec = 0;
 596   //  Double_t psumRestGen = 0;
 597   //Pythia_________________________________________________________________________________________________________________
 598
 599   for(int ip = 0;ip < nPythiaGenJets;++ip)
 600     {
 601       if(ip>=kMaxJets)continue;
 602       Float_t p[4];
 603       pythiaGenHeader->TriggerJet(ip,p);
 604       genJetsPythia[ip].SetPxPyPzE(p[0],p[1],p[2],p[3]);
 605     }
 606
 607 //_________________________________________________________________________________________________________________________
 608
 609
 610 //________histos for MC___________________________________________________________________________________________________________
 611
 612   Int_t nGenSel = 0;
 613   Int_t counter = 0;
 614   Int_t tag = 0;
 615
 616   AliAODJet selJets[kMaxJets];
 617
 618   for(Int_t i = 0; i < nGenJets; i++)
 619     {
 620       if(nGenJets == 1)
 621         {
 622           selJets[nGenSel] = genJets[i];
 623           nGenSel++;
 624         }
 625       else
 626         {
 627           counter = 0;
 628           tag = 0;
 629           for(Int_t j = 0; j < nGenJets; j++)
 630             {
 631               if(i!=j)
 632                 {
 633                   Double_t dRij = genJets[i].DeltaR(&genJets[j]);
 634                   counter++;
 635                   if(dRij > 2*fR) tag++;
 636                 }
 637             }
 638           if(counter!=0)
 639             {
 640               if(tag/counter == 1)
 641                 {
 642                   selJets[nGenSel] = genJets[i];
 643                   nGenSel++;
 644                 }
 645             }
 646         }
 647     }
 648
 649   if(nGenSel == 0) return;
 650
 651   for (Int_t gj = 0; gj < nGenSel; gj++)
 652     {
 653       eGen[gj] = selJets[gj].E();
 654     }
 655
 656   TMath::Sort(nGenSel, eGen, idxGen);
 657   for (Int_t ig = 0; ig < nGenSel; ig++)
 658     {
 659       jetGen[ig] = selJets[idxGen[ig]];
 660     }
 661   AliStack * stack = mcEvent->Stack();
 662   Int_t nMCtracks = stack->GetNprimary();
 663   Double_t * eTracksMC = new Double_t[kTracks];
 664   Double_t pTracksMC[kTracks];
 665   Int_t * idxTracksMC = new Int_t[kTracks];
 666   TLorentzVector jetTracksMC[kTracks];
 667   TLorentzVector jetTracksSortMC[kTracks];
 668   TVector3 pTrackMC[kTracks];
 669   TLorentzVector vTrackMCAll[kTracks];
 670   Double_t pTrackMCAll[kTracks];
 671   TLorentzVector vTrackMC[kTracks];
 672   TVector3 pTrackMCBoost[kTracks];
 673   Double_t eventShapes[4];
 674
 675   Int_t nAccTr = 0;
 676   Int_t nInJet[kMaxJets];
 677   TLorentzVector inJetPartV[kMaxJets][kTracks];
 678   Int_t nAllTracksMC = 0;
 679
 680   for(Int_t iTrack = 0; iTrack < nMCtracks; iTrack++)
 681     {
 682       TParticle * part = (TParticle*)stack->Particle(iTrack);
 683       if (!part) continue;
 684       Double_t fEta = part->Eta();
 685       if(TMath::Abs(fEta) > .9) continue;
 686       if(!IsPrimChar(part, nMCtracks, 0)) continue;
 687       vTrackMCAll[nAllTracksMC].SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 688       pTrackMCAll[nAllTracksMC] = part->Pt();
 689       nAllTracksMC++;
 690     }
 691   if(nAllTracksMC == 0) return;
 692   for(Int_t iJet = 0; iJet < nGenSel; iJet++)
 693     {
 694       Int_t nJetTracks = 0;
 695       for(Int_t i = 0; i < nAllTracksMC; i++)
 696         {
 697           Double_t dPhi = (jetGen[iJet].Phi()-vTrackMCAll[i].Phi());
 698           if(dPhi > TMath::Pi()) dPhi = dPhi - 2.*TMath::Pi();
 699           if(dPhi < (-1.*TMath::Pi())) dPhi = dPhi + 2.*TMath::Pi();
 700           Double_t dEta = (jetGen[iJet].Eta()-vTrackMCAll[i].Eta());
 701           Double_t deltaR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 702           if(deltaR < fR && vTrackMCAll[i].Pt() > 1.5)
 703             {
 704               jetTracksMC[nAccTr] = vTrackMCAll[i];
 705               eTracksMC[nAccTr] = vTrackMCAll[i].E();
 706               pTracksMC[nAccTr] = vTrackMCAll[i].Pt();
 707               inJetPartV[iJet][nJetTracks].SetPxPyPzE(vTrackMCAll[i].Px(), vTrackMCAll[i].Py(), vTrackMCAll[i].Pz(),vTrackMCAll[i].E());
 708               nAccTr++;
 709               nJetTracks++;
 710             }
 711         }
 712       nInJet[iJet] = nJetTracks;
 713     }
 714
 715   if(nAccTr == 0) return;
 716   Printf("*********** Number of Jets : %d ***************\n", nGenSel);
 717   Double_t pTav[kMaxJets];
 718   for(Int_t i = 0; i < nGenSel; i++)
 719     {
 720       Double_t pTsum = 0;
 721       Printf("*********** Number of particles in Jet %d = %d *******************\n", i+3, nInJet[i]);
 722       for(Int_t iT = 0; iT < nInJet[i]; iT++)
 723         {
 724           Double_t pt = inJetPartV[i][iT].Pt();
 725           pTsum += pt;
 726         }
 727       pTav[i] = pTsum/nInJet[i];
 728     }
 729
 730   TMath::Sort(nAccTr, pTracksMC, idxTracksMC);
 731   for(Int_t i = 0; i < nAccTr; i++)
 732     {
 733       jetTracksSortMC[i] = jetTracksMC[idxTracksMC[i]];
 734       pTrackMC[i].SetXYZ(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz());
 735       vTrackMC[i].SetPxPyPzE(jetTracksSortMC[i].Px(), jetTracksSortMC[i].Py(), jetTracksSortMC[i].Pz(), jetTracksSortMC[i].E());
 736     }
 737
 738   TVector3 n01MC = pTrackMC[0].Unit();
 739   //Thrust calculation, iterative method
 740   if(nGenSel > 1)
 741     {
 742 //       if(fGlobVar == 1)
 743 //      {
 744       GetEventShapes(n01MC, pTrackMC, nAccTr, eventShapes);
 745 //      }
 746
 747       Double_t s = eventShapes[1];
 748       Double_t a = eventShapes[2];
 749       Double_t c = eventShapes[3];
 750
 751       switch(nGenSel)
 752         {
 753         case 2:
 754           {
 755             fhAGen2->Fill(a);
 756             fhSGen2->Fill(s);
 757             fhCGen2->Fill(c);
 758           }
 759           break;
 760         case 3:
 761           {
 762             fhAGen3->Fill(a);
 763             fhSGen3->Fill(s);
 764             fhCGen3->Fill(c);
 765           }
 766           break;
 767         }
 768       Double_t thrust01MC = eventShapes[0];
 769
 770       switch(nGenSel)
 771         {
 772         case 2:
 773           fhThrustGen2->Fill(thrust01MC, fXsection);
 774           break;
 775         case 3:
 776           fhThrustGen3->Fill(thrust01MC, fXsection);
 777           break;
 778         }
 779     }
 780
 781
 782   //rest frame MC jets
 783   for (Int_t i = 0; i < nGenSel; ++i)
 784     {
 785       vGen[i].SetPxPyPzE(jetGen[i].Px(), jetGen[i].Py(), jetGen[i].Pz(), jetGen[i].E());
 786       pGen[i].SetXYZ(vGen[i].Px(), vGen[i].Py(), vGen[i].Pz());
 787       vsumGen += vGen[i];
 788     }
 789   if(eventShapes[0] >0.8 )
 790     {
 791       for(Int_t i = 0; i < nGenSel; i++)
 792         fhdPhiThrustGen->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
 793     }
 794
 795   if(eventShapes[0] <= 0.8)
 796     {
 797       for(Int_t i = 0; i < nGenSel; i++)
 798         fhdPhiThrustGenALL->Fill(n01MC.DeltaPhi(pGen[i]), jetGen[i].E());
 799     }
 800
 801   Double_t fPxGen = vsumGen.Px();
 802   Double_t fPyGen = vsumGen.Py();
 803   Double_t fPzGen = vsumGen.Pz();
 804   Double_t fEGen = vsumGen.E();
 805
 806   Double_t eRestGen[kMaxJets];
 807   for (Int_t j = 0; j < nGenSel; j++)
 808     {
 809       vGen[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
 810       eRestGen[j] = vGen[j].E();
 811     }
 812
 813   for (Int_t j = 0; j < nAccTr; j++)
 814     {
 815       vTrackMC[j].Boost(-fPxGen/fEGen, -fPyGen/fEGen, -fPzGen/fEGen);
 816       pTrackMCBoost[j].SetXYZ(vTrackMC[j].Px(),vTrackMC[j].Py(),vTrackMC[j].Pz());
 817     }
 818
 819   Int_t idxRestGen[kMaxJets];
 820   TMath::Sort(nGenSel, eRestGen, idxRestGen);
 821   for(Int_t j = 0; j < nGenSel; j++)
 822     {
 823       vRestGen[j] = vGen[idxRestGen[j]];
 824       eSumGen += vRestGen[j].E();
 825     }
 826
 827   if (nGenSel == 3)
 828     {
 829       //      if(nInJet[0] < 3 || nInJet[1] < 3 || nInJet[2] < 3) return;
 830       //       if(pRestGen[1].DeltaPhi(pRestGen[2]) > 0.95 && pRestGen[1].DeltaPhi(pRestGen[2]) < 1.15)
 831       //        {
 832
 833       for(Int_t i = 0; i < nGenSel; i++)
 834         {
 835           xGen[i] = 2*vRestGen[i].E()/eSumGen;
 836         }
 837
 838       Printf("***************** Values of Dalitz variables are : %f, %f, %f ****************\n", xGen[0], xGen[1], xGen[2]);
 839
 840       Printf("***************** fXSection = %f ******************\n", fXsection);
 841       if(eSumGen <= 60)
 842         fhX3X4Gen60->Fill(xGen[0], xGen[1], fXsection);
 843
 844       if(eSumGen > 60 && eSumGen <= 100)
 845         fhX3X4Gen60100->Fill(xGen[0], xGen[1], fXsection);
 846
 847       if(eSumGen > 100)
 848         fhX3X4Gen100->Fill(xGen[0], xGen[1], fXsection);
 849
 850       FillTopology(fhX3X4Gen, fhMu34Gen, fhMu45Gen, fhMu35Gen, xGen, pRestGen, fXsection);
 851     }
 852
 853
 854
 855 //_______________________________________________histos for MC_____________________________________________________
 856
 857
 858 //_______________________________________________histos AOD________________________________________________________
 859
 860 //   Printf("Event Number : %d, Number of gen jets : %d ", fEntry, nGenJets);
 861
 862   Int_t nRecSel = 0;
 863   Int_t counter1 = 0;
 864   Int_t tag1 = 0;
 865
 866   AliAODJet recSelJets[kMaxJets];
 867
 868   for(Int_t i = 0; i < nRecJets; i++)
 869     {
 870       if(nRecJets == 1)
 871         {
 872           recSelJets[nRecSel] = recJets[i];
 873           nRecSel++;
 874         }
 875       else
 876         {
 877           counter1 = 0;
 878           tag1 = 0;
 879           for(Int_t j = 0; j < nRecJets; j++)
 880             {
 881               if(i!=j)
 882                 {
 883                   Double_t dRij = recJets[i].DeltaR(&recJets[j]);
 884                   counter1++;
 885                   if(dRij > 2*fR) tag1++;
 886                 }
 887             }
 888           if(counter1!=0)
 889             {
 890               if(tag1/counter1 == 1)
 891                 {
 892                   recSelJets[nRecSel] = recJets[i];
 893                   nRecSel++;
 894                 }
 895             }
 896         }
 897     }
 898
 899   if(nRecSel == 0) return;
 900
 901   //sort rec/gen jets by energy in C.M.S
 902   for (Int_t rj = 0; rj < nRecSel; rj++)
 903     {
 904       eRec[rj] = recSelJets[rj].E();
 905     }
 906
 907   //  Int_t nAODtracks = fAOD->GetNumberOfTracks();
 908   Int_t nTracks = 0; //tracks accepted in the whole event
 909   //  Int_t nTracksALL = 0;
 910   TLorentzVector jetTracks[kTracks];
 911   TLorentzVector jetTracksSort[kTracks];
 912   Double_t * eTracks = new Double_t[kTracks];
 913   Double_t pTracks[kTracks];
 914   Int_t * idxTracks = new Int_t[kTracks];
 915   Double_t eventShapesRec[4];
 916   Int_t jetMult[kMaxJets];
 917   //  TLorentzVector vTracksAll[kTracks];
 918   //  Double_t pTracksAll[kTracks];
 919   Int_t nAccJets = 0;
 920   AliAODJet jetRecAcc[kMaxJets];
 921   Int_t nJetTracks = 0;
 922
 923   AliAODTrack jetTrack[kTracks];
 924   Double_t * cv = new Double_t[21];
 925   TMath::Sort(nRecSel, eRec, idxRec);
 926   for (Int_t rj = 0; rj < nRecSel; rj++)
 927     {
 928       nJetTracks = 0;
 929       eJetRec[rj] = eRec[idxRec[rj]];
 930       jetRec[rj] = recSelJets[idxRec[rj]];
 931       TRefArray * jetTracksAOD = dynamic_cast<TRefArray*>(jetRec[rj].GetRefTracks());
 932       if(!jetTracksAOD) continue;
 933       if(jetTracksAOD->GetEntries() < 3) continue;
 934       for(Int_t i = 0; i < jetTracksAOD->GetEntries(); i++)
 935         {
 936           AliAODTrack * track = (AliAODTrack*)jetTracksAOD->At(i);
 937           track->GetCovarianceXYZPxPyPz(cv);
 938           if(cv[14] > 1000.) continue;
 939           jetTrack[nTracks] = *track;
 940           jetTracks[nTracks].SetPxPyPzE(jetTrack[nTracks].Px(), jetTrack[nTracks].Py(), jetTrack[nTracks].Pz(), jetTrack[nTracks].E());
 941           eTracks[nTracks] = jetTracks[nTracks].E();
 942           pTracks[nTracks] = jetTracks[nTracks].Pt();
 943           nTracks++;
 944           nJetTracks++;
 945         }
 946       if(nJetTracks < 3) continue;
 947       jetRecAcc[nAccJets] = jetRec[rj];
 948       jetMult[nAccJets] = jetTracksAOD->GetEntries();
 949       nAccJets++;
 950     }
 951
 952   if (nAccJets == 0) return;
 953
 954   TLorentzVector vTrack[kTracks];
 955   TMath::Sort(nTracks, pTracks, idxTracks);
 956   for(Int_t i = 0; i < nTracks; i++)
 957     {
 958       jetTracksSort[i] = jetTracks[idxTracks[i]];
 959       pTrack[i].SetXYZ(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz());
 960       vTrack[i].SetPxPyPzE(jetTracksSort[i].Px(), jetTracksSort[i].Py(), jetTracksSort[i].Pz(), jetTracksSort[i].E());
 961     }
 962
 963   for (Int_t i = 0; i < nAccJets; ++i)
 964     {
 965       vRec[i].SetPxPyPzE(jetRecAcc[i].Px(), jetRecAcc[i].Py(), jetRecAcc[i].Pz(), jetRecAcc[i].E());
 966       pRec[i].SetXYZ(vRec[i].Px(), vRec[i].Py(), vRec[i].Pz());
 967       vsumRec += vRec[i];
 968     }
 969
 970   //Thrust, iterative method, AODs
 971   TVector3 n01 = pTrack[0].Unit();
 972   if(nAccJets > 1)
 973     {
 974 //       if(fGlobVar == 1)
 975 //      {
 976       GetEventShapes(n01, pTrack, nTracks, eventShapesRec);
 977 //      }
 978 //       fGlobVar = 0;
 979 //       Double_t Max3 = TMath::Max(eventShapesRec0[0],eventShapesRec1[0]);
 980 //       Double_t Max4 = TMath::Max(eventShapesRec3[0],eventShapesRec2[0]);
 981
 982       Double_t thrust = eventShapesRec[0];
 983
 984       if(eventShapesRec[0] > 0.8)
 985         {
 986           for(Int_t i = 0; i < nAccJets; i++)
 987             fhdPhiThrustRec->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
 988
 989         }
 990
 991       if(eventShapesRec[0] <= 0.8)
 992         {
 993           for(Int_t i = 0; i < nAccJets; i++)
 994             fhdPhiThrustRecALL->Fill(n01.DeltaPhi(pRec[i]), jetRecAcc[i].E());
 995         }
 996
 997       switch(nAccJets)
 998         {
 999         case 2:
1000           fhThrustRec2->Fill(thrust, fXsection);
1001           break;
1002         case 3:
1003           fhThrustRec3->Fill(thrust, fXsection);
1004           break;
1005         }
1006
1007       switch(nAccJets)
1008         {
1009         case 2:
1010           {
1011             fhARec2->Fill(eventShapesRec[2], fXsection);
1012             fhSRec2->Fill(eventShapesRec[1], fXsection);
1013             fhCRec2->Fill(eventShapesRec[3], fXsection);
1014           }
1015           break;
1016         case 3:
1017           {
1018             fhARec3->Fill(eventShapesRec[2], fXsection);
1019             fhSRec3->Fill(eventShapesRec[1], fXsection);
1020             fhCRec3->Fill(eventShapesRec[3], fXsection);
1021           }
1022           break;
1023         }
1024
1025     }
1026
1027   //rest frame for reconstructed jets
1028   Double_t fPx = vsumRec.Px();
1029   Double_t fPy = vsumRec.Py();
1030   Double_t fPz = vsumRec.Pz();
1031   Double_t fE = vsumRec.E();
1032
1033   TVector3 pTrackRest[kTracks];
1034   for(Int_t j = 0; j < nTracks; j++)
1035     {
1036       vTrack[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1037       pTrackRest[j].SetXYZ(vTrack[j].Px(), vTrack[j].Py(),vTrack[j].Pz());
1038     }
1039
1040   Double_t eRestRec[kMaxJets];
1041   Int_t idxRestRec[kMaxJets];
1042   for (Int_t j = 0; j < nAccJets; j++)
1043     {
1044       vRec[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1045       eRestRec[j] = vRec[j].E();
1046       eSumRec += vRec[j].E();
1047     }
1048
1049   TMath::Sort(nAccJets, eRestRec, idxRestRec);
1050   for (Int_t i = 0; i < nAccJets; i++)
1051     {
1052       vRestRec[i] = vRec[idxRestRec[i]];
1053       pRestRec[i].SetXYZ(vRestRec[i].Px(), vRestRec[i].Py(), vRestRec[i].Pz());
1054       psumRestRec += pRestRec[i].Perp();
1055     }
1056
1057   if(nAccJets == 3)
1058     {
1059 //       if(pRest[1].DeltaPhi(pRest[2]) > 0.95 && pRest[1].DeltaPhi(pRest[2]) < 1.15)
1060 //      {
1061           fhInOut->Fill(nGenSel);
1062 //        for(Int_t j = 0; j < nTracksALL; j++)
1063 //          {
1064 //            vTracksAll[j].Boost(-fPx/fE, -fPy/fE, -fPz/fE);
1065 //            pTracksAll[j].SetXYZ(vTracksAll[j].Px(), vTracksAll[j].Py(),vTracksAll[j].Pz());
1066 //            fhdPhiRec->Fill(pRest[0].DeltaPhi(pTracksAll[j]), pTracksAll[j].Perp(), fXsection);
1067 //          }
1068           //and the Dalitz variables and Energy distributions in the rest frame
1069           for (Int_t i = 0; i < nAccJets; i++)
1070             xRec[i] = 2*vRestRec[i].E()/eSumRec;
1071
1072           if(eSumRec <= 60)
1073             fhX3X4Rec60->Fill(xRec[0], xRec[1], fXsection);
1074
1075           if(eSumRec > 60 && eSumRec <= 100)
1076             fhX3X4Rec60100->Fill(xRec[0], xRec[1], fXsection);
1077
1078           if(eSumRec > 100)
1079             fhX3X4Rec100->Fill(xRec[0], xRec[1], fXsection);
1080
1081           if(nAccJets == 3 && nAccJets == nGenJets)
1082             {
1083               fhX3->Fill(xGen[0], TMath::Abs(xGen[0]-xRec[0])/xGen[0], fXsection);
1084               fhX4->Fill(xGen[1], TMath::Abs(xGen[1]-xRec[1])/xGen[1], fXsection);
1085               fhX5->Fill(xGen[2], TMath::Abs(xGen[2]-xRec[2])/xGen[2], fXsection);
1086             }
1087
1088           FillTopology(fhX3X4Rec, fhMu34Rec, fhMu45Rec, fhMu35Rec, xRec, pRestRec, fXsection);
1089     }
1090   Printf("%s:%d",(char*)__FILE__,__LINE__);
1091
1092   PostData(1, fList);
1093
1094   Printf("%s:%d Data Posted",(char*)__FILE__,__LINE__);
1095
1096 }
1097
1098 //__________________________________________________________________________________________________________________________________________________
1099 void AliAnalysisTaskThreeJets::Terminate(Option_t *)
1100 {
1101   printf("AnalysisJetCorrelation::Terminate()");
1102
1103 }
1104
1105 //_______________________________________User defined functions_____________________________________________________________________________________
1106 void AliAnalysisTaskThreeJets::FillTopology(TH2F * Dalitz, TH1F * fhMu34, TH1F * fhMu45, TH1F * fhMu35, Double_t * x, TVector3 * pRest, Double_t xsection)
1107 {
1108   Dalitz->Fill(x[0], x[1], xsection);
1109   fhMu35->Fill(TMath::Sqrt(x[0]*x[2]*(1-(pRest[0].Unit()).Dot(pRest[2].Unit()))/2), xsection);
1110   fhMu34->Fill(TMath::Sqrt(x[0]*x[1]*(1-(pRest[0].Unit()).Dot(pRest[1].Unit()))/2), xsection);
1111   fhMu45->Fill(TMath::Sqrt(x[1]*x[2]*(1-(pRest[1].Unit()).Dot(pRest[2].Unit()))/2), xsection);
1112 }
1113
1114 //_____________________________________________________________________________________________________________________________
1115
1116 Bool_t AliAnalysisTaskThreeJets::IsPrimChar(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug)
1117 {
1118   //
1119   // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack)
1120   // shall be counted as a primary particle
1121   //
1122   // This function or a equivalent should be available in some common place of AliRoot
1123   //
1124   // WARNING: Call this function only for particles that are among the particles from the event generator!
1125   // --> stack->Particle(id) with id < stack->GetNprimary()
1126
1127   // if the particle has a daughter primary, we do not want to count it
1128   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
1129   {
1130     if (adebug)
1131       printf("Dropping particle because it has a daughter among the primaries.\n");
1132     return kFALSE;
1133   }
1134
1135   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
1136
1137
1138   // skip quarks and gluon
1139   if (pdgCode <= 10 || pdgCode == 21)
1140   {
1141     if (adebug)
1142       printf("Dropping particle because it is a quark or gluon.\n");
1143     return kFALSE;
1144   }
1145
1146   Int_t status = aParticle->GetStatusCode();
1147   // skip non final state particles..
1148   if(status!=1){
1149     if (adebug)
1150       printf("Dropping particle because it is not a final state particle.\n");
1151     return kFALSE;
1152   }
1153
1154   if (strcmp(aParticle->GetName(),"XXX") == 0)
1155   {
1156     Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode);
1157     return kFALSE;
1158   }
1159
1160   TParticlePDG* pdgPart = aParticle->GetPDG();
1161
1162   if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
1163   {
1164     Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode);
1165     return kFALSE;
1166   }
1167
1168   if (pdgPart->Charge() == 0)
1169   {
1170     if (adebug)
1171       printf("Dropping particle because it is not charged.\n");
1172     return kFALSE;
1173   }
1174
1175   return kTRUE;
1176 }
1177
1178 //______________________________________________________________________________________________________
1179
1180 void AliAnalysisTaskThreeJets::GetThrustAxis(TVector3 &n01, TVector3 * pTrack, Int_t &nTracks)
1181 {
1182   TVector3 psum;
1183   Double_t psum1 = 0;
1184   Double_t psum2 = 0;
1185   Double_t thrust[kTracks];
1186   Double_t th = -3;
1187
1188   for(Int_t j = 0; j < nTracks; j++)
1189     {
1190       psum.SetXYZ(0., 0., 0.);
1191       psum1 = 0;
1192       psum2 = 0;
1193       for(Int_t i = 0; i < nTracks; i++)
1194         {
1195           psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
1196           psum2 += pTrack[i].Mag();
1197
1198           if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
1199           if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
1200         }
1201
1202       thrust[j] = psum1/psum2;
1203
1204       if(th == thrust[j])
1205         break;
1206
1207       th = thrust[j];
1208
1209       n01 = psum.Unit();
1210     }
1211 }
1212
1213 //___________________________________________________________________________________________________________
1214
1215 void AliAnalysisTaskThreeJets::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
1216 {
1217   TVector3 psum;
1218   Double_t psum1 = 0;
1219   Double_t psum2 = 0;
1220   Double_t thrust[kTracks];
1221   Double_t th = -3;
1222
1223   //Sphericity calculation
1224   TMatrixDSym m(3);
1225   Double_t s00 = 0;
1226   Double_t s01 = 0;
1227   Double_t s02 = 0;
1228
1229   Double_t s10 = 0;
1230   Double_t s11 = 0;
1231   Double_t s12 = 0;
1232
1233   Double_t s20 = 0;
1234   Double_t s21 = 0;
1235   Double_t s22 = 0;
1236
1237   Double_t ptot = 0;
1238
1239   Double_t c = -10;
1240
1241   for(Int_t j = 0; j < nTracks; j++)
1242     {
1243       psum.SetXYZ(0., 0., 0.);
1244       psum1 = 0;
1245       psum2 = 0;
1246       for(Int_t i = 0; i < nTracks; i++)
1247         {
1248           psum1 += (TMath::Abs(n01.Dot(pTrack[i])));
1249           psum2 += pTrack[i].Mag();
1250
1251           if (n01.Dot(pTrack[i]) > 0) psum += pTrack[i];
1252           if (n01.Dot(pTrack[i]) < 0) psum -= pTrack[i];
1253         }
1254
1255       thrust[j] = psum1/psum2;
1256       if(th == thrust[j])
1257         break;
1258
1259       th = thrust[j];
1260
1261       n01 = psum.Unit();
1262     }
1263
1264   eventShapes[0] = th;
1265
1266   for(Int_t j = 0; j < nTracks; j++)
1267     {
1268       s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
1269       s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
1270       s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
1271
1272       s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
1273       s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
1274       s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
1275
1276       s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
1277       s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
1278       s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
1279
1280       ptot += pTrack[j].Mag();
1281     }
1282
1283   if(ptot !=0)
1284     {
1285       m(0,0) = s00/ptot;
1286       m(0,1) = s01/ptot;
1287       m(0,2) = s02/ptot;
1288
1289       m(1,0) = s10/ptot;
1290       m(1,1) = s11/ptot;
1291       m(1,2) = s12/ptot;
1292
1293       m(2,0) = s20/ptot;
1294       m(2,1) = s21/ptot;
1295       m(2,2) = s22/ptot;
1296
1297       TMatrixDSymEigen eigen(m);
1298       TVectorD eigenVal = eigen.GetEigenValues();
1299
1300       Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
1301       eventShapes[1] = sphericity;
1302
1303       Double_t aplaarity = (3/2)*(eigenVal(2));
1304       eventShapes[2] = aplaarity;
1305
1306       c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
1307       eventShapes[3] = c;
1308     }
1309 }
1310
1311 //__________________________________________________________________________________________________________________________
1312
1313 Double_t AliAnalysisTaskThreeJets::Exponent(Double_t x, Double_t * par)
1314 {
1315   return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2])+0.5*TMath::Power((x-par[3])/par[0], 2))+par[4]*x;
1316 }
1317
1318 Double_t AliAnalysisTaskThreeJets::Exponent2(Double_t x, Double_t * par)
1319 {
1320   return par[0]*TMath::Power(1/TMath::E(), TMath::Power(par[1]/x, par[2]))+par[3]*x;
1321 }
1322
1323 Double_t AliAnalysisTaskThreeJets::Gauss(Double_t x, Double_t * par)
1324 {
1325   return 1/(par[1])*TMath::Power(1/TMath::E(), 0.5*(x-par[0])*(x-par[0])/(par[1]*par[1]));
1326 }
1327
1328 Double_t AliAnalysisTaskThreeJets::Total(Double_t x, Double_t * par)
1329 {
1330   return Exponent(x, par)+Gauss(x, &par[4]);
1331 }
1332
1333