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