PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx

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