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