4c112b152a999d06ff71966f82b6fbfbac5c1546
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisHelperJetTasks.cxx
1
2 #include "TROOT.h"
3 #include "TDirectory.h"
4 #include "TKey.h"
5 #include "TList.h"
6 #include "TSystem.h"
7 #include "TH1F.h"
8 #include "TProfile.h"
9 #include "THnSparse.h"
10 #include "TSeqCollection.h"
11 #include "TMethodCall.h"
12 #include "TFile.h"
13 #include "TString.h"
14 #include "TArrayI.h"
15 #include "TArrayS.h"
16
17 #include "AliMCEvent.h"
18 #include "AliLog.h"
19 #include "AliESDEvent.h"
20 #include "AliAODJet.h"
21 #include "AliAODEvent.h"
22 #include "AliStack.h"
23 #include "AliGenEventHeader.h"
24 #include "AliGenCocktailEventHeader.h"
25 #include <AliGenDPMjetEventHeader.h>
26 #include "AliGenPythiaEventHeader.h"
27 #include <fstream>
28 #include <iostream>
29 #include "AliAnalysisHelperJetTasks.h"
30 #include "TMatrixDSym.h"
31 #include "TMatrixDSymEigen.h"
32 #include "TVector.h"
33
34 ClassImp(AliAnalysisHelperJetTasks)
35
36 Int_t AliAnalysisHelperJetTasks::fgLastProcessType = -1;
37
38  
39 AliGenPythiaEventHeader*  AliAnalysisHelperJetTasks::GetPythiaEventHeader(AliMCEvent *mcEvent){
40   
41   if(!mcEvent)return 0;
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) {
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__));
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){
60       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
61       return 0;
62     }
63   }
64   return pythiaGenHeader;
65
66 }
67
68
69 void 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
97
98
99 void AliAnalysisHelperJetTasks::GetClosestJets(AliAODJet *genJets,const Int_t &kGenJets,
100                                                AliAODJet *recJets,const Int_t &kRecJets,
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
115
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
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
143   
144   const int kMode = 3;
145
146   const Int_t nGenJets = TMath::Min(kMaxJets,kGenJets);
147   const Int_t nRecJets = TMath::Min(kMaxJets,kRecJets);
148
149   if(nRecJets==0||nGenJets==0)return;
150
151   // UShort_t *iFlag = new UShort_t[nGenJets*nRecJets];
152   UShort_t iFlag[kMaxJets*kMaxJets] = {0}; // all values set to zero
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){
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       }    
170       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
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     }
177     if(iRecIndex[ig]>=0)iFlag[ig*nRecJets+iRecIndex[ig]]+=1;
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     }
191     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]*nRecJets+ir]+=2;
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
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
228 void 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]);
333
334       if(kMode==3){
335         // we have a uniqie correlation
336         if(iFlag[ig*nRecJets+ir]==3){
337           iGenIndex[ir] = ig;
338           iRecIndex[ig] = ir;
339         }
340       }
341       else{
342         // we just take the correlation from on side
343         if((iFlag[ig*nRecJets+ir]&2)==2){
344           iGenIndex[ir] = ig;
345         }
346         if((iFlag[ig*nRecJets+ir]&1)==1){
347           iRecIndex[ig] = ir;
348         }
349       }
350     }
351     if(iDebug>1)printf("\n");
352   }
353 }
354
355 void AliAnalysisHelperJetTasks::GetJetMatching(TList *genJetsList, const Int_t &kGenJets,
356                                                TList *recJetsList, const Int_t &kRecJets,
357                                                TArrayI &iMatchIndex, TArrayF &fPtFraction,
358                                                Int_t iDebug, Float_t maxDist){
359
360                                             
361     // Matching jets from two lists
362     // Start with highest energetic jet from first list (generated/embedded)
363     // Calculate distance (\Delta R) to all jets from second list (reconstructed)
364     // Select N closest jets = kClosestJetsN
365     // Check energy fraction from jets from first list in jets from second list
366     // Matched jets = jet with largest energy fraction
367     // Store index of matched jet in TArrayI iMatchIndex
368                                             
369     // reset index
370     iMatchIndex.Reset(-1);
371     fPtFraction.Reset(-1.);
372     
373     // N closest jets: store list with index and \Delta R
374     const Int_t kClosestJetsN = 4; 
375     Double_t closestJets[kClosestJetsN][2]; //[][0] = index, [][1] = \Delta R
376     
377     for(Int_t i=0; i<kClosestJetsN; ++i){
378         for(Int_t j=0; j<2; ++j){
379             closestJets[i][j] = 1e6;
380         }
381     }
382
383     
384     const Int_t nGenJets = TMath::Min(genJetsList->GetEntries(),kGenJets);
385     const Int_t nRecJets = TMath::Min(recJetsList->GetEntries(),kRecJets);
386     if(nRecJets==0||nGenJets==0) return;
387     
388     AliAODJet *genJet = 0x0;
389     AliAODJet *recJet = 0x0;
390     
391     // loop over generated/embedded jets
392     for(Int_t ig=0; ig<nGenJets; ++ig){
393         genJet = (AliAODJet*)genJetsList->At(ig);
394         //if(!genJet || !JetSelected(genJet)) continue;
395         if(!genJet) continue;
396         
397         // find N closest reconstructed jets
398         Double_t deltaR = 0.;
399         for(Int_t ir=0; ir<nRecJets; ++ir){
400             recJet = (AliAODJet*)recJetsList->At(ir);
401             //if(!recJet || !JetSelected(recJet)) continue;
402             if(!recJet) continue;
403             
404             deltaR = genJet->DeltaR(recJet);
405             
406             Int_t i=kClosestJetsN-1;
407             if(deltaR<closestJets[i][1] && deltaR<maxDist){
408                 closestJets[i][0] = (Double_t) ir; // index
409                 closestJets[i][1] = deltaR;
410                 
411                 // sort array (closest at the beginning)
412                 while(i>=1 && closestJets[i][1]<closestJets[i-1][1]){
413                     Double_t tmpArr[2];
414                     for(Int_t j=0; j<2; j++){
415                        tmpArr[j] = closestJets[i-1][j];
416                        closestJets[i-1][j]   = closestJets[i][j];
417                        closestJets[i][j] = tmpArr[j];
418                     }
419                     i--;
420                 }
421             } 
422         } // end: loop over reconstructed jets
423         
424         // calculate fraction for the N closest jets
425         Double_t maxFraction = 0.; // maximum found fraction in one jets
426         Double_t cumFraction = 0.; // cummulated fraction of closest jets (for break condition)
427         Double_t fraction = 0.;
428         Int_t ir = -1;  // index of close reconstruced jet
429         
430         for(Int_t irc=0; irc<kClosestJetsN; irc++){
431             ir = (Int_t)(closestJets[irc][0]);
432             recJet = (AliAODJet*)recJetsList->At(ir);
433             if(!(recJet)) continue;
434             
435             fraction = GetFractionOfJet(recJet,genJet);
436             
437             cumFraction += fraction;
438             
439             // check if jet fullfills current matching condition
440             if(fraction>maxFraction){
441                 // avoid multiple links
442                 for(Int_t ij=0; ij<ig; ++ij){
443                     if(iMatchIndex[ij]==ir) continue;
444                 }
445                 // set index
446                 maxFraction = fraction;
447                 fPtFraction[ig] = fraction;                
448                 iMatchIndex[ig] = ir;
449             }
450             // break condition: less energy left as already found in one jet or
451             // as required for positiv matching
452             if(1-cumFraction<maxFraction) break;
453         } // end: loop over closest jets
454         
455         if(iMatchIndex[ig]<0){
456             if(iDebug) Printf("Matching failed for (gen) jet #%d", ig);
457         }
458     }
459 }
460
461 Double_t AliAnalysisHelperJetTasks::GetFractionOfJet(AliAODJet *recJet, AliAODJet *genJet){
462
463     Double_t ptGen = genJet->Pt();
464     if(ptGen==0.) return 999.;
465     
466     Double_t ptAssocTracks = 0.; // sum of pT of tracks found in both jets
467     
468     // look at tracks inside jet
469     TRefArray *genTrackList = genJet->GetRefTracks();
470     TRefArray *recTrackList = recJet->GetRefTracks();
471     Int_t nTracksGenJet = genTrackList->GetEntriesFast();
472     Int_t nTracksRecJet = recTrackList->GetEntriesFast();
473     
474     AliAODTrack* recTrack;
475     AliAODTrack* genTrack;
476     for(Int_t ir=0; ir<nTracksRecJet; ++ir){
477         recTrack = (AliAODTrack*)(recTrackList->At(ir));
478         if(!recTrack) continue;
479         
480         for(Int_t ig=0; ig<nTracksGenJet; ++ig){
481             genTrack = (AliAODTrack*)(genTrackList->At(ig));
482             if(!genTrack) continue;
483             
484             // look if it points to the same track
485             if(genTrack==recTrack){
486                 ptAssocTracks += genTrack->Pt();
487                 break;
488             }
489         }
490     }
491     
492     // calculate fraction
493     Double_t fraction = ptAssocTracks/ptGen;
494     
495     return fraction;
496 }
497
498
499 void  AliAnalysisHelperJetTasks::MergeOutputDirs(const char* cFiles,const char* cPattern,const char *cOutFile,Bool_t bUpdate){
500   // Routine to merge only directories containing the pattern
501   //
502   TString outFile(cOutFile);
503   if(outFile.Length()==0)outFile = Form("%s.root",cPattern);
504   ifstream in1;
505   in1.open(cFiles);
506   // open all files
507   TList fileList;
508   char cFile[240];
509   while(in1>>cFile){// only open the first file
510     Printf("Adding %s",cFile);
511     TFile *f1 = TFile::Open(cFile);
512     fileList.Add(f1);
513   }
514
515   TFile *fOut = 0;
516   if(fileList.GetEntries()){// open the first file
517     TFile* fIn = dynamic_cast<TFile*>(fileList.At(0));
518     if(!fIn){
519       Printf("Input File not Found");
520       return;
521     }
522     // fetch the keys for the directories
523     TList *ldKeys = fIn->GetListOfKeys();
524     for(int id = 0;id < ldKeys->GetEntries();id++){
525       // check if it is a directory
526       TKey *dKey = (TKey*)ldKeys->At(id);
527       TDirectory *dir = dynamic_cast<TDirectory*>(dKey->ReadObj());
528       if(dir){
529         TString dName(dir->GetName());
530         if(dName.Contains(cPattern)){
531           // open new file if first match
532           if(!fOut){
533             if(bUpdate)fOut = new TFile(outFile.Data(),"UPDATE");
534             else fOut = new TFile(outFile.Data(),"RECREATE");
535           }
536           // merge all the stuff that is in this directory
537           TList *llKeys = dir->GetListOfKeys();
538           TList *tmplist;
539           TMethodCall callEnv;
540
541           fOut->cd();
542           TDirectory *dOut = fOut->mkdir(dir->GetName());
543
544           for(int il = 0;il < llKeys->GetEntries();il++){
545             TKey *lKey = (TKey*)llKeys->At(il);
546             TObject *object = dynamic_cast<TObject*>(lKey->ReadObj());
547             //  from TSeqCollections::Merge
548             if(!object)continue;
549             // If current object is not mergeable just skip it
550             if (!object->IsA()) {
551               continue;
552             }
553             callEnv.InitWithPrototype(object->IsA(), "Merge", "TCollection*");
554             if (!callEnv.IsValid()) {
555               continue;
556             }
557             // Current object mergeable - get corresponding objects in input lists
558             tmplist = new TList();
559             for(int i = 1; i < fileList.GetEntries();i++){
560               TDirectory *fInTmp = dynamic_cast<TDirectory*>(fileList.At(i)); 
561               if(!fInTmp){
562                 Printf("File %d not found",i);
563                 continue;
564               }
565               TDirectory *dTmp = (TDirectory*)fInTmp->Get(dName.Data());
566               if(!dTmp){
567                 Printf("Directory %s not found",dName.Data());
568                 continue;
569               }
570               TObject* oTmp = (TObject*)dTmp->Get(object->GetName());
571               if(!oTmp){
572                 Printf("Object %s not found",object->GetName());
573                 continue;
574               }
575               tmplist->Add(oTmp);
576             }
577             callEnv.SetParam((Long_t) tmplist);
578             callEnv.Execute(object);
579             delete tmplist;tmplist = 0;
580             dOut->cd();
581             object->Write(object->GetName(),TObject::kSingleKey);
582           }
583         }
584       }
585     }
586     if(fOut){
587       fOut->Close();
588     }
589   }
590 }
591
592 void  AliAnalysisHelperJetTasks::MergeOutput(char* cFiles, char* cDir, char *cList,char *cOutFile,Bool_t bUpdate){
593
594   // This is used to merge the analysis-output from different 
595   // data samples/pt_hard bins
596   // in case the eventweigth was set to xsection/ntrials already, this
597   // is not needed. Both methods only work in case we do not mix different 
598   // pt_hard bins, and do not have overlapping bins
599
600   const Int_t nMaxBins = 12;
601   // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
602   // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
603   // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
604   // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
605
606   Float_t xsection[nMaxBins];
607   Float_t nTrials[nMaxBins];
608   Float_t sf[nMaxBins];
609   TList *lIn[nMaxBins];
610   TDirectory *dIn[nMaxBins];
611   TFile *fIn[nMaxBins];
612
613   ifstream in1;
614   in1.open(cFiles);
615
616   char cFile[120];
617   Int_t ibTotal = 0;
618   while(in1>>cFile){
619     fIn[ibTotal] = TFile::Open(cFile);
620     if(strlen(cDir)==0){
621       dIn[ibTotal] = gDirectory;
622     }
623     else{
624       dIn[ibTotal] = (TDirectory*)fIn[ibTotal]->Get(cDir);
625     }
626     if(!dIn[ibTotal]){
627       Printf("%s:%d No directory %s found, exiting...",__FILE__,__LINE__,cDir);
628       fIn[ibTotal]->ls();
629       return;
630     }
631
632     lIn[ibTotal] = (TList*)dIn[ibTotal]->Get(cList);
633     Printf("Merging file %s %s",cFile, cDir);
634     if(!lIn[ibTotal]){
635       Printf("%s:%d No list %s found, exiting...",__FILE__,__LINE__,cList);
636       fIn[ibTotal]->ls();
637       return;
638     }
639     TH1* hTrials = (TH1F*)lIn[ibTotal]->FindObject("fh1Trials");
640     if(!hTrials){
641       Printf("%s:%d fh1PtHard_Trials not found in list, exiting...",__FILE__,__LINE__);
642       return;
643     }
644     TProfile* hXsec = (TProfile*)lIn[ibTotal]->FindObject("fh1Xsec");
645     if(!hXsec){
646       Printf("%s:%d fh1Xsec  not found in list, exiting...",__FILE__,__LINE__);
647       return;
648     }
649     xsection[ibTotal] = hXsec->GetBinContent(1);
650     nTrials[ibTotal] = hTrials->Integral();
651     sf[ibTotal] = xsection[ibTotal]/ nTrials[ibTotal];
652     ibTotal++;
653   }
654
655   if(ibTotal==0){
656     Printf("%s:%d No files found for mergin, exiting",__FILE__,__LINE__);
657     return;
658   }
659
660   TFile *fOut = 0;
661   if(bUpdate)fOut = new TFile(cOutFile,"UPDATE");
662   else fOut = new TFile(cOutFile,"RECREATE");
663   TDirectory *dOut = fOut->mkdir(dIn[0]->GetName());
664   dOut->cd();
665   TList *lOut = new TList();
666   lOut->SetName(lIn[0]->GetName());
667
668   // for the start scale all...
669   for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
670     TH1 *h1Add = 0;
671     THnSparse *hnAdd = 0;
672     for(int ib = 0;ib < ibTotal;++ib){
673       // dynamic cast does not work with cint
674       TObject *h = lIn[ib]->At(ie);
675       if(h->InheritsFrom("TH1")){
676         TH1 *h1 = (TH1*)h;
677         if(ib==0){
678           h1Add = (TH1*)h1->Clone(h1->GetName());
679           h1Add->Scale(sf[ib]);
680         }
681         else{
682           h1Add->Add(h1,sf[ib]);
683         }
684       }
685       else if(h->InheritsFrom("THnSparse")){
686         THnSparse *hn = (THnSparse*)h;
687         if(ib==0){
688           hnAdd = (THnSparse*)hn->Clone(hn->GetName());
689           hnAdd->Scale(sf[ib]);
690         }
691         else{
692           hnAdd->Add(hn,sf[ib]);
693         }
694       }
695       
696
697     }// ib
698     if(h1Add)lOut->Add(h1Add);
699     else if(hnAdd)lOut->Add(hnAdd);
700   }
701   dOut->cd();
702   lOut->Write(lOut->GetName(),TObject::kSingleKey);
703   fOut->Close();
704 }
705
706 Bool_t AliAnalysisHelperJetTasks::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
707   //
708   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
709   // This is to called in Notify and should provide the path to the AOD/ESD file
710
711   TString file(currFile);  
712   fXsec = 0;
713   fTrials = 1;
714
715   if(file.Contains("root_archive.zip#")){
716     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
717     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
718     file.Replace(pos+1,20,"");
719   }
720   else {
721     // not an archive take the basename....
722     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
723   }
724   Printf("%s",file.Data());
725   
726  
727    
728
729   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
730   if(!fxsec){
731     // next trial fetch the histgram file
732     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
733     if(!fxsec){
734         // not a severe condition but inciate that we have no information
735       return kFALSE;
736     }
737     else{
738       // find the tlist we want to be independtent of the name so use the Tkey
739       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
740       if(!key){
741         fxsec->Close();
742         return kFALSE;
743       }
744       TList *list = dynamic_cast<TList*>(key->ReadObj());
745       if(!list){
746         fxsec->Close();
747         return kFALSE;
748       }
749       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
750       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
751       fxsec->Close();
752     }
753   } // no tree pyxsec.root
754   else {
755     TTree *xtree = (TTree*)fxsec->Get("Xsection");
756     if(!xtree){
757       fxsec->Close();
758       return kFALSE;
759     }
760     UInt_t   ntrials  = 0;
761     Double_t  xsection  = 0;
762     xtree->SetBranchAddress("xsection",&xsection);
763     xtree->SetBranchAddress("ntrials",&ntrials);
764     xtree->GetEntry(0);
765     fTrials = ntrials;
766     fXsec = xsection;
767     fxsec->Close();
768   }
769   return kTRUE;
770 }
771
772 Bool_t AliAnalysisHelperJetTasks::PrintDirectorySize(const char* currFile){
773
774   TFile *fIn = TFile::Open(currFile);
775   if(!fIn){
776     // not a severe condition but inciate that we have no information
777     return kFALSE;
778   }
779   // find the tlists we want to be independtent of the name so use the Tkey
780   TList* keyList = fIn->GetListOfKeys();
781   Float_t memorySize = 0;
782   Float_t diskSize = 0;
783
784   for(int i = 0;i < keyList->GetEntries();i++){
785     TKey* ikey = (TKey*)keyList->At(i); 
786     
787     //    TList *list = dynamic_cast<TList*>(key->ReadObj());
788     //    TNamed *name = dynamic_cast<TNamed*>(ikey->ReadObj());
789     TDirectory *dir =  dynamic_cast<TDirectory*>(ikey->ReadObj());
790     
791
792
793
794     if(dir){
795       Printf("%03d    : %60s %8d %8d ",i,dir->GetName(),ikey->GetObjlen(),ikey->GetNbytes());
796       TList * dirKeyList = dir->GetListOfKeys();
797       for(int j = 0;j<dirKeyList->GetEntries();j++){
798         TKey* jkey = (TKey*)dirKeyList->At(j); 
799         TList *list =  dynamic_cast<TList*>(jkey->ReadObj());
800
801         memorySize += (Float_t)jkey->GetObjlen()/1024./1024.;
802         diskSize +=  (Float_t)jkey->GetNbytes()/1024./1024.;
803         if(list){
804           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.);
805         }
806         else{
807           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.);
808         }
809       }
810     }
811   }
812   Printf("Total %5.2f MB %5.2f MB",memorySize,diskSize);
813   return kTRUE;
814 }
815
816
817 Bool_t  AliAnalysisHelperJetTasks::Selected(Bool_t bSet,Bool_t bNew){
818   static Bool_t bSelected = kTRUE; // if service task is not run we acccpet all
819   if(bSet){
820     bSelected = bNew;
821   }
822   return bSelected;
823 }
824
825 Bool_t  AliAnalysisHelperJetTasks::IsCosmic(){
826   return ((SelectInfo()&kIsCosmic)==kIsCosmic);
827 }
828
829 Bool_t  AliAnalysisHelperJetTasks::IsPileUp(){
830   return ((SelectInfo()&kIsPileUp)==kIsPileUp);
831 }
832
833
834 Bool_t  AliAnalysisHelperJetTasks::TestSelectInfo(UInt_t iMask){
835   return ((SelectInfo()&iMask)==iMask);
836 }
837
838
839 Bool_t  AliAnalysisHelperJetTasks::TestEventClass(Int_t iMask){
840   if(iMask==0)return kTRUE;
841   return (EventClass()==iMask);
842 }
843
844
845 UInt_t  AliAnalysisHelperJetTasks::SelectInfo(Bool_t bSet,UInt_t iNew){
846   static UInt_t iSelectInfo = 0; //
847   if(bSet){
848     iSelectInfo = iNew;
849   }
850   return iSelectInfo;
851 }
852
853
854 Int_t  AliAnalysisHelperJetTasks::EventClass(Bool_t bSet,Int_t iNew){
855   static Int_t iEventClass = 0; //
856   if(bSet){
857     iEventClass = iNew;
858   }
859   return iEventClass;
860 }
861
862
863
864
865 //___________________________________________________________________________________________________________
866
867 Bool_t AliAnalysisHelperJetTasks::GetEventShapes(TVector3 &n01, TVector3 * pTrack, Int_t nTracks, Double_t * eventShapes)
868 {       
869   // ***
870   // Event shape calculation
871   // sona.pochybova@cern.ch
872
873   const Int_t kTracks = 1000;
874   if(nTracks>kTracks)return kFALSE;
875
876   //variables for thrust calculation
877   TVector3 pTrackPerp[kTracks];
878   Double_t psum2 = 0;
879
880   TVector3 psum;
881   TVector3 psum02;
882   TVector3 psum03;
883
884   Double_t psum1 = 0;
885   Double_t psum102 = 0;
886   Double_t psum103 = 0;
887
888   Double_t thrust[kTracks] = {0.0};
889   Double_t th = -3;
890   Double_t thrust02[kTracks] = {0.0};
891   Double_t th02 = -4;
892   Double_t thrust03[kTracks] = {0.0};
893   Double_t th03 = -5;
894
895   //Sphericity calculation variables
896   TMatrixDSym m(3);
897   Double_t s00 = 0;
898   Double_t s01 = 0;
899   Double_t s02 = 0;
900   
901   Double_t s10 = 0;
902   Double_t s11 = 0;
903   Double_t s12 = 0;
904   
905   Double_t s20 = 0;
906   Double_t s21 = 0;
907   Double_t s22 = 0;
908   
909   Double_t ptot = 0;
910   
911   Double_t c = -10;
912
913 //
914 //loop for thrust calculation  
915 //
916     
917   for(Int_t i = 0; i < nTracks; i++)
918     {
919       pTrackPerp[i].SetXYZ(pTrack[i].X(), pTrack[i].Y(), 0);
920       psum2 += pTrackPerp[i].Mag();
921     }
922
923   //additional starting axis    
924   TVector3 n02;
925   n02 = pTrack[1].Unit();
926   n02.SetZ(0.);   
927   TVector3 n03;
928   n03 = pTrack[2].Unit();
929   n03.SetZ(0.);
930
931   //switches for calculating thrust for different starting points
932   Int_t switch1 = 1;
933   Int_t switch2 = 1;
934   Int_t switch3 = 1;
935
936   //indexes for iteration of different starting points
937   Int_t l1 = 0;
938   Int_t l2 = 0;
939   Int_t l3 = 0;
940
941   //maximal number of iterations
942   //  Int_t nMaxIter = 100;
943  
944   for(Int_t k = 0; k < nTracks; k++)
945     {  
946       
947       if(switch1 == 1){
948         psum.SetXYZ(0., 0., 0.);
949         psum1 = 0;
950         for(Int_t i = 0; i < nTracks; i++)
951           {
952             psum1 += (TMath::Abs(n01.Dot(pTrackPerp[i])));
953             if (n01.Dot(pTrackPerp[i]) > 0) psum += pTrackPerp[i];
954             if (n01.Dot(pTrackPerp[i]) < 0) psum -= pTrackPerp[i];
955           }
956         thrust[l1] = psum1/psum2;
957       }
958
959       if(switch2 == 1){
960         psum02.SetXYZ(0., 0., 0.);
961         psum102 = 0;
962         for(Int_t i = 0; i < nTracks; i++)
963           {
964             psum102 += (TMath::Abs(n02.Dot(pTrackPerp[i])));
965             if (n02.Dot(pTrackPerp[i]) > 0) psum02 += pTrackPerp[i];
966             if (n02.Dot(pTrackPerp[i]) < 0) psum02 -= pTrackPerp[i];
967           }
968         thrust02[l2] = psum102/psum2;
969       }
970       
971       if(switch3 == 1){
972         psum03.SetXYZ(0., 0., 0.);
973         psum103 = 0;
974         for(Int_t i = 0; i < nTracks; i++)
975           {
976             psum103 += (TMath::Abs(n03.Dot(pTrackPerp[i])));
977             if (n03.Dot(pTrackPerp[i]) > 0) psum03 += pTrackPerp[i];
978             if (n03.Dot(pTrackPerp[i]) < 0) psum03 -= pTrackPerp[i];
979           }
980         thrust03[l3] = psum103/psum2;
981       }
982
983       //check whether thrust value converged    
984       if(TMath::Abs(th-thrust[l1]) < 10e-7){
985         switch1 = 0;
986       }
987       
988       if(TMath::Abs(th02-thrust02[l2]) < 10e-7){
989         switch2 = 0;
990       }
991
992       if(TMath::Abs(th03-thrust03[l3]) < 10e-7){
993         switch3 = 0;
994       }
995
996       //if it didn't, continue with the calculation
997       if(switch1 == 1){
998         th = thrust[l1]; 
999         n01 = psum.Unit();
1000         l1++;
1001       }
1002
1003       if(switch2 == 1){
1004         th02 = thrust02[l2];
1005         n02 = psum02.Unit();
1006         l2++;
1007       }
1008
1009       if(switch3 == 1){
1010         th03 = thrust03[l3];
1011         n03 = psum03.Unit();
1012         l3++;
1013       }
1014
1015       //if thrust values for all starting direction converged check if to the same value
1016       if(switch2 == 0 && switch1 == 0 && switch3 == 0){
1017         if(TMath::Abs(th-th02) < 10e-7 && TMath::Abs(th-th03) < 10e-7 && TMath::Abs(th02-th03) < 10e-7){
1018           eventShapes[0] = th;
1019           AliInfoGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),Form("===== THRUST VALUE FOUND AT %d :: %f\n", k, th));
1020           break;
1021         }
1022         //if they did not, reset switches
1023         else{
1024           switch1 = 1;
1025           //      th = -1.;
1026           switch2 = 1;
1027           //      th02 = -2.;
1028           switch3 = 1;
1029           //      th03 = -4.;
1030         }
1031       }
1032
1033       //      Printf("========== %d +++ th :: %f=============\n", l1, th);
1034       //      Printf("========== %d +++ th2 :: %f=============\n", l2, th02);
1035       //      Printf("========== %d +++ th3 :: %f=============\n", l3, th03);
1036       
1037     }
1038
1039   //if no common limitng value was found, take the maximum and take the corresponding thrust axis
1040   if(switch1 == 1 && switch2 == 1 && switch3 == 1){
1041     eventShapes[0] = TMath::Max(thrust[l1-1], thrust02[l2-1]);
1042     eventShapes[0] = TMath::Max(eventShapes[0], thrust03[l3-1]);
1043     if(TMath::Abs(eventShapes[0]-thrust[l1-1]) < 10e-7)
1044       n01 = n01;
1045     if(TMath::Abs(eventShapes[0]-thrust02[l2-1]) < 10e-7)
1046       n01 = n02;
1047     if(TMath::Abs(eventShapes[0]-thrust03[l3-1]) < 10e-7)
1048       n01 = n03;
1049     Printf("NO LIMITING VALUE FOUND :: MAXIMUM = %f\n", eventShapes[0]);
1050   }
1051
1052 //
1053 //other event shapes variables
1054 //
1055   for(Int_t j = 0; j < nTracks; j++)
1056     {  
1057       s00 = s00 + (pTrack[j].Px()*pTrack[j].Px())/pTrack[j].Mag();
1058       s01 = s01 + (pTrack[j].Px()*pTrack[j].Py())/pTrack[j].Mag();
1059       s02 = s02 + (pTrack[j].Px()*pTrack[j].Pz())/pTrack[j].Mag();
1060       
1061       s10 = s10 + (pTrack[j].Py()*pTrack[j].Px())/pTrack[j].Mag();
1062       s11 = s11 + (pTrack[j].Py()*pTrack[j].Py())/pTrack[j].Mag();
1063       s12 = s12 + (pTrack[j].Py()*pTrack[j].Pz())/pTrack[j].Mag();
1064       
1065       s20 = s20 + (pTrack[j].Pz()*pTrack[j].Px())/pTrack[j].Mag();
1066       s21 = s21 + (pTrack[j].Pz()*pTrack[j].Py())/pTrack[j].Mag();
1067       s22 = s22 + (pTrack[j].Pz()*pTrack[j].Pz())/pTrack[j].Mag();
1068       
1069       ptot += pTrack[j].Mag();
1070     }
1071
1072   if(ptot > 0.)
1073     {
1074       m(0,0) = s00/ptot;
1075       m(0,1) = s01/ptot;
1076       m(0,2) = s02/ptot;
1077
1078       m(1,0) = s10/ptot;
1079       m(1,1) = s11/ptot;
1080       m(1,2) = s12/ptot;
1081
1082       m(2,0) = s20/ptot;
1083       m(2,1) = s21/ptot;
1084       m(2,2) = s22/ptot;
1085
1086       TMatrixDSymEigen eigen(m);
1087       TVectorD eigenVal = eigen.GetEigenValues();
1088
1089       Double_t sphericity = (3/2)*(eigenVal(2)+eigenVal(1));
1090       eventShapes[1] = sphericity;
1091
1092       Double_t aplanarity = (3/2)*(eigenVal(2));
1093       eventShapes[2] = aplanarity;
1094
1095       c = 3*(eigenVal(0)*eigenVal(1)+eigenVal(0)*eigenVal(2)+eigenVal(1)*eigenVal(2));
1096       eventShapes[3] = c;
1097     }
1098   return kTRUE;
1099 }
1100   
1101
1102
1103  //__________________________________________________________________________________________________________________________
1104 // Trigger Decisions copid from PWG0/AliTriggerAnalysis
1105
1106
1107 Bool_t AliAnalysisHelperJetTasks::IsTriggerFired(const AliVEvent* aEv, Trigger trigger)
1108 {
1109   // checks if an event has been triggered
1110   // no usage of ofline trigger here yet
1111   
1112   // here we do a dirty hack to take also into account the
1113   // missing trigger bits and Bunch crossing paatern for real data 
1114
1115
1116   if(aEv->InheritsFrom("AliESDEvent")){
1117     const AliESDEvent *esd = (AliESDEvent*)aEv;
1118     switch (trigger)
1119       {
1120       case kAcceptAll:
1121         {
1122           return kTRUE;
1123           break;
1124         }
1125       case kMB1:
1126         {
1127           if(esd->GetFiredTriggerClasses().Contains("CINT1B-"))return kTRUE;
1128           // does the same but without or'ed V0s
1129           if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;  
1130           if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE; 
1131           // this is for simulated data
1132           if(esd->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;   
1133           break;
1134         }
1135       case kMB2:
1136         {
1137           if(esd->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;   
1138           break;
1139         }
1140       case kMB3:
1141         {
1142           if(esd->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;   
1143           break;
1144         }
1145       case kSPDGFO:
1146         {
1147           if(esd->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;
1148           if(esd->GetFiredTriggerClasses().Contains("CINT6B-"))return kTRUE; 
1149           if(esd->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;
1150           break;
1151         }
1152       default:
1153         {
1154           Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1155           break;
1156         }
1157       }
1158   }
1159   else if(aEv->InheritsFrom("AliAODEvent")){
1160     const AliAODEvent *aod = (AliAODEvent*)aEv;
1161     switch (trigger)
1162       {
1163       case kAcceptAll:
1164         {
1165           return kTRUE;
1166           break;
1167         }
1168       case kMB1:
1169         {
1170           if(aod->GetFiredTriggerClasses().Contains("CINT1B"))return kTRUE;
1171           // does the same but without or'ed V0s
1172           if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;   
1173           // for sim data
1174           if(aod->GetFiredTriggerClasses().Contains("MB1"))return kTRUE;   
1175           break;
1176         }
1177       case kMB2:
1178         {
1179           if(aod->GetFiredTriggerClasses().Contains("MB2"))return kTRUE;   
1180           break;
1181         }
1182       case kMB3:
1183         {
1184           if(aod->GetFiredTriggerClasses().Contains("MB3"))return kTRUE;   
1185           break;
1186         }
1187       case kSPDGFO:
1188         {
1189           if(aod->GetFiredTriggerClasses().Contains("CSMBB"))return kTRUE;        
1190           if(aod->GetFiredTriggerClasses().Contains("GFO"))return kTRUE;   
1191           break;
1192         }
1193       default:
1194         {
1195           Printf("IsEventTriggered: ERROR: Trigger type %d not implemented in this method", (Int_t) trigger);
1196           break;
1197         }
1198       }
1199   }
1200     return kFALSE;
1201 }
1202
1203
1204  AliAnalysisHelperJetTasks::MCProcessType  AliAnalysisHelperJetTasks::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
1205
1206   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader);
1207
1208   if (!pythiaGenHeader) {
1209     //    printf(" AliAnalysisHelperJetTasks::GetProcessType : Unknown gen Header type). \n");
1210     return kInvalidProcess;
1211   }
1212
1213
1214   Int_t pythiaType = pythiaGenHeader->ProcessType();
1215   fgLastProcessType = pythiaType;
1216   MCProcessType globalType = kInvalidProcess;  
1217
1218
1219   if (adebug) {
1220     printf(" AliAnalysisHelperJetTasks::GetProcessType : Pythia process type found: %d \n",pythiaType);
1221   }
1222
1223
1224   if(pythiaType==92||pythiaType==93){
1225     globalType = kSD;
1226   }
1227   else if(pythiaType==94){
1228     globalType = kDD;
1229   }
1230   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
1231   else {
1232     globalType = kND;
1233   }
1234   return globalType;
1235 }
1236
1237
1238  AliAnalysisHelperJetTasks::MCProcessType  AliAnalysisHelperJetTasks::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) {
1239   //
1240   // get the process type of the event.
1241   //
1242
1243   // can only read pythia headers, either directly or from cocktalil header
1244   AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader);
1245
1246   if (!dpmJetGenHeader) {
1247     printf(" AliAnalysisHelperJetTasks::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n");
1248     return kInvalidProcess;
1249   }
1250
1251   Int_t dpmJetType = dpmJetGenHeader->ProcessType();
1252   fgLastProcessType = dpmJetType;
1253   MCProcessType globalType = kInvalidProcess;  
1254
1255
1256   if (adebug) {
1257     printf(" AliAnalysisHelperJetTasks::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType);
1258   }
1259
1260
1261   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
1262     globalType = kND;
1263   }  
1264   else if (dpmJetType==5 || dpmJetType==6) {
1265     globalType = kSD;
1266   }
1267   else if (dpmJetType==7) {
1268     globalType = kDD;
1269   }
1270   return globalType;
1271 }
1272