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