PWG4/JetTasks/AliAnalysisHelperJetTasks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //
  17 // Helper Class that contains a lot of
  18 // usefull static functions jet matchin pythia access etc.
  19 //
  20 // Author: C. Klein-Boesing IKP Muenster
  21
  22 #include "TROOT.h"
  23 #include "TDirectory.h"
  24 #include "TKey.h"
  25 #include "TList.h"
  26 #include "TSystem.h"
  27 #include "TH1F.h"
  28 #include "TProfile.h"
  29 #include "THnSparse.h"
  30 #include "TSeqCollection.h"
  31 #include "TMethodCall.h"
  32 #include "TFile.h"
  33 #include "TString.h"
  34 #include "TArrayI.h"
  35 #include "TArrayS.h"
  36
  37 #include "AliMCEvent.h"
  38 #include "AliLog.h"
  39 #include "AliESDEvent.h"
  40 #include "AliAODJet.h"
  41 #include "AliAODEvent.h"
  42 #include "AliStack.h"
  43 #include "AliGenEventHeader.h"
  44 #include "AliGenCocktailEventHeader.h"
  45 #include <AliGenDPMjetEventHeader.h>
  46 #include "AliGenPythiaEventHeader.h"
  47 #include <fstream>
  48 #include <iostream>
  49 #include "AliAnalysisHelperJetTasks.h"
  50 #include "TMatrixDSym.h"
  51 #include "TMatrixDSymEigen.h"
  52 #include "TVector.h"
  53
  54 ClassImp(AliAnalysisHelperJetTasks)
  55
  56 Int_t AliAnalysisHelperJetTasks::fgLastProcessType = -1;
  57
  58
  59 AliGenPythiaEventHeader*  AliAnalysisHelperJetTasks::GetPythiaEventHeader(const AliMCEvent *mcEvent){
  60   //
  61   // Fetch  the pythia event header
  62   //
  63   if(!mcEvent)return 0;
  64   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
  65   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
  66   if(!pythiaGenHeader){
  67     // cocktail ??
  68     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
  69
  70     if (!genCocktailHeader) {
  71       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
  72       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
  73       return 0;
  74     }
  75     TList* headerList = genCocktailHeader->GetHeaders();
  76     for (Int_t i=0; i<headerList->GetEntries(); i++) {
  77       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
  78       if (pythiaGenHeader)
  79         break;
  80     }
  81     if(!pythiaGenHeader){
  82       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
  83       return 0;
  84     }
  85   }
  86   return pythiaGenHeader;
  87
  88 }
  89
  90
  91 void AliAnalysisHelperJetTasks::PrintStack(AliMCEvent *mcEvent,Int_t iFirst,Int_t iLast,Int_t iMaxPrint){
  92   //
  93   // Print the stack informatuin up to the iMaxPrint event
  94   //
  95
  96   AliStack *stack = mcEvent->Stack();
  97   if(!stack){
  98     Printf("%s%d No Stack available",(char*)__FILE__,__LINE__);
  99     return;
 100   }
 101
 102   static Int_t iCount = 0;
 103   if(iCount>iMaxPrint)return;
 104   Int_t nStack = stack->GetNtrack();
 105   if(iLast == 0)iLast = nStack;
 106   else if(iLast > nStack)iLast = nStack;
 107
 108
 109   Printf("####################################################################");
 110   for(Int_t np = iFirst;np<iLast;++np){
 111     TParticle *p = stack->Particle(np);
 112     Printf("Nr.%d --- Status %d ---- Mother1 %d Mother2 %d Daughter1 %d Daughter2 %d ",
 113            np,p->GetStatusCode(),p->GetMother(0),p->GetMother(1),p->GetDaughter(0),p->GetDaughter(1));
 114     Printf("Eta %3.3f Phi %3.3f ",p->Eta(),p->Phi());
 115     p->Print();
 116     Printf("---------------------------------------");
 117   }
 118   iCount++;
 119 }
 120
 121
 122
 123
 124 void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &kGenJets,
 125                                                const AliAODJet *recJets,const Int_t &kRecJets,
 126                                                Int_t *iGenIndex,Int_t *iRecIndex,
 127                                                Int_t iDebug,Float_t maxDist){
 128
 129   //
 130   // Relate the two input jet Arrays
 131   //
 132
 133   //
 134   // The association has to be unique
 135   // So check in two directions
 136   // find the closest rec to a gen
 137   // and check if there is no other rec which is closer
 138   // Caveat: Close low energy/split jets may disturb this correlation
 139
 140
 141   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
 142   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
 143   // in the end we have something like this
 144   //    1   2   3
 145   // ------------
 146   // a| 3   2   0
 147   // b| 0   1   0
 148   // c| 0   0   3
 149   // d| 0   0   1
 150   // e| 0   0   1
 151   // Topology
 152   //   1     2
 153   //     a         b
 154   //
 155   //  d      c
 156   //        3     e
 157   // Only entries with "3" match from both sides
 158
 159   // In case we have more jets than kmaxjets only the
 160   // first kmaxjets are searched
 161   // all other are -1
 162   // use kMaxJets for a test not to fragemnt the memory...
 163
 164   for(int i = 0;i < kGenJets;++i)iGenIndex[i] = -1;
 165   for(int j = 0;j < kRecJets;++j)iRecIndex[j] = -1;
 166
 167
 168
 169   const int kMode = 3;
 170
 171   const Int_t nGenJets = TMath::Min(kMaxJets,kGenJets);
 172   const Int_t nRecJets = TMath::Min(kMaxJets,kRecJets);
 173
 174   if(nRecJets==0||nGenJets==0)return;
 175
 176   // UShort_t *iFlag = new UShort_t[nGenJets*nRecJets];
 177   UShort_t iFlag[kMaxJets*kMaxJets] = {0}; // all values set to zero
 178   for(int i = 0;i < nGenJets;++i){
 179     for(int j = 0;j < nRecJets;++j){
 180       iFlag[i*nGenJets+j] = 0;
 181     }
 182   }
 183
 184
 185
 186   // find the closest distance to the generated
 187   for(int ig = 0;ig<nGenJets;++ig){
 188     Float_t dist = maxDist;
 189     if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
 190     for(int ir = 0;ir<nRecJets;++ir){
 191       if(iDebug>1){
 192         printf("Rec (%d) ",ir);
 193         Printf("p_T %3.3f eta %3.3f ph %3.3f ",recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
 194       }
 195       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 196       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
 197       if(dR<dist){
 198         iRecIndex[ig] = ir;
 199         dist = dR;
 200       }
 201     }
 202     if(iRecIndex[ig]>=0)iFlag[ig*nRecJets+iRecIndex[ig]]+=1;
 203     // reset...
 204     iRecIndex[ig] = -1;
 205   }
 206   // other way around
 207   for(int ir = 0;ir<nRecJets;++ir){
 208     Float_t dist = maxDist;
 209     for(int ig = 0;ig<nGenJets;++ig){
 210       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
 211       if(dR<dist){
 212         iGenIndex[ir] = ig;
 213         dist = dR;
 214       }
 215     }
 216     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nRecJets+ir]+=2;
 217     // reset...
 218     iGenIndex[ir] = -1;
 219   }
 220
 221   // check for "true" correlations
 222
 223   if(iDebug>1)Printf(">>>>>> Matrix");
 224
 225   for(int ig = 0;ig<nGenJets;++ig){
 226     for(int ir = 0;ir<nRecJets;++ir){
 227       // Print
 228       if(iDebug>1)printf("Flag[%d][%d] %d ",ig,ir,iFlag[ig*nRecJets+ir]);
 229
 230       if(kMode==3){
 231         // we have a uniqie correlation
 232         if(iFlag[ig*nRecJets+ir]==3){
 233           iGenIndex[ir] = ig;
 234           iRecIndex[ig] = ir;
 235         }
 236       }
 237       else{
 238         // we just take the correlation from on side
 239         if((iFlag[ig*nRecJets+ir]&2)==2){
 240           iGenIndex[ir] = ig;
 241         }
 242         if((iFlag[ig*nRecJets+ir]&1)==1){
 243           iRecIndex[ig] = ir;
 244         }
 245       }
 246     }
 247     if(iDebug>1)printf("\n");
 248   }
 249 }
 250
 251
 252
 253 void AliAnalysisHelperJetTasks::GetClosestJets(const TList *genJetsList,const Int_t &kGenJets,
 254                                                const TList *recJetsList,const Int_t &kRecJets,
 255                                                TArrayI &iGenIndex,TArrayI &iRecIndex,
 256                                                Int_t iDebug,Float_t maxDist){
 257
 258   // Size indepnedendentt Implemenation of jet matching
 259   // Thepassed TArrayI should be static in the user function an only increased if needed
 260
 261   //
 262   // Relate the two input jet Arrays
 263   //
 264
 265   //
 266   // The association has to be unique
 267   // So check in two directions
 268   // find the closest rec to a gen
 269   // and check if there is no other rec which is closer
 270   // Caveat: Close low energy/split jets may disturb this correlation
 271
 272
 273   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
 274   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
 275   // in the end we have something like this
 276   //    1   2   3
 277   // ------------
 278   // a| 3   2   0
 279   // b| 0   1   0
 280   // c| 0   0   3
 281   // d| 0   0   1
 282   // e| 0   0   1
 283   // Topology
 284   //   1     2
 285   //     a         b
 286   //
 287   //  d      c
 288   //        3     e
 289   // Only entries with "3" match from both sides
 290
 291   iGenIndex.Reset(-1);
 292   iRecIndex.Reset(-1);
 293
 294   const int kMode = 3;
 295   const Int_t nGenJets = TMath::Min(genJetsList->GetEntries(),kGenJets);
 296   const Int_t nRecJets = TMath::Min(recJetsList->GetEntries(),kRecJets);
 297   if(nRecJets==0||nGenJets==0)return;
 298
 299   static TArrayS iFlag(nGenJets*nRecJets);
 300   if(iFlag.GetSize()<(nGenJets*nRecJets)){
 301     iFlag.Set(nGenJets*nRecJets);
 302   }
 303   iFlag.Reset(0);
 304
 305   // find the closest distance to the generated
 306   for(int ig = 0;ig<nGenJets;++ig){
 307     AliAODJet *genJet = (AliAODJet*)genJetsList->At(ig);
 308     if(!genJet)continue;
 309
 310     Float_t dist = maxDist;
 311     if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJet->Pt(),genJet->Eta(),genJet->Phi());
 312     for(int ir = 0;ir<nRecJets;++ir){
 313       AliAODJet *recJet = (AliAODJet*)recJetsList->At(ir);
 314       if(!recJet)continue;
 315       if(iDebug>1){
 316         printf("Rec (%d) ",ir);
 317         Printf("p_T %3.3f eta %3.3f ph %3.3f ",recJet->Pt(),recJet->Eta(),recJet->Phi());
 318       }
 319       Double_t dR = genJet->DeltaR(recJet);
 320       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
 321       if(dR<dist){
 322         iRecIndex[ig] = ir;
 323         dist = dR;
 324       }
 325     }
 326     if(iRecIndex[ig]>=0)iFlag[ig*nRecJets+iRecIndex[ig]]+=1;
 327     // reset...
 328     iRecIndex[ig] = -1;
 329   }
 330   // other way around
 331
 332   for(int ir = 0;ir<nRecJets;++ir){
 333       AliAODJet *recJet = (AliAODJet*)recJetsList->At(ir);
 334       if(!recJet)continue;
 335       Float_t dist = maxDist;
 336       for(int ig = 0;ig<nGenJets;++ig){
 337         AliAODJet *genJet = (AliAODJet*)genJetsList->At(ig);
 338         if(!genJet)continue;
 339         Double_t dR = genJet->DeltaR(recJet);
 340         if(dR<dist){
 341         iGenIndex[ir] = ig;
 342         dist = dR;
 343       }
 344     }
 345     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nRecJets+ir]+=2;
 346     // reset...
 347     iGenIndex[ir] = -1;
 348   }
 349
 350   // check for "true" correlations
 351
 352   if(iDebug>1)Printf(">>>>>> Matrix Size %d",iFlag.GetSize());
 353
 354   for(int ig = 0;ig<nGenJets;++ig){
 355     for(int ir = 0;ir<nRecJets;++ir){
 356       // Print
 357       if(iDebug>1)printf("Flag2[%d][%d] %d ",ig,ir,iFlag[ig*nRecJets+ir]);
 358
 359       if(kMode==3){
 360         // we have a uniqie correlation
 361         if(iFlag[ig*nRecJets+ir]==3){
 362           iGenIndex[ir] = ig;
 363           iRecIndex[ig] = ir;
 364         }
 365       }
 366       else{
 367         // we just take the correlation from on side
 368         if((iFlag[ig*nRecJets+ir]&2)==2){
 369           iGenIndex[ir] = ig;
 370         }
 371         if((iFlag[ig*nRecJets+ir]&1)==1){
 372           iRecIndex[ig] = ir;
 373         }
 374       }
 375     }
 376     if(iDebug>1)printf("\n");
 377   }
 378 }
 379
 380 void AliAnalysisHelperJetTasks::GetJetMatching(const TList *genJetsList, const Int_t &kGenJets,
 381                                                const TList *recJetsList, const Int_t &kRecJets,
 382                                                TArrayI &iMatchIndex, TArrayF &fPtFraction,
 383                                                Int_t iDebug, Float_t maxDist){
 384
 385
 386     // Matching jets from two lists
 387     // Start with highest energetic jet from first list (generated/embedded)
 388     // Calculate distance (\Delta R) to all jets from second list (reconstructed)
 389     // Select N closest jets = kClosestJetsN
 390     // Check energy fraction from jets from first list in jets from second list
 391     // Matched jets = jet with largest energy fraction
 392     // Store index of matched jet in TArrayI iMatchIndex
 393
 394     // reset index
 395     iMatchIndex.Reset(-1);
 396     fPtFraction.Reset(-1.);
 397
 398     // N closest jets: store list with index and \Delta R
 399     const Int_t kClosestJetsN = 4;
 400     Double_t closestJets[kClosestJetsN][2]; //[][0] = index, [][1] = \Delta R
 401
 402     for(Int_t i=0; i<kClosestJetsN; ++i){
 403         for(Int_t j=0; j<2; ++j){
 404             closestJets[i][j] = 1e6;
 405         }
 406     }
 407
 408
 409     const Int_t nGenJets = TMath::Min(genJetsList->GetEntries(),kGenJets);
 410     const Int_t nRecJets = TMath::Min(recJetsList->GetEntries(),kRecJets);
 411     if(nRecJets==0||nGenJets==0) return;
 412
 413     AliAODJet *genJet = 0x0;
 414     AliAODJet *recJet = 0x0;
 415
 416     // loop over generated/embedded jets
 417     for(Int_t ig=0; ig<nGenJets; ++ig){
 418         genJet = (AliAODJet*)genJetsList->At(ig);
 419         //if(!genJet || !JetSelected(genJet)) continue;
 420         if(!genJet) continue;
 421
 422         // find N closest reconstructed jets
 423         Double_t deltaR = 0.;
 424         for(Int_t ir=0; ir<nRecJets; ++ir){
 425             recJet = (AliAODJet*)recJetsList->At(ir);
 426             //if(!recJet || !JetSelected(recJet)) continue;
 427             if(!recJet) continue;
 428
 429             deltaR = genJet->DeltaR(recJet);
 430
 431             Int_t i=kClosestJetsN-1;
 432             if(deltaR<closestJets[i][1] && deltaR<maxDist){
 433                 closestJets[i][0] = (Double_t) ir; // index
 434                 closestJets[i][1] = deltaR;
 435
 436                 // sort array (closest at the beginning)
 437                 while(i>=1 && closestJets[i][1]<closestJets[i-1][1]){
 438                     Double_t tmpArr[2];
 439                     for(Int_t j=0; j<2; j++){
 440                        tmpArr[j] = closestJets[i-1][j];
 441                        closestJets[i-1][j]   = closestJets[i][j];
 442                        closestJets[i][j] = tmpArr[j];
 443                     }
 444                     i--;
 445                 }
 446             }
 447         } // end: loop over reconstructed jets
 448
 449         // calculate fraction for the N closest jets
 450         Double_t maxFraction = 0.; // maximum found fraction in one jets
 451         Double_t cumFraction = 0.; // cummulated fraction of closest jets (for break condition)
 452         Double_t fraction = 0.;
 453         Int_t ir = -1;  // index of close reconstruced jet
 454
 455         for(Int_t irc=0; irc<kClosestJetsN; irc++){
 456             ir = (Int_t)(closestJets[irc][0]);
 457             recJet = (AliAODJet*)recJetsList->At(ir);
 458             if(!(recJet)) continue;
 459
 460             fraction = GetFractionOfJet(recJet,genJet);
 461
 462             cumFraction += fraction;
 463
 464             // check if jet fullfills current matching condition
 465             if(fraction>maxFraction){
 466                 // avoid multiple links
 467                 for(Int_t ij=0; ij<ig; ++ij){
 468                     if(iMatchIndex[ij]==ir) continue;
 469                 }
 470                 // set index
 471                 maxFraction = fraction;
 472                 fPtFraction[ig] = fraction;
 473                 iMatchIndex[ig] = ir;
 474             }
 475             // break condition: less energy left as already found in one jet or
 476             // as required for positiv matching
 477             if(1-cumFraction<maxFraction) break;
 478         } // end: loop over closest jets
 479
 480         if(iMatchIndex[ig]<0){
 481             if(iDebug) Printf("Matching failed for (gen) jet #%d", ig);
 482         }
 483     }
 484 }
 485
 486 Double_t AliAnalysisHelperJetTasks::GetFractionOfJet(const AliAODJet *recJet, const AliAODJet *genJet){
 487   //
 488   // get the fraction of hte signal jet in the full jt
 489   //
 490     Double_t ptGen = genJet->Pt();
 491     if(ptGen==0.) return 999.;
 492
 493     Double_t ptAssocTracks = 0.; // sum of pT of tracks found in both jets
 494
 495     // look at tracks inside jet
 496     TRefArray *genTrackList = genJet->GetRefTracks();
 497     TRefArray *recTrackList = recJet->GetRefTracks();
 498     Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 499     Int_t nTracksRecJet = recTrackList->GetEntriesFast();
 500
 501     AliAODTrack* recTrack;
 502     AliAODTrack* genTrack;
 503     for(Int_t ir=0; ir<nTracksRecJet; ++ir){
 504         recTrack = (AliAODTrack*)(recTrackList->At(ir));
 505         if(!recTrack) continue;
 506
 507         for(Int_t ig=0; ig<nTracksGenJet; ++ig){
 508             genTrack = (AliAODTrack*)(genTrackList->At(ig));
 509             if(!genTrack) continue;
 510
 511             // look if it points to the same track
 512             if(genTrack==recTrack){
 513                 ptAssocTracks += genTrack->Pt();
 514                 break;
 515             }
 516         }
 517     }
 518
 519     // calculate fraction
 520     Double_t fraction = ptAssocTracks/ptGen;
 521
 522     return fraction;
 523 }
 524
 525
 526 void  AliAnalysisHelperJetTasks::MergeOutputDirs(const char* cFiles,const char* cPattern,const char *cOutFile,Bool_t bUpdate){
 527   // Routine to merge only directories containing the pattern
 528   //
 529   TString outFile(cOutFile);
 530   if(outFile.Length()==0)outFile = Form("%s.root",cPattern);
 531   ifstream in1;
 532   in1.open(cFiles);
 533   // open all files
 534   TList fileList;
 535   char cFile[240];
 536   while(in1>>cFile){// only open the first file
 537     Printf("Adding %s",cFile);
 538     TFile *f1 = TFile::Open(cFile);
 539     fileList.Add(f1);
 540   }
 541
 542   TFile *fOut = 0;
 543   if(fileList.GetEntries()){// open the first file
 544     TFile* fIn = dynamic_cast<TFile*>(fileList.At(0));
 545     if(!fIn){
 546       Printf("Input File not Found");
 547       return;
 548     }
 549     // fetch the keys for the directories
 550     TList *ldKeys = fIn->GetListOfKeys();
 551     for(int id = 0;id < ldKeys->GetEntries();id++){
 552       // check if it is a directory
 553       TKey *dKey = (TKey*)ldKeys->At(id);
 554       TDirectory *dir = dynamic_cast<TDirectory*>(dKey->ReadObj());
 555       if(dir){
 556         TString dName(dir->GetName());
 557         if(dName.Contains(cPattern)){
 558           // open new file if first match
 559           if(!fOut){
 560             if(bUpdate)fOut = new TFile(outFile.Data(),"UPDATE");
 561             else fOut = new TFile(outFile.Data(),"RECREATE");
 562           }
 563           // merge all the stuff that is in this directory
 564           TList *llKeys = dir->GetListOfKeys();
 565           TList *tmplist;
 566           TMethodCall callEnv;
 567
 568           fOut->cd();
 569           TDirectory *dOut = fOut->mkdir(dir->GetName());
 570
 571           for(int il = 0;il < llKeys->GetEntries();il++){
 572             TKey *lKey = (TKey*)llKeys->At(il);
 573             TObject *object = dynamic_cast<TObject*>(lKey->ReadObj());
 574             //  from TSeqCollections::Merge
 575             if(!object)continue;
 576             // If current object is not mergeable just skip it
 577             if (!object->IsA()) {
 578               continue;
 579             }
 580             callEnv.InitWithPrototype(object->IsA(), "Merge", "TCollection*");
 581             if (!callEnv.IsValid()) {
 582               continue;
 583             }
 584             // Current object mergeable - get corresponding objects in input lists
 585             tmplist = new TList();
 586             for(int i = 1; i < fileList.GetEntries();i++){
 587               TDirectory *fInTmp = dynamic_cast<TDirectory*>(fileList.At(i));
 588               if(!fInTmp){
 589                 Printf("File %d not found",i);
 590                 continue;
 591               }
 592               TDirectory *dTmp = (TDirectory*)fInTmp->Get(dName.Data());
 593               if(!dTmp){
 594                 Printf("Directory %s not found",dName.Data());
 595                 continue;
 596               }
 597               TObject* oTmp = (TObject*)dTmp->Get(object->GetName());
 598               if(!oTmp){
 599                 Printf("Object %s not found",object->GetName());
 600                 continue;
 601               }
 602               tmplist->Add(oTmp);
 603             }
 604             callEnv.SetParam((Long_t) tmplist);
 605             callEnv.Execute(object);
 606             delete tmplist;tmplist = 0;
 607             dOut->cd();
 608             object->Write(object->GetName(),TObject::kSingleKey);
 609           }
 610         }
 611       }
 612     }
 613     if(fOut){
 614       fOut->Close();
 615     }
 616   }
 617 }
 618
 619 void  AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cDir, char *cList,char *cOutFile,Bool_t bUpdate){
 620
 621   // This is used to merge the analysis-output from different
 622   // data samples/pt_hard bins
 623   // in case the eventweigth was set to xsection/ntrials already, this
 624   // is not needed. Both methods only work in case we do not mix different
 625   // pt_hard bins, and do not have overlapping bins
 626
 627   const Int_t nMaxBins = 12;
 628   // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
 629   // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
 630   // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
 631   // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
 632
 633   Float_t xsection[nMaxBins];
 634   Float_t nTrials[nMaxBins];
 635   Float_t sf[nMaxBins];
 636   TList *lIn[nMaxBins];
 637   TDirectory *dIn[nMaxBins];
 638   TFile *fIn[nMaxBins];
 639
 640   ifstream in1;
 641   in1.open(cFiles);
 642
 643   char cFile[120];
 644   Int_t ibTotal = 0;
 645   while(in1>>cFile){
 646     fIn[ibTotal] = TFile::Open(cFile);
 647     if(strlen(cDir)==0){
 648       dIn[ibTotal] = gDirectory;
 649     }
 650     else{
 651       dIn[ibTotal] = (TDirectory*)fIn[ibTotal]->Get(cDir);
 652     }
 653     if(!dIn[ibTotal]){
 654       Printf("%s:%d No directory %s found, exiting...",__FILE__,__LINE__,cDir);
 655       fIn[ibTotal]->ls();
 656       return;
 657     }
 658
 659     lIn[ibTotal] = (TList*)dIn[ibTotal]->Get(cList);
 660     Printf("Merging file %s %s",cFile, cDir);
 661     if(!lIn[ibTotal]){
 662       Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList);
 663       fIn[ibTotal]->ls();
 664       return;
 665     }
 666     TH1* hTrials = (TH1F*)lIn[ibTotal]->FindObject("fh1Trials");
 667     if(!hTrials){
 668       Printf("%s:%d fh1PtHard_Trials not found in list, exiting...",__FILE__,__LINE__);
 669       return;
 670     }
 671     TProfile* hXsec = (TProfile*)lIn[ibTotal]->FindObject("fh1Xsec");
 672     if(!hXsec){
 673       Printf("%s:%d fh1Xsec  not found in list, exiting...",__FILE__,__LINE__);
 674       return;
 675     }
 676     xsection[ibTotal] = hXsec->GetBinContent(1);
 677     nTrials[ibTotal] = hTrials->Integral();
 678     sf[ibTotal] = xsection[ibTotal]/ nTrials[ibTotal];
 679     ibTotal++;
 680   }
 681
 682   if(ibTotal==0){
 683     Printf("%s:%d No files found for mergin, exiting",__FILE__,__LINE__);
 684     return;
 685   }
 686
 687   TFile *fOut = 0;
 688   if(bUpdate)fOut = new TFile(cOutFile,"UPDATE");
 689   else fOut = new TFile(cOutFile,"RECREATE");
 690   TDirectory *dOut = fOut->mkdir(dIn[0]->GetName());
 691   dOut->cd();
 692   TList *lOut = new TList();
 693   lOut->SetName(lIn[0]->GetName());
 694
 695   // for the start scale all...
 696   for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
 697     TH1 *h1Add = 0;
 698     THnSparse *hnAdd = 0;
 699     for(int ib = 0;ib < ibTotal;++ib){
 700       // dynamic cast does not work with cint
 701       TObject *h = lIn[ib]->At(ie);
 702       if(h->InheritsFrom("TH1")){
 703         TH1 *h1 = (TH1*)h;
 704         if(ib==0){
 705           h1Add = (TH1*)h1->Clone(h1->GetName());
 706           h1Add->Scale(sf[ib]);
 707         }
 708         else{
 709           h1Add->Add(h1,sf[ib]);
 710         }
 711       }
 712       else if(h->InheritsFrom("THnSparse")){
 713         THnSparse *hn = (THnSparse*)h;
 714         if(ib==0){
 715           hnAdd = (THnSparse*)hn->Clone(hn->GetName());
 716           hnAdd->Scale(sf[ib]);
 717         }
 718         else{
 719           hnAdd->Add(hn,sf[ib]);
 720         }
 721       }
 722
 723
 724     }// ib
 725     if(h1Add)lOut->Add(h1Add);
 726     else if(hnAdd)lOut->Add(hnAdd);
 727   }
 728   dOut->cd();
 729   lOut->Write(lOut->GetName(),TObject::kSingleKey);
 730   fOut->Close();
 731 }
 732
 733 Bool_t AliAnalysisHelperJetTasks::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 734   //
 735   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 736   // This is to called in Notify and should provide the path to the AOD/ESD file
 737
 738   TString file(currFile);
 739   fXsec = 0;
 740   fTrials = 1;
 741
 742   if(file.Contains("root_archive.zip#")){
 743     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 744     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 745     file.Replace(pos+1,20,"");
 746   }
 747   else {
 748     // not an archive take the basename....
 749     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 750   }
 751   Printf("%s",file.Data());
 752
 753
 754
 755
 756   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 757   if(!fxsec){
 758     // next trial fetch the histgram file
 759     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 760     if(!fxsec){
 761         // not a severe condition but inciate that we have no information
 762       return kFALSE;
 763     }
 764     else{
 765       // find the tlist we want to be independtent of the name so use the Tkey
 766       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 767       if(!key){
 768         fxsec->Close();
 769         return kFALSE;
 770       }
 771       TList *list = dynamic_cast<TList*>(key->ReadObj());
 772       if(!list){
 773         fxsec->Close();
 774         return kFALSE;
 775       }
 776       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 777       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 778       fxsec->Close();
 779     }
 780   } // no tree pyxsec.root
 781   else {
 782     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 783     if(!xtree){
 784       fxsec->Close();
 785       return kFALSE;
 786     }
 787     UInt_t   ntrials  = 0;
 788     Double_t  xsection  = 0;
 789     xtree->SetBranchAddress("xsection",&xsection);
 790     xtree->SetBranchAddress("ntrials",&ntrials);
 791     xtree->GetEntry(0);
 792     fTrials = ntrials;
 793     fXsec = xsection;
 794     fxsec->Close();
 795   }
 796   return kTRUE;
 797 }
 798
 799 Bool_t AliAnalysisHelperJetTasks::PrintDirectorySize(const char* currFile,Int_t iDetail){
 800
 801   //
 802   // Print the size on disk and in memory occuppied by a directory
 803   //
 804
 805   TFile *fDummy = 0;
 806   if(iDetail>=0)fDummy = TFile::Open("/tmp/dummy.root","RECREATE");
 807
 808   TFile *fIn = TFile::Open(currFile);
 809   if(!fIn){
 810     // not a severe condition but inciate that we have no information
 811     return kFALSE;
 812   }
 813   // find the tlists we want to be independtent of the name so use the Tkey
 814   TList* keyList = fIn->GetListOfKeys();
 815   Float_t memorySize = 0;
 816   Float_t diskSize = 0;
 817
 818   for(int i = 0;i < keyList->GetEntries();i++){
 819     TKey* ikey = (TKey*)keyList->At(i);
 820
 821     //    TList *list = dynamic_cast<TList*>(key->ReadObj());
 822     //    TNamed *name = dynamic_cast<TNamed*>(ikey->ReadObj());
 823     TDirectory *dir =  dynamic_cast<TDirectory*>(ikey->ReadObj());
 824
 825
 826
 827
 828     if(dir){
 829       Printf("%03d    : %60s %8d %8d ",i,dir->GetName(),ikey->GetObjlen(),ikey->GetNbytes());
 830       TList * dirKeyList = dir->GetListOfKeys();
 831       for(int j = 0;j<dirKeyList->GetEntries();j++){
 832         TKey* jkey = (TKey*)dirKeyList->At(j);
 833         TList *list =  dynamic_cast<TList*>(jkey->ReadObj());
 834
 835         memorySize += (Float_t)jkey->GetObjlen()/1024./1024.;
 836         diskSize +=  (Float_t)jkey->GetNbytes()/1024./1024.;
 837         if(list){
 838           Printf("%03d/%03d: %60s %5.2f MB %5.2f MB",i,j,list->GetName(),(Float_t)jkey->GetObjlen()/1024./1024.,(Float_t)jkey->GetNbytes()/1024./1024.);
 839           if(iDetail==i){
 840             for(int il = 0;il<list->GetEntries();il++){
 841               TObject *ob = list->At(il);
 842               if(fDummy){
 843                 fDummy->cd();
 844                 Int_t nBytesWrite = ob->Write();
 845                 Printf("%03d/%03d/%03d: %60s  %5.2f kB",i,j,il,ob->GetName(),(Float_t)nBytesWrite/1024.);
 846               }
 847             }
 848           }
 849         }
 850         else{
 851           Printf("%03d/%03d: %60s %5.2f MB %5.2f MB",i,j,jkey->GetName(),(Float_t)jkey->GetObjlen()/1024./1024.,(Float_t)jkey->GetNbytes()/1024./1024.);
 852         }
 853       }
 854     }
 855   }
 856   Printf("Total %5.2f MB %5.2f MB",memorySize,diskSize);
 857   return kTRUE;
 858 }
 859
 860
 861 Bool_t  AliAnalysisHelperJetTasks::Selected(Bool_t bSet,Bool_t bNew){
 862   //
 863   // Static helper task, (re)set event by event
 864   //
 865
 866
 867   static Bool_t bSelected = kTRUE; // if service task is not run we acccpet all
 868   if(bSet){
 869     bSelected = bNew;
 870   }
 871   return bSelected;
 872 }
 873
 874 Bool_t  AliAnalysisHelperJetTasks::IsCosmic(){
 875   return ((SelectInfo()&kIsCosmic)==kIsCosmic);
 876 }
 877
 878 Bool_t  AliAnalysisHelperJetTasks::IsPileUp(){
 879   return ((SelectInfo()&kIsPileUp)==kIsPileUp);
 880 }
 881
 882
 883 Bool_t  AliAnalysisHelperJetTasks::TestSelectInfo(UInt_t iMask){
 884   return ((SelectInfo()&iMask)==iMask);
 885 }
 886
 887
 888 Bool_t  AliAnalysisHelperJetTasks::TestEventClass(Int_t iMask){
 889   if(iMask==0)return kTRUE;
 890   return (EventClass()==iMask);
 891 }
 892
 893
 894 UInt_t  AliAnalysisHelperJetTasks::SelectInfo(Bool_t bSet,UInt_t iNew){
 895   //
 896   // set event by event the slection info
 897   //
 898
 899   static UInt_t iSelectInfo = 0; //
 900   if(bSet){
 901     iSelectInfo = iNew;
 902   }
 903   return iSelectInfo;
 904 }
 905
 906
 907 Int_t  AliAnalysisHelperJetTasks::EventClass(Bool_t bSet,Int_t iNew){
 908   //
 909   // gloab storage/access to event class
 910   //
 911
 912   static Int_t iEventClass = 0; //
 913   if(bSet){
 914     iEventClass = iNew;
 915   }
 916   return iEventClass;
 917 }
 918
 919
 920
 921
 922 //___________________________________________________________________________________________________________
 923
 924 Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01,const  TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
 925 {
 926   // ***
 927   // Event shape calculation
 928   // sona.pochybova@cern.ch
 929
 930   const Int_t kTracks = 1000;
 931   if(nTracks>kTracks)return kFALSE;
 932
 933   //variables for thrust calculation
 934   TVector3 pTrackPerp[kTracks];
 935   Double_t psum2 = 0;
 936
 937   TVector3 psum;
 938   TVector3 psum02;
 939   TVector3 psum03;
 940
 941   Double_t psum1 = 0;
 942   Double_t psum102 = 0;
 943   Double_t psum103 = 0;
 944
 945   Double_t thrust[kTracks] = {0.0};
 946   Double_t th = -3;
 947   Double_t thrust02[kTracks] = {0.0};
 948   Double_t th02 = -4;
 949   Double_t thrust03[kTracks] = {0.0};
 950   Double_t th03 = -5;
 951
 952   //Sphericity calculation variables
 953   TMatrixDSym m(3);
 954   Double_t s00 = 0;
 955   Double_t s01 = 0;
 956   Double_t s02 = 0;
 957
 958   Double_t s10 = 0;
 959   Double_t s11 = 0;
 960   Double_t s12 = 0;
 961
 962   Double_t s20 = 0;
 963   Double_t s21 = 0;
 964   Double_t s22 = 0;
 965
 966   Double_t ptot = 0;
 967
 968   Double_t c = -10;
 969
 970 //
 971 //loop for thrust calculation
 972 //
 973
 974   for(Int_t i = 0; i < nTracks; i++)
 975     {
 976       pTrackPerp[i].SetXYZ(pTrack[i].X(), pTrack[i].Y(), 0);
 977       psum2 += pTrackPerp[i].Mag();
 978     }
 979
 980   //additional starting axis
 981   TVector3 n02;
 982   n02 = pTrack[1].Unit();
 983   n02.SetZ(0.);
 984   TVector3 n03;
 985   n03 = pTrack[2].Unit();
 986   n03.SetZ(0.);
 987
 988   //switches for calculating thrust for different starting points
 989   Int_t switch1 = 1;
 990   Int_t switch2 = 1;
 991   Int_t switch3 = 1;
 992
 993   //indexes for iteration of different starting points
 994   Int_t l1 = 0;
 995   Int_t l2 = 0;
 996   Int_t l3 = 0;
 997
 998   //maximal number of iterations
 999   //  Int_t nMaxIter = 100;
1000
1001   for(Int_t k = 0; k < nTracks; k++)
1002     {
1003
1004       if(switch1 == 1){
1005         psum.SetXYZ(0., 0., 0.);
1006         psum1 = 0;
1007         for(Int_t i = 0; i < nTracks; i++)
1008           {
1009             psum1 += (TMath::Abs(n01.Dot(pTrackPerp[i])));
1010             if (n01.Dot(pTrackPerp[i]) > 0) psum += pTrackPerp[i];
1011             if (n01.Dot(pTrackPerp[i]) < 0) psum -= pTrackPerp[i];
1012           }
1013         thrust[l1] = psum1/psum2;
1014       }
1015
1016       if(switch2 == 1){
1017         psum02.SetXYZ(0., 0., 0.);
1018         psum102 = 0;
1019         for(Int_t i = 0; i < nTracks; i++)
1020           {
1021             psum102 += (TMath::Abs(n02.Dot(pTrackPerp[i])));
1022             if (n02.Dot(pTrackPerp[i]) > 0) psum02 += pTrackPerp[i];
1023             if (n02.Dot(pTrackPerp[i]) < 0) psum02 -= pTrackPerp[i];
1024           }
1025         thrust02[l2] = psum102/psum2;
1026       }
1027
1028       if(switch3 == 1){
1029         psum03.SetXYZ(0., 0., 0.);
1030         psum103 = 0;
1031         for(Int_t i = 0; i < nTracks; i++)
1032           {
1033             psum103 += (TMath::Abs(n03.Dot(pTrackPerp[i])));
1034             if (n03.Dot(pTrackPerp[i]) > 0) psum03 += pTrackPerp[i];
1035             if (n03.Dot(pTrackPerp[i]) < 0) psum03 -= pTrackPerp[i];
1036           }
1037         thrust03[l3] = psum103/psum2;
1038       }
1039
1040       //check whether thrust value converged
1041       if(TMath::Abs(th-thrust[l1]) < 10e-7){
1042         switch1 = 0;
1043       }
1044
1045       if(TMath::Abs(th02-thrust02[l2]) < 10e-7){
1046         switch2 = 0;
1047       }
1048
1049       if(TMath::Abs(th03-thrust03[l3]) < 10e-7){
1050         switch3 = 0;
1051       }
1052
1053       //if it didn't, continue with the calculation
1054       if(switch1 == 1){
1055         th = thrust[l1];
1056         n01 = psum.Unit();
1057         l1++;
1058       }
1059
1060       if(switch2 == 1){
1061         th02 = thrust02[l2];
1062         n02 = psum02.Unit();
1063         l2++;
1064       }
1065
1066       if(switch3 == 1){
1067         th03 = thrust03[l3];
1068         n03 = psum03.Unit();
1069         l3++;
1070       }
1071
1072       //if thrust values for all starting direction converged check if to the same value
1073       if(switch2 == 0 && switch1 == 0 && switch3 == 0){
1074         if(TMath::Abs(th-th02) < 10e-7 && TMath::Abs(th-th03) < 10e-7 && TMath::Abs(th02-th03) < 10e-7){
1075           eventShapes[0] = th;
1076           AliInfoGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),Form("===== THRUST VALUE FOUND AT %d :: %f\n", k, th));
1077           break;
1078         }
1079         //if they did not, reset switches
1080         else{
1081           switch1 = 1;
1082           //      th = -1.;
1083           switch2 = 1;
1084           //      th02 = -2.;
1085           switch3 = 1;
1086           //      th03 = -4.;
1087         }
1088       }
1089
1090       //      Printf("========== %d +++ th :: %f=============\n", l1, th);
1091       //      Printf("========== %d +++ th2 :: %f=============\n", l2, th02);
1092       //      Printf("========== %d +++ th3 :: %f=============\n", l3, th03);
1093
1094     }
1095
1096   //if no common limitng value was found, take the maximum and take the corresponding thrust axis
1097   if(switch1 == 1 && switch2 == 1 && switch3 == 1){
1098     eventShapes[0] = TMath::Max(thrust[l1-1], thrust02[l2-1]);
1099     eventShapes[0] = TMath::Max(eventShapes[0], thrust03[l3-1]);
1100     if(TMath::Abs(eventShapes[0]-thrust[l1-1]) < 10e-7)
1101       n01 = n01;
1102     if(TMath::Abs(eventShapes[0]-thrust02[l2-1]) < 10e-7)
1103       n01 = n02;
1104     if(TMath::Abs(eventShapes[0]-thrust03[l3-1]) < 10e-7)
1105       n01 = n03;
1106     Printf("NO LIMITING VALUE FOUND :: MAXIMUM = %f\n", eventShapes[0]);
1107   }
1108
1109 //
1110 //other event shapes variables
1111 //
1112   for(Int_t j = 0; j < nTracks; j++)
1113     {
1114       s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
1115       s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
1116       s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
1117
1118       s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
1119       s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
1120       s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
1121
1122       s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
1123       s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
1124       s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
1125
1126       ptot += pTrack[j].Mag();
1127     }
1128
1129   if(ptot > 0.)
1130     {
1131       m(0,0) = s00/ptot;
1132       m(0,1) = s01/ptot;
1133       m(0,2) = s02/ptot;
1134
1135       m(1,0) = s10/ptot;
1136       m(1,1) = s11/ptot;
1137       m(1,2) = s12/ptot;
1138
1139       m(2,0) = s20/ptot;
1140       m(2,1) = s21/ptot;
1141       m(2,2) = s22/ptot;
1142
1143       TMatrixDSymEigen eigen(m);
1144       TVectorD eigenVal = eigen.GetEigenValues();
1145
1146       Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
1147       eventShapes[1] = sphericity;
1148
1149       Double_t aplanarity = (3/2)*(eigenVal(2));
1150       eventShapes[2] = aplanarity;
1151
1152       c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
1153       eventShapes[3] = c;
1154     }
1155   return kTRUE;
1156 }
1157
1158
1159
1160  //__________________________________________________________________________________________________________________________
1161 // Trigger Decisions copid from PWG0/AliTriggerAnalysis
1162
1163
1164 Bool_t AliAnalysisHelperJetTasks::IsTriggerFired(const AliVEvent* aEv, Trigger trigger)
1165 {
1166   // checks if an event has been triggered
1167   // no usage of ofline trigger here yet
1168
1169   // here we do a dirty hack to take also into account the
1170   // missing trigger bits and Bunch crossing paatern for real data
1171
1172
1173   if(aEv->InheritsFrom("AliESDEvent")){
1174     const AliESDEvent *esd = (AliESDEvent*)aEv;
1175     switch (trigger)
1176       {
1177       case kAcceptAll:
1178         {
1179           return kTRUE;
1180           break;
1181         }
1182       case kMB1:
1183         {
1184           if(esd->GetFiredTriggerClasses().Contains("CINT1B-"))return kTRUE;
1185           // does the same but without or'ed V0s
1186           if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1187           if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE;
1188           // this is for simulated data
1189           if(esd->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;
1190           break;
1191         }
1192       case kMB2:
1193         {
1194           if(esd->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;
1195           break;
1196         }
1197       case kMB3:
1198         {
1199           if(esd->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;
1200           break;
1201         }
1202       case kSPDGFO:
1203         {
1204           if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1205           if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE;
1206           if(esd->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;
1207           break;
1208         }
1209       default:
1210         {
1211           Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1212           break;
1213         }
1214       }
1215   }
1216   else if(aEv->InheritsFrom("AliAODEvent")){
1217     const AliAODEvent *aod = (AliAODEvent*)aEv;
1218     switch (trigger)
1219       {
1220       case kAcceptAll:
1221         {
1222           return kTRUE;
1223           break;
1224         }
1225       case kMB1:
1226         {
1227           if(aod->GetFiredTriggerClasses().Contains("CINT1B"))return kTRUE;
1228           // does the same but without or'ed V0s
1229           if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1230           // for sim data
1231           if(aod->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;
1232           break;
1233         }
1234       case kMB2:
1235         {
1236           if(aod->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;
1237           break;
1238         }
1239       case kMB3:
1240         {
1241           if(aod->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;
1242           break;
1243         }
1244       case kSPDGFO:
1245         {
1246           if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1247           if(aod->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;
1248           break;
1249         }
1250       default:
1251         {
1252           Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1253           break;
1254         }
1255       }
1256   }
1257     return kFALSE;
1258 }
1259
1260
1261  AliAnalysisHelperJetTasks::MCProcessType  AliAnalysisHelperJetTasks::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
1262    //
1263    // fetch the process type from the mc header
1264    //
1265
1266
1267   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
1268
1269   if (!pythiaGenHeader) {
1270     //    printf(" AliAnalysisHelperJetTasks::GetProcessType : Unknown gen Header type). \n");
1271     return kInvalidProcess;
1272   }
1273
1274
1275   Int_t pythiaType = pythiaGenHeader->ProcessType();
1276   fgLastProcessType = pythiaType;
1277   MCProcessType globalType = kInvalidProcess;
1278
1279
1280   if (adebug) {
1281     printf(" AliAnalysisHelperJetTasks::GetProcessType : Pythia process type found: %d \n",pythiaType);
1282   }
1283
1284
1285   if(pythiaType==92||pythiaType==93){
1286     globalType = kSD;
1287   }
1288   else if(pythiaType==94){
1289     globalType = kDD;
1290   }
1291   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
1292   else {
1293     globalType = kND;
1294   }
1295   return globalType;
1296 }
1297
1298
1299  AliAnalysisHelperJetTasks::MCProcessType  AliAnalysisHelperJetTasks::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
1300   //
1301   // get the process type of the event.
1302   //
1303
1304   // can only read pythia headers, either directly or from cocktalil header
1305   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
1306
1307   if (!dpmJetGenHeader) {
1308     printf(" AliAnalysisHelperJetTasks::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
1309     return kInvalidProcess;
1310   }
1311
1312   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
1313   fgLastProcessType = dpmJetType;
1314   MCProcessType globalType = kInvalidProcess;
1315
1316
1317   if (adebug) {
1318     printf(" AliAnalysisHelperJetTasks::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
1319   }
1320
1321
1322   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
1323     globalType = kND;
1324   }
1325   else if (dpmJetType==5 || dpmJetType==6) {
1326     globalType = kSD;
1327   }
1328   else if (dpmJetType==7) {
1329     globalType = kDD;
1330   }
1331   return globalType;
1332 }
1333
1334
1335 Int_t AliAnalysisHelperJetTasks::GetPhiBin(Double_t phi,Int_t fNRPBins)
1336 {
1337     Int_t phibin=-1;
1338     if(!(TMath::Abs(phi)<=2*TMath::Pi()))return -1;
1339     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1340     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1341     return phibin;
1342 }
1343
1344 Double_t  AliAnalysisHelperJetTasks::ReactionPlane(Bool_t bSet,Double_t fNew){
1345   //
1346   // Static helper task, (re)set event by event
1347   //
1348
1349
1350   static Double_t fRP = 0; // if service task is not run we acccpet all
1351   if(bSet){
1352     fRP = fNew;
1353   }
1354   return fRP;
1355 }