PWG4/JetTasks/AliAnalysisTaskThreeJets.cxx

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