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