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