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