AliPhysics  97344c9 (97344c9)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisHelperJetTasks.cxx
Go to the documentation of this file.
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>
51 #include "TMatrixDSym.h"
52 #include "TMatrixDSymEigen.h"
53 #include "TVector.h"
54 
55 using std::ifstream;
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) {
407  if(iDebug>10) Printf("No jets nRecJets %d nGenJets %d\n",nRecJets,nGenJets);
408  return;
409  }
410  AliAODJet *genJet = 0x0;
411  AliAODJet *recJet = 0x0;
412 
413  // loop over generated/embedded jets
414  for(Int_t ig=0; ig<nGenJets; ++ig){
415 
416  for(Int_t i=0; i<kClosestJetsN; ++i){
417  closestJets[i][0] = -1; // index
418  closestJets[i][1] = 1e6; // delta R
419  }
420 
421  genJet = (AliAODJet*)genJetsList->At(ig);
422  //if(!genJet || !JetSelected(genJet)) continue;
423  if(!genJet) {
424  if(iDebug>10) Printf("genJet %d doesnot exist",ig);
425  continue;
426  }
427 
428  // find N closest reconstructed jets
429  Double_t deltaR = 0.;
430  for(Int_t ir=0; ir<nRecJets; ++ir){
431  recJet = (AliAODJet*)recJetsList->At(ir);
432  //if(!recJet || !JetSelected(recJet)) continue;
433  if(!recJet) {
434  if(iDebug>10) Printf("recJet %d doesnot exist",ir);
435  continue;
436  }
437  deltaR = genJet->DeltaR(recJet);
438 
439  Int_t i=kClosestJetsN-1;
440  if(deltaR<closestJets[i][1] && deltaR<maxDist){
441  closestJets[i][0] = (Double_t) ir; // index
442  closestJets[i][1] = deltaR;
443 
444  // sort array (closest at the beginning)
445  while(i>=1 && closestJets[i][1]<closestJets[i-1][1]){
446  Double_t tmpArr[2];
447  for(Int_t j=0; j<2; j++){
448  tmpArr[j] = closestJets[i-1][j];
449  closestJets[i-1][j] = closestJets[i][j];
450  closestJets[i][j] = tmpArr[j];
451  }
452  i--;
453  }
454  }
455  } // end: loop over reconstructed jets
456 
457  // calculate fraction for the N closest jets
458  Double_t maxFraction = -1.; // maximum found fraction in one jets
459  Double_t cumFraction = 0.; // cummulated fraction of closest jets (for break condition)
460  Double_t fraction = 0.;
461  Int_t ir = -1; // index of close reconstruced jet
462 
463  for(Int_t irc=0; irc<kClosestJetsN; irc++){
464  ir = (Int_t)(closestJets[irc][0]);
465  if(ir<0 || ir>nRecJets-1) continue;
466  recJet = (AliAODJet*)recJetsList->At(ir);
467  if(!(recJet)) continue;
468 
469  fraction = GetFractionOfJet(recJet,genJet,mode);
470 
471  cumFraction += fraction;
472 
473  // check if jet fullfills current matching condition
474  if(fraction>maxFraction){
475  // avoid multiple links
476  for(Int_t ij=0; ij<ig; ++ij){
477  if(iMatchIndex[ij]==ir) continue;
478  }
479  // set index
480  maxFraction = fraction;
481  fPtFraction[ig] = fraction;
482  iMatchIndex[ig] = ir;
483  }
484  // break condition: less energy left as already found in one jet or
485  // as required for positiv matching
486  if(1-cumFraction<maxFraction) break;
487  } // end: loop over closest jets
488 
489  if(iMatchIndex[ig]<0){
490  if(iDebug) Printf("Matching failed for (gen) jet #%d", ig);
491  }
492  }
493 }
494 
495 Double_t AliAnalysisHelperJetTasks::GetFractionOfJet(const AliAODJet *recJet, const AliAODJet *genJet, Int_t mode){
496  //
497  // get the fraction of hte signal jet in the full jt
498  //
499  Double_t ptGen = genJet->Pt();
500  if(ptGen==0.) return 999.;
501 
502  Double_t ptAssocTracks = 0.; // sum of pT of tracks found in both jets
503 
504  // look at tracks inside jet
505  TRefArray *genTrackList = genJet->GetRefTracks();
506  TRefArray *recTrackList = recJet->GetRefTracks();
507  Int_t nTracksGenJet = genTrackList->GetEntriesFast();
508  Int_t nTracksRecJet = recTrackList->GetEntriesFast();
509 
510  // AliAODTrack* recTrack;
511  // AliAODTrack* genTrack;
512  AliVParticle* recTrack;
513  AliVParticle* genTrack;
514  for(Int_t ir=0; ir<nTracksRecJet; ++ir){
515  // recTrack = (AliAODTrack*)(recTrackList->At(ir));
516  recTrack = dynamic_cast<AliVParticle*>(recTrackList->At(ir));
517  if(!recTrack) continue;
518 
519  for(Int_t ig=0; ig<nTracksGenJet; ++ig){
520  // genTrack = (AliAODTrack*)(genTrackList->At(ig));
521  genTrack = dynamic_cast<AliVParticle*>(genTrackList->At(ig));
522  if(!genTrack) continue;
523 
524  // look if it points to the same track
525  if( (mode&1)!=0 && genTrack==recTrack){
526  ptAssocTracks += genTrack->Pt();
527  break;
528  }
529 
530  if( (mode&2)!=0
531  && genTrack->GetLabel()>-1
532  && recTrack->GetLabel()>-1
533  && genTrack->GetLabel()==recTrack->GetLabel()){
534 
535  ptAssocTracks += genTrack->Pt();
536  break;
537  }
538  }
539  }
540 
541  // calculate fraction
542  Double_t fraction = ptAssocTracks/ptGen;
543 
544  return fraction;
545 }
546 
547 
548 void AliAnalysisHelperJetTasks::MergeOutputDirs(const char* cFiles,const char* cPattern,const char *cOutFile,Bool_t bUpdate){
549  // Routine to merge only directories containing the pattern
550  //
551  TString outFile(cOutFile);
552  if(outFile.Length()==0)outFile = Form("%s.root",cPattern);
553  ifstream in1;
554  in1.open(cFiles);
555  // open all files
556  TList fileList;
557  char cFile[240];
558  while(in1>>cFile){// only open the first file
559  Printf("Adding %s",cFile);
560  TFile *f1 = TFile::Open(cFile);
561  fileList.Add(f1);
562  }
563 
564  TFile *fOut = 0;
565  if(fileList.GetEntries()){// open the first file
566  TFile* fIn = dynamic_cast<TFile*>(fileList.At(0));
567  if(!fIn){
568  Printf("Input File not Found");
569  return;
570  }
571  // fetch the keys for the directories
572  TList *ldKeys = fIn->GetListOfKeys();
573  for(int id = 0;id < ldKeys->GetEntries();id++){
574  // check if it is a directory
575  TKey *dKey = (TKey*)ldKeys->At(id);
576  TDirectory *dir = dynamic_cast<TDirectory*>(dKey->ReadObj());
577  if(dir){
578  TString dName(dir->GetName());
579  if(dName.Contains(cPattern)){
580  // open new file if first match
581  if(!fOut){
582  if(bUpdate)fOut = new TFile(outFile.Data(),"UPDATE");
583  else fOut = new TFile(outFile.Data(),"RECREATE");
584  }
585  // merge all the stuff that is in this directory
586  TList *llKeys = dir->GetListOfKeys();
587  TList *tmplist;
588  TMethodCall callEnv;
589 
590  fOut->cd();
591  TDirectory *dOut = fOut->mkdir(dir->GetName());
592 
593  for(int il = 0;il < llKeys->GetEntries();il++){
594  TKey *lKey = (TKey*)llKeys->At(il);
595  TObject *object = dynamic_cast<TObject*>(lKey->ReadObj());
596  // from TSeqCollections::Merge
597  if(!object)continue;
598  // If current object is not mergeable just skip it
599  if (!object->IsA()) {
600  continue;
601  }
602  callEnv.InitWithPrototype(object->IsA(), "Merge", "TCollection*");
603  if (!callEnv.IsValid()) {
604  continue;
605  }
606  // Current object mergeable - get corresponding objects in input lists
607  tmplist = new TList();
608  for(int i = 1; i < fileList.GetEntries();i++){
609  TDirectory *fInTmp = dynamic_cast<TDirectory*>(fileList.At(i));
610  if(!fInTmp){
611  Printf("File %d not found",i);
612  continue;
613  }
614  TDirectory *dTmp = (TDirectory*)fInTmp->Get(dName.Data());
615  if(!dTmp){
616  Printf("Directory %s not found",dName.Data());
617  continue;
618  }
619  TObject* oTmp = (TObject*)dTmp->Get(object->GetName());
620  if(!oTmp){
621  Printf("Object %s not found",object->GetName());
622  continue;
623  }
624  tmplist->Add(oTmp);
625  }
626  callEnv.SetParam((Long_t) tmplist);
627  callEnv.Execute(object);
628  delete tmplist;tmplist = 0;
629  dOut->cd();
630  object->Write(object->GetName(),TObject::kSingleKey);
631  }
632  }
633  }
634  }
635  if(fOut){
636  fOut->Close();
637  }
638  }
639 }
640 
641 void AliAnalysisHelperJetTasks::MergeOutput(const char* cFiles,const char* cDir,const char *cList,const char *cOutFile,Bool_t bUpdate){
642 
643  // This is used to merge the analysis-output from different
644  // data samples/pt_hard bins
645  // in case the eventweigth was set to xsection/ntrials already, this
646  // is not needed. Both methods only work in case we do not mix different
647  // pt_hard bins, and do not have overlapping bins
648 
649  const Int_t nMaxBins = 12;
650  // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
651  // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
652  // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
653  // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
654 
655  Float_t xsection[nMaxBins];
656  Float_t nTrials[nMaxBins];
657  Float_t sf[nMaxBins];
658  TList *lIn[nMaxBins];
659  TDirectory *dIn[nMaxBins];
660  TFile *fIn[nMaxBins];
661 
662  ifstream in1;
663  in1.open(cFiles);
664 
665  char cFile[120];
666  Int_t ibTotal = 0;
667  while(in1>>cFile){
668  fIn[ibTotal] = TFile::Open(cFile);
669  if(strlen(cDir)==0){
670  dIn[ibTotal] = gDirectory;
671  }
672  else{
673  dIn[ibTotal] = (TDirectory*)fIn[ibTotal]->Get(cDir);
674  }
675  if(!dIn[ibTotal]){
676  Printf("%s:%d No directory %s found, exiting...",__FILE__,__LINE__,cDir);
677  fIn[ibTotal]->ls();
678  return;
679  }
680 
681  lIn[ibTotal] = (TList*)dIn[ibTotal]->Get(cList);
682  Printf("Merging file %s %s",cFile, cDir);
683  if(!lIn[ibTotal]){
684  Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList);
685  fIn[ibTotal]->ls();
686  return;
687  }
688  TH1* hTrials = (TH1F*)lIn[ibTotal]->FindObject("fh1Trials");
689  if(!hTrials){
690  Printf("%s:%d fh1PtHard_Trials not found in list, exiting...",__FILE__,__LINE__);
691  return;
692  }
693  TProfile* hXsec = (TProfile*)lIn[ibTotal]->FindObject("fh1Xsec");
694  if(!hXsec){
695  Printf("%s:%d fh1Xsec not found in list, exiting...",__FILE__,__LINE__);
696  return;
697  }
698  xsection[ibTotal] = hXsec->GetBinContent(1);
699  nTrials[ibTotal] = hTrials->Integral();
700  sf[ibTotal] = xsection[ibTotal]/ nTrials[ibTotal];
701  ibTotal++;
702  }
703 
704  if(ibTotal==0){
705  Printf("%s:%d No files found for mergin, exiting",__FILE__,__LINE__);
706  return;
707  }
708 
709  TFile *fOut = 0;
710  if(bUpdate)fOut = new TFile(cOutFile,"UPDATE");
711  else fOut = new TFile(cOutFile,"RECREATE");
712  TDirectory *dOut = fOut->mkdir(dIn[0]->GetName());
713  dOut->cd();
714  TList *lOut = new TList();
715  lOut->SetName(lIn[0]->GetName());
716 
717  // for the start scale all...
718  for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
719  TH1 *h1Add = 0;
720  THnSparse *hnAdd = 0;
721  for(int ib = 0;ib < ibTotal;++ib){
722  // dynamic cast does not work with cint
723  TObject *h = lIn[ib]->At(ie);
724  if(h->InheritsFrom("TH1")){
725  TH1 *h1 = (TH1*)h;
726  if(ib==0){
727  h1Add = (TH1*)h1->Clone(h1->GetName());
728  h1Add->Scale(sf[ib]);
729  }
730  else{
731  h1Add->Add(h1,sf[ib]);
732  }
733  }
734  else if(h->InheritsFrom("THnSparse")){
735  THnSparse *hn = (THnSparse*)h;
736  if(ib==0){
737  hnAdd = (THnSparse*)hn->Clone(hn->GetName());
738  hnAdd->Scale(sf[ib]);
739  }
740  else{
741  hnAdd->Add(hn,sf[ib]);
742  }
743  }
744 
745 
746  }// ib
747  if(h1Add)lOut->Add(h1Add);
748  else if(hnAdd)lOut->Add(hnAdd);
749  }
750  dOut->cd();
751  lOut->Write(lOut->GetName(),TObject::kSingleKey);
752  fOut->Close();
753 }
754 
756  //
757  // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
758  // This is to called in Notify and should provide the path to the AOD/ESD file
759 
760  TString file(currFile);
761  fXsec = 0;
762  fTrials = 1;
763 
764  if(file.Contains("root_archive.zip#")){
765  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
766  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
767  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
768  file.Replace(pos+1,pos2-pos1,"");
769  }
770  else {
771  // not an archive take the basename....
772  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
773  }
774  Printf("%s",file.Data());
775 
776 
777 
778 
779  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
780  if(!fxsec){
781  // next trial fetch the histgram file
782  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
783  if(!fxsec){
784  // not a severe condition but inciate that we have no information
785  return kFALSE;
786  }
787  else{
788  // find the tlist we want to be independtent of the name so use the Tkey
789  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
790  if(!key){
791  fxsec->Close();
792  return kFALSE;
793  }
794  TList *list = dynamic_cast<TList*>(key->ReadObj());
795  if(!list){
796  fxsec->Close();
797  return kFALSE;
798  }
799  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
800  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
801  fxsec->Close();
802  }
803  } // no tree pyxsec.root
804  else {
805  TTree *xtree = (TTree*)fxsec->Get("Xsection");
806  if(!xtree){
807  fxsec->Close();
808  return kFALSE;
809  }
810  UInt_t ntrials = 0;
811  Double_t xsection = 0;
812  xtree->SetBranchAddress("xsection",&xsection);
813  xtree->SetBranchAddress("ntrials",&ntrials);
814  xtree->GetEntry(0);
815  fTrials = ntrials;
816  fXsec = xsection;
817  fxsec->Close();
818  }
819  return kTRUE;
820 }
821 
823 
824  //
825  // Print the size on disk and in memory occuppied by a directory
826  //
827 
828  TFile *fDummy = 0;
829  if(iDetail>=0)fDummy = TFile::Open("/tmp/dummy.root","RECREATE");
830 
831  TFile *fIn = TFile::Open(currFile);
832  if(!fIn){
833  // not a severe condition but inciate that we have no information
834  return kFALSE;
835  }
836  // find the tlists we want to be independtent of the name so use the Tkey
837  TList* keyList = fIn->GetListOfKeys();
838  Float_t memorySize = 0;
839  Float_t diskSize = 0;
840 
841  for(int i = 0;i < keyList->GetEntries();i++){
842  TKey* ikey = (TKey*)keyList->At(i);
843 
844  // TList *list = dynamic_cast<TList*>(key->ReadObj());
845  // TNamed *name = dynamic_cast<TNamed*>(ikey->ReadObj());
846  TDirectory *dir = dynamic_cast<TDirectory*>(ikey->ReadObj());
847 
848 
849 
850 
851  if(dir){
852  Printf("%03d : %60s %8d %8d ",i,dir->GetName(),ikey->GetObjlen(),ikey->GetNbytes());
853  TList * dirKeyList = dir->GetListOfKeys();
854  for(int j = 0;j<dirKeyList->GetEntries();j++){
855  TKey* jkey = (TKey*)dirKeyList->At(j);
856  TList *list = dynamic_cast<TList*>(jkey->ReadObj());
857 
858  memorySize += (Float_t)jkey->GetObjlen()/1024./1024.;
859  diskSize += (Float_t)jkey->GetNbytes()/1024./1024.;
860  if(list){
861  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.);
862  if(iDetail==i){
863  for(int il = 0;il<list->GetEntries();il++){
864  TObject *ob = list->At(il);
865  if(fDummy){
866  fDummy->cd();
867  Int_t nBytesWrite = ob->Write();
868  Printf("%03d/%03d/%03d: %60s %5.2f kB",i,j,il,ob->GetName(),(Float_t)nBytesWrite/1024.);
869  }
870  }
871  }
872  }
873  else{
874  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.);
875  }
876  }
877  }
878  }
879  Printf("Total %5.2f MB %5.2f MB",memorySize,diskSize);
880  return kTRUE;
881 }
882 
883 
885  //
886  // Static helper task, (re)set event by event
887  //
888 
889 
890  static Bool_t bSelected = kTRUE; // if service task is not run we acccpet all
891  if(bSet){
892  bSelected = bNew;
893  }
894  return bSelected;
895 }
896 
898  return ((SelectInfo()&kIsCosmic)==kIsCosmic);
899 }
900 
902  return ((SelectInfo()&kIsPileUp)==kIsPileUp);
903 }
904 
905 
907  return ((SelectInfo()&iMask)==iMask);
908 }
909 
910 
912  if(iMask==0)return kTRUE;
913  return (EventClass()==iMask);
914 }
915 
916 
918  //
919  // set event by event the slection info
920  //
921 
922  static UInt_t iSelectInfo = 0; //
923  if(bSet){
924  iSelectInfo = iNew;
925  }
926  return iSelectInfo;
927 }
928 
929 
931  //
932  // gloab storage/access to event class
933  //
934 
935  static Int_t iEventClass = 0; //
936  if(bSet){
937  iEventClass = iNew;
938  }
939  return iEventClass;
940 }
941 
942 
943 
944 
945 //___________________________________________________________________________________________________________
946 
947 Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01,const TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
948 {
949  // ***
950  // Event shape calculation
951  // sona.pochybova@cern.ch
952 
953  const Int_t kTracks = 1000;
954  if(nTracks>kTracks)return kFALSE;
955 
956  //variables for thrust calculation
957  TVector3 pTrackPerp[kTracks];
958  Double_t psum2 = 0;
959 
960  TVector3 psum;
961  TVector3 psum02;
962  TVector3 psum03;
963 
964  Double_t psum1 = 0;
965  Double_t psum102 = 0;
966  Double_t psum103 = 0;
967 
968  Double_t thrust[kTracks] = {0.0};
969  Double_t th = -3;
970  Double_t thrust02[kTracks] = {0.0};
971  Double_t th02 = -4;
972  Double_t thrust03[kTracks] = {0.0};
973  Double_t th03 = -5;
974 
975  //Sphericity calculation variables
976  TMatrixDSym m(3);
977  Double_t s00 = 0;
978  Double_t s01 = 0;
979  Double_t s02 = 0;
980 
981  Double_t s10 = 0;
982  Double_t s11 = 0;
983  Double_t s12 = 0;
984 
985  Double_t s20 = 0;
986  Double_t s21 = 0;
987  Double_t s22 = 0;
988 
989  Double_t ptot = 0;
990 
991  Double_t c = -10;
992 
993 //
994 //loop for thrust calculation
995 //
996 
997  for(Int_t i = 0; i < nTracks; i++)
998  {
999  pTrackPerp[i].SetXYZ(pTrack[i].X(), pTrack[i].Y(), 0);
1000  psum2 += pTrackPerp[i].Mag();
1001  }
1002 
1003  //additional starting axis
1004  TVector3 n02;
1005  n02 = pTrack[1].Unit();
1006  n02.SetZ(0.);
1007  TVector3 n03;
1008  n03 = pTrack[2].Unit();
1009  n03.SetZ(0.);
1010 
1011  //switches for calculating thrust for different starting points
1012  Int_t switch1 = 1;
1013  Int_t switch2 = 1;
1014  Int_t switch3 = 1;
1015 
1016  //indexes for iteration of different starting points
1017  Int_t l1 = 0;
1018  Int_t l2 = 0;
1019  Int_t l3 = 0;
1020 
1021  //maximal number of iterations
1022  // Int_t nMaxIter = 100;
1023 
1024  for(Int_t k = 0; k < nTracks; k++)
1025  {
1026 
1027  if(switch1 == 1){
1028  psum.SetXYZ(0., 0., 0.);
1029  psum1 = 0;
1030  for(Int_t i = 0; i < nTracks; i++)
1031  {
1032  psum1 += (TMath::Abs(n01.Dot(pTrackPerp[i])));
1033  if (n01.Dot(pTrackPerp[i]) > 0) psum += pTrackPerp[i];
1034  if (n01.Dot(pTrackPerp[i]) < 0) psum -= pTrackPerp[i];
1035  }
1036  thrust[l1] = psum1/psum2;
1037  }
1038 
1039  if(switch2 == 1){
1040  psum02.SetXYZ(0., 0., 0.);
1041  psum102 = 0;
1042  for(Int_t i = 0; i < nTracks; i++)
1043  {
1044  psum102 += (TMath::Abs(n02.Dot(pTrackPerp[i])));
1045  if (n02.Dot(pTrackPerp[i]) > 0) psum02 += pTrackPerp[i];
1046  if (n02.Dot(pTrackPerp[i]) < 0) psum02 -= pTrackPerp[i];
1047  }
1048  thrust02[l2] = psum102/psum2;
1049  }
1050 
1051  if(switch3 == 1){
1052  psum03.SetXYZ(0., 0., 0.);
1053  psum103 = 0;
1054  for(Int_t i = 0; i < nTracks; i++)
1055  {
1056  psum103 += (TMath::Abs(n03.Dot(pTrackPerp[i])));
1057  if (n03.Dot(pTrackPerp[i]) > 0) psum03 += pTrackPerp[i];
1058  if (n03.Dot(pTrackPerp[i]) < 0) psum03 -= pTrackPerp[i];
1059  }
1060  thrust03[l3] = psum103/psum2;
1061  }
1062 
1063  //check whether thrust value converged
1064  if(TMath::Abs(th-thrust[l1]) < 10e-7){
1065  switch1 = 0;
1066  }
1067 
1068  if(TMath::Abs(th02-thrust02[l2]) < 10e-7){
1069  switch2 = 0;
1070  }
1071 
1072  if(TMath::Abs(th03-thrust03[l3]) < 10e-7){
1073  switch3 = 0;
1074  }
1075 
1076  //if it didn't, continue with the calculation
1077  if(switch1 == 1){
1078  th = thrust[l1];
1079  n01 = psum.Unit();
1080  l1++;
1081  }
1082 
1083  if(switch2 == 1){
1084  th02 = thrust02[l2];
1085  n02 = psum02.Unit();
1086  l2++;
1087  }
1088 
1089  if(switch3 == 1){
1090  th03 = thrust03[l3];
1091  n03 = psum03.Unit();
1092  l3++;
1093  }
1094 
1095  //if thrust values for all starting direction converged check if to the same value
1096  if(switch2 == 0 && switch1 == 0 && switch3 == 0){
1097  if(TMath::Abs(th-th02) < 10e-7 && TMath::Abs(th-th03) < 10e-7 && TMath::Abs(th02-th03) < 10e-7){
1098  eventShapes[0] = th;
1099  AliInfoGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),Form("===== THRUST VALUE FOUND AT %d :: %f\n", k, th));
1100  break;
1101  }
1102  //if they did not, reset switches
1103  else{
1104  switch1 = 1;
1105  // th = -1.;
1106  switch2 = 1;
1107  // th02 = -2.;
1108  switch3 = 1;
1109  // th03 = -4.;
1110  }
1111  }
1112 
1113  // Printf("========== %d +++ th :: %f=============\n", l1, th);
1114  // Printf("========== %d +++ th2 :: %f=============\n", l2, th02);
1115  // Printf("========== %d +++ th3 :: %f=============\n", l3, th03);
1116 
1117  }
1118 
1119  //if no common limitng value was found, take the maximum and take the corresponding thrust axis
1120  if(switch1 == 1 && switch2 == 1 && switch3 == 1){
1121  eventShapes[0] = TMath::Max(thrust[l1-1], thrust02[l2-1]);
1122  eventShapes[0] = TMath::Max(eventShapes[0], thrust03[l3-1]);
1123  if(TMath::Abs(eventShapes[0]-thrust[l1-1]) < 10e-7)
1124  n01 = n01;
1125  if(TMath::Abs(eventShapes[0]-thrust02[l2-1]) < 10e-7)
1126  n01 = n02;
1127  if(TMath::Abs(eventShapes[0]-thrust03[l3-1]) < 10e-7)
1128  n01 = n03;
1129  Printf("NO LIMITING VALUE FOUND :: MAXIMUM = %f\n", eventShapes[0]);
1130  }
1131 
1132 //
1133 //other event shapes variables
1134 //
1135  for(Int_t j = 0; j < nTracks; j++)
1136  {
1137  s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
1138  s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
1139  s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
1140 
1141  s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
1142  s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
1143  s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
1144 
1145  s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
1146  s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
1147  s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
1148 
1149  ptot += pTrack[j].Mag();
1150  }
1151 
1152  if(ptot > 0.)
1153  {
1154  m(0,0) = s00/ptot;
1155  m(0,1) = s01/ptot;
1156  m(0,2) = s02/ptot;
1157 
1158  m(1,0) = s10/ptot;
1159  m(1,1) = s11/ptot;
1160  m(1,2) = s12/ptot;
1161 
1162  m(2,0) = s20/ptot;
1163  m(2,1) = s21/ptot;
1164  m(2,2) = s22/ptot;
1165 
1166  TMatrixDSymEigen eigen(m);
1167  TVectorD eigenVal = eigen.GetEigenValues();
1168 
1169  Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
1170  eventShapes[1] = sphericity;
1171 
1172  Double_t aplanarity = (3/2)*(eigenVal(2));
1173  eventShapes[2] = aplanarity;
1174 
1175  c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
1176  eventShapes[3] = c;
1177  }
1178  return kTRUE;
1179 }
1180 
1181 
1182 
1183  //__________________________________________________________________________________________________________________________
1184 // Trigger Decisions copid from PWG0/AliTriggerAnalysis
1185 
1186 
1188 {
1189  // checks if an event has been triggered
1190  // no usage of ofline trigger here yet
1191 
1192  // here we do a dirty hack to take also into account the
1193  // missing trigger bits and Bunch crossing paatern for real data
1194 
1195 
1196  if(aEv->InheritsFrom("AliESDEvent")){
1197  const AliESDEvent *esd = (AliESDEvent*)aEv;
1198  switch (trigger)
1199  {
1200  case kAcceptAll:
1201  {
1202  return kTRUE;
1203  break;
1204  }
1205  case kMB1:
1206  {
1207  if(esd->GetFiredTriggerClasses().Contains("CINT1B-"))return kTRUE;
1208  // does the same but without or'ed V0s
1209  if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1210  if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE;
1211  // this is for simulated data
1212  if(esd->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;
1213  break;
1214  }
1215  case kMB2:
1216  {
1217  if(esd->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;
1218  break;
1219  }
1220  case kMB3:
1221  {
1222  if(esd->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;
1223  break;
1224  }
1225  case kSPDGFO:
1226  {
1227  if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1228  if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE;
1229  if(esd->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;
1230  break;
1231  }
1232  default:
1233  {
1234  Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1235  break;
1236  }
1237  }
1238  }
1239  else if(aEv->InheritsFrom("AliAODEvent")){
1240  const AliAODEvent *aod = (AliAODEvent*)aEv;
1241  switch (trigger)
1242  {
1243  case kAcceptAll:
1244  {
1245  return kTRUE;
1246  break;
1247  }
1248  case kMB1:
1249  {
1250  if(aod->GetFiredTriggerClasses().Contains("CINT1B"))return kTRUE;
1251  // does the same but without or'ed V0s
1252  if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1253  // for sim data
1254  if(aod->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;
1255  break;
1256  }
1257  case kMB2:
1258  {
1259  if(aod->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;
1260  break;
1261  }
1262  case kMB3:
1263  {
1264  if(aod->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;
1265  break;
1266  }
1267  case kSPDGFO:
1268  {
1269  if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1270  if(aod->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;
1271  break;
1272  }
1273  default:
1274  {
1275  Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1276  break;
1277  }
1278  }
1279  }
1280  return kFALSE;
1281 }
1282 
1283 
1285  //
1286  // fetch the process type from the mc header
1287  //
1288 
1289 
1290  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
1291 
1292  if (!pythiaGenHeader) {
1293  // printf(" AliAnalysisHelperJetTasks::GetProcessType : Unknown gen Header type). \n");
1294  return kInvalidProcess;
1295  }
1296 
1297 
1298  Int_t pythiaType = pythiaGenHeader->ProcessType();
1299  fgLastProcessType = pythiaType;
1300  MCProcessType globalType = kInvalidProcess;
1301 
1302 
1303  if (adebug) {
1304  printf(" AliAnalysisHelperJetTasks::GetProcessType : Pythia process type found: %d \n",pythiaType);
1305  }
1306 
1307 
1308  if(pythiaType==92||pythiaType==93){
1309  globalType = kSD;
1310  }
1311  else if(pythiaType==94){
1312  globalType = kDD;
1313  }
1314  //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
1315  else {
1316  globalType = kND;
1317  }
1318  return globalType;
1319 }
1320 
1321 
1323  //
1324  // get the process type of the event.
1325  //
1326 
1327  // can only read pythia headers, either directly or from cocktalil header
1328  AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
1329 
1330  if (!dpmJetGenHeader) {
1331  printf(" AliAnalysisHelperJetTasks::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
1332  return kInvalidProcess;
1333  }
1334 
1335  Int_t dpmJetType = dpmJetGenHeader->ProcessType();
1336  fgLastProcessType = dpmJetType;
1337  MCProcessType globalType = kInvalidProcess;
1338 
1339 
1340  if (adebug) {
1341  printf(" AliAnalysisHelperJetTasks::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
1342  }
1343 
1344 
1345  if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
1346  globalType = kND;
1347  }
1348  else if (dpmJetType==5 || dpmJetType==6) {
1349  globalType = kSD;
1350  }
1351  else if (dpmJetType==7) {
1352  globalType = kDD;
1353  }
1354  return globalType;
1355 }
1356 
1357 
1359 {
1360  //
1361  // Method to get phi bin o reaction plane
1362  //
1363  Int_t phibin=-1;
1364  if(!(TMath::Abs(phi)<=2*TMath::Pi()))return -1;
1365  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1366  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1367  return phibin;
1368 }
1369 
1371  //
1372  // Static helper task, (re)set event by event
1373  //
1374 
1375 
1376  static Double_t fRP = 0; // if service task is not run we acccpet all
1377  if(bSet){
1378  fRP = fNew;
1379  }
1380  return fRP;
1381 }
double Double_t
Definition: External.C:58
static void GetJetMatching(const TList *genJetsList, const Int_t &kGenJets, const TList *recJetsList, const Int_t &kRecJets, TArrayI &iMatchIndex, TArrayF &fPtFraction, Int_t iDebug=0, Float_t maxDist=0.3, Int_t mode=1)
static UInt_t SelectInfo(Bool_t bSet=kFALSE, UInt_t iNew=0)
TSystem * gSystem
TList * list
TCanvas * c
Definition: TestFitELoss.C:172
static Int_t GetPhiBin(Double_t phi, Int_t fNRPbins)
static Bool_t Selected(Bool_t bSet=kFALSE, Bool_t bNew=kTRUE)
static Bool_t IsTriggerFired(const AliVEvent *aEsd, Trigger trigger)
static Double_t ReactionPlane(Bool_t bSet=kFALSE, Double_t fNew=0)
static Double_t GetFractionOfJet(const AliAODJet *recJet, const AliAODJet *genJet, Int_t mode=1)
int Int_t
Definition: External.C:63
static Int_t EventClass(Bool_t bSet=kFALSE, Int_t iNew=0)
static Bool_t TestSelectInfo(UInt_t iMask)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
static MCProcessType GetDPMjetEventProcessType(AliGenEventHeader *aHeader, Bool_t adebug=kFALSE)
static void PrintStack(AliMCEvent *mcEvent, Int_t iFirst=0, Int_t iLast=0, Int_t iMaxPrint=10)
Int_t mode
Definition: anaM.C:40
Double_t Mag(const AliVParticle &trk)
TObject * FindObject(int bin, const char *nameH, const TList *lst, Bool_t normPerEvent=kTRUE)
static void MergeOutputDirs(const char *cFiles, const char *cPattern, const char *cOutFile, Bool_t bUpdate=false)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
static void MergeOutput(const char *cFiles, const char *cDir="", const char *cList="", const char *cOutFile="allpt.root", Bool_t bUpdate=false)
TFile * file
static Bool_t PrintDirectorySize(const char *currFile, Int_t iDetail=-1)
unsigned short UShort_t
Definition: External.C:28
static Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials)
bool Bool_t
Definition: External.C:53
static Bool_t GetEventShapes(TVector3 &n01, const TVector3 *pTrack, Int_t nTracks, Double_t *eventShapes)
static void GetClosestJets(const AliAODJet *genJets, const Int_t &kGenJets, const AliAODJet *recJets, const Int_t &kRecJets, Int_t *iGenIndex, Int_t *iRecIndex, Int_t iDebug=0, Float_t maxDist=0.3)
static MCProcessType GetPythiaEventProcessType(AliGenEventHeader *aHeader, Bool_t adebug=kFALSE)
Definition: External.C:196
static Bool_t TestEventClass(Int_t iClass)