ee12d84191e2f8ca48d4f3b8ee000da85eaf2358
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCore.cxx
1
2 // ******************************************
3 // This task computes several jet observables like 
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and 
6 // correlation strength distribution of particles inside jets.    
7 // Author: lcunquei@cern.ch
8 // *******************************************
9
10
11 /**************************************************************************
12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13  *                                                                        *
14  * Author: The ALICE Off-line Project.                                    *
15  * Contributors are mentioned in the code where appropriate.              *
16  *                                                                        *
17  * Permission to use, copy, modify and distribute this software and its   *
18  * documentation strictly for non-commercial purposes is hereby granted   *
19  * without fee, provided that the above copyright notice appears in all   *
20  * copies and that both the copyright notice and this permission notice   *
21  * appear in the supporting documentation. The authors make no claims     *
22  * about the suitability of this software for any purpose. It is          *
23  * provided "as is" without express or implied warranty.                  *
24  **************************************************************************/
25
26
27 #include "TChain.h"
28 #include "TTree.h"
29 #include "TMath.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH3F.h"
33 #include "THnSparse.h"
34 #include "TCanvas.h"
35
36 #include "AliLog.h"
37
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliVEvent.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliCentrality.h"
45 #include "AliAnalysisHelperJetTasks.h"
46 #include "AliInputEventHandler.h"
47 #include "AliAODJetEventBackground.h"
48 #include "AliAnalysisTaskFastEmbedding.h"
49 #include "AliAODEvent.h"
50 #include "AliAODHandler.h"
51 #include "AliAODJet.h"
52
53 #include "AliAnalysisTaskJetCore.h"
54
55 ClassImp(AliAnalysisTaskJetCore)
56
57 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
58 AliAnalysisTaskSE(),
59 fESD(0x0),
60 fAOD(0x0),
61 fAODExtension(0x0),
62 fBackgroundBranch(""),
63 fNonStdFile(""),
64 fIsPbPb(kTRUE),
65 fOfflineTrgMask(AliVEvent::kAny),
66 fMinContribVtx(1),
67 fVtxZMin(-10.),
68 fVtxZMax(10.),
69 fEvtClassMin(0),
70 fEvtClassMax(4),
71 fFilterMask(0),
72 fRadioFrac(0.2),
73 fMinDist(0.1),
74 fCentMin(0.),
75 fCentMax(100.),
76 fNInputTracksMin(0),
77 fNInputTracksMax(-1),
78 fAngStructCloseTracks(0),
79 fCheckMethods(0),
80 fDoEventMixing(0), 
81 fJetEtaMin(-.5),
82 fJetEtaMax(.5),
83 fNevents(0x0),
84 fTindex(0x0),
85 fTrigBufferIndex(0x0),
86 fJetPtMin(20.),
87 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
88 fJetPtFractionMin(0.5),
89 fNMatchJets(4),
90 fMatchMaxDist(0.8),
91 fKeepJets(kFALSE),
92 fkNbranches(2),
93 fkEvtClasses(12),
94 fOutputList(0x0),
95 fbEvent(kTRUE),
96 fHistEvtSelection(0x0),
97 fhnDeltaR(0x0),
98 fhnMixedEvents(0x0),
99 fh2JetCoreMethod1C10(0x0),
100 fh2JetCoreMethod2C10(0x0),
101 fh2JetCoreMethod1C20(0x0),
102 fh2JetCoreMethod2C20(0x0),
103 fh2JetCoreMethod1C30(0x0),
104 fh2JetCoreMethod2C30(0x0),
105 fh2JetCoreMethod1C60(0x0),
106 fh2JetCoreMethod2C60(0x0),
107 fh2AngStructpt1C10(0x0),
108 fh2AngStructpt2C10(0x0),
109 fh2AngStructpt3C10(0x0),
110 fh2AngStructpt4C10(0x0),
111 fh2AngStructpt1C20(0x0),
112 fh2AngStructpt2C20(0x0),
113 fh2AngStructpt3C20(0x0),
114 fh2AngStructpt4C20(0x0),    
115 fh2AngStructpt1C30(0x0),
116 fh2AngStructpt2C30(0x0),
117 fh2AngStructpt3C30(0x0),
118 fh2AngStructpt4C30(0x0),   
119 fh2AngStructpt1C60(0x0),
120 fh2AngStructpt2C60(0x0),
121 fh2AngStructpt3C60(0x0),
122 fh2AngStructpt4C60(0x0),
123 fh3spectriggered(0x0),
124 fh3specbiased(0x0),
125 fh3spectot(0x0)
126
127  
128 {
129    // default Constructor
130
131
132  // Trigger buffer.
133    for(Int_t i=0; i<10; i++) {
134                 for(Int_t j=0; j<7; j++) {
135                         fTrigBuffer[i][j]=0;
136                 }                               
137         }       
138
139
140
141
142
143    fJetBranchName[0] = "";
144    fJetBranchName[1] = "";
145
146    fListJets[0] = new TList;
147    fListJets[1] = new TList;
148 }
149
150 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
151 AliAnalysisTaskSE(name),
152 fESD(0x0),
153 fAOD(0x0),
154 fAODExtension(0x0),
155 fBackgroundBranch(""),
156 fNonStdFile(""),
157 fIsPbPb(kTRUE),
158 fOfflineTrgMask(AliVEvent::kAny),
159 fMinContribVtx(1),
160 fVtxZMin(-10.),
161 fVtxZMax(10.),
162 fEvtClassMin(0),
163 fEvtClassMax(4),
164 fFilterMask(0),
165 fRadioFrac(0.2),
166 fMinDist(0.1),
167 fCentMin(0.),
168 fCentMax(100.),
169 fNInputTracksMin(0),
170 fNInputTracksMax(-1),
171 fAngStructCloseTracks(0),
172 fCheckMethods(0),
173 fDoEventMixing(0),
174 fJetEtaMin(-.5),
175 fJetEtaMax(.5),
176 fNevents(0x0),
177 fTindex(0x0),
178 fTrigBufferIndex(0x0),
179 fJetPtMin(20.),
180 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
181 fJetPtFractionMin(0.5),
182 fNMatchJets(4),
183 fMatchMaxDist(0.8),
184 fKeepJets(kFALSE),
185 fkNbranches(2),
186 fkEvtClasses(12),
187 fOutputList(0x0),
188 fbEvent(kTRUE),
189 fHistEvtSelection(0x0),
190 fhnDeltaR(0x0),
191 fhnMixedEvents(0x0),
192 fh2JetCoreMethod1C10(0x0),
193 fh2JetCoreMethod2C10(0x0),
194 fh2JetCoreMethod1C20(0x0),
195 fh2JetCoreMethod2C20(0x0),
196 fh2JetCoreMethod1C30(0x0),
197 fh2JetCoreMethod2C30(0x0),
198 fh2JetCoreMethod1C60(0x0),
199 fh2JetCoreMethod2C60(0x0),
200 fh2AngStructpt1C10(0x0),
201 fh2AngStructpt2C10(0x0),
202 fh2AngStructpt3C10(0x0),
203 fh2AngStructpt4C10(0x0),
204 fh2AngStructpt1C20(0x0),
205 fh2AngStructpt2C20(0x0),
206 fh2AngStructpt3C20(0x0),
207 fh2AngStructpt4C20(0x0),    
208 fh2AngStructpt1C30(0x0),
209 fh2AngStructpt2C30(0x0),
210 fh2AngStructpt3C30(0x0),
211 fh2AngStructpt4C30(0x0),   
212 fh2AngStructpt1C60(0x0),
213 fh2AngStructpt2C60(0x0),
214 fh2AngStructpt3C60(0x0),
215 fh2AngStructpt4C60(0x0),    
216 fh3spectriggered(0x0),
217 fh3specbiased(0x0),
218 fh3spectot(0x0)
219
220  {
221    // Constructor
222
223
224     for(Int_t i=0; i<10; i++) {
225        for(Int_t j=0; j<7; j++) {
226             fTrigBuffer[i][j]=0;
227                 }                               
228     }   
229
230
231
232    fJetBranchName[0] = "";
233    fJetBranchName[1] = "";
234
235    fListJets[0] = new TList;
236    fListJets[1] = new TList;
237
238    DefineOutput(1, TList::Class());
239 }
240
241 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
242 {
243    delete fListJets[0];
244    delete fListJets[1];
245 }
246
247 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
248 {
249    fJetBranchName[0] = branch1;
250    fJetBranchName[1] = branch2;
251 }
252
253 void AliAnalysisTaskJetCore::Init()
254 {
255
256    // check for jet branches
257    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
258       AliError("Jet branch name not set.");
259    }
260
261 }
262
263 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
264 {
265    // Create histograms
266    // Called once
267    OpenFile(1);
268    if(!fOutputList) fOutputList = new TList;
269    fOutputList->SetOwner(kTRUE);
270
271    Bool_t oldStatus = TH1::AddDirectoryStatus();
272    TH1::AddDirectory(kFALSE);
273
274
275    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
276    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
277    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
278    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
279    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
280    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
281    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
282
283      UInt_t entries = 0; // bit coded, see GetDimParams() below 
284      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7; 
285      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
286    
287     
288      UInt_t cifras = 0; // bit coded, see GetDimParams() below 
289      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7; 
290      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
291
292     if(fCheckMethods){
293
294
295     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
296     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
297     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
298     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
299     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
300     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
301     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
302     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
303
304    
305     if(fAngStructCloseTracks>0){
306     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
307     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
308     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
309     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
310     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
311     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
312     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
313     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
314     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
315     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
316     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
317     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
318     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
319     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
320     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
321     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
322     fh3spectriggered = new TH3F("Triggered spectrum","",10,0,100,50,0.,200,50,0.,50.);
323     fh3specbiased = new TH3F("Biased spectrum","",10,0,100,50,0.,200.,50,0.,50.);
324     fh3spectot = new TH3F("Total spectrum","",50,0.,200.,50,0.,50.,50,0.,50.);
325
326    fOutputList->Add(fHistEvtSelection);
327
328    fOutputList->Add(fhnDeltaR);
329    
330    fOutputList->Add(fhnMixedEvents);
331          
332      
333   
334       if(fCheckMethods){
335
336       fOutputList->Add(fh2JetCoreMethod1C10);
337       fOutputList->Add(fh2JetCoreMethod2C10);
338       fOutputList->Add(fh2JetCoreMethod1C20);
339       fOutputList->Add(fh2JetCoreMethod2C20);
340       fOutputList->Add(fh2JetCoreMethod1C30);
341       fOutputList->Add(fh2JetCoreMethod2C30);
342       fOutputList->Add(fh2JetCoreMethod1C60);
343       fOutputList->Add(fh2JetCoreMethod2C60);}
344       
345       
346      
347
348
349         if(fAngStructCloseTracks>0){
350        fOutputList->Add(fh2AngStructpt1C10);
351        fOutputList->Add(fh2AngStructpt2C10);
352        fOutputList->Add(fh2AngStructpt3C10);
353        fOutputList->Add(fh2AngStructpt4C10); 
354        fOutputList->Add(fh2AngStructpt1C20);
355        fOutputList->Add(fh2AngStructpt2C20);
356        fOutputList->Add(fh2AngStructpt3C20);
357        fOutputList->Add(fh2AngStructpt4C20); 
358        fOutputList->Add(fh2AngStructpt1C30);
359        fOutputList->Add(fh2AngStructpt2C30);
360        fOutputList->Add(fh2AngStructpt3C30);
361        fOutputList->Add(fh2AngStructpt4C30);
362        fOutputList->Add(fh2AngStructpt1C60);
363        fOutputList->Add(fh2AngStructpt2C60);
364        fOutputList->Add(fh2AngStructpt3C60);
365        fOutputList->Add(fh2AngStructpt4C60);}  
366        fOutputList->Add(fh3spectriggered);
367        fOutputList->Add(fh3specbiased);
368        fOutputList->Add(fh3spectot);
369
370    // =========== Switch on Sumw2 for all histos ===========
371    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
372       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
373       if (h1){
374          h1->Sumw2();
375          continue;
376       }
377     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
378       if (hn){
379          hn->Sumw2();
380       }   
381    }
382    TH1::AddDirectory(oldStatus);
383
384    PostData(1, fOutputList);
385 }
386
387 void AliAnalysisTaskJetCore::UserExec(Option_t *)
388 {
389    
390
391    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
392       AliError("Jet branch name not set.");
393       return;
394    }
395
396    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
397    if (!fESD) {
398       AliError("ESD not available");
399       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
400    } else {
401       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
402    }
403  
404     if(fNonStdFile.Length()!=0){
405     // case that we have an AOD extension we need can fetch the jets from the extended output
406     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
407     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
408     if(!fAODExtension){
409       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
410     }}
411     
412
413
414
415
416    // -- event selection --
417    fHistEvtSelection->Fill(1); // number of events before event selection
418
419    // physics selection
420    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
421    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
422    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
423    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
424       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
425       fHistEvtSelection->Fill(2);
426       PostData(1, fOutputList);
427       return;
428    }
429
430    // vertex selection
431    if(!fAOD){
432      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
433      fHistEvtSelection->Fill(3);
434       PostData(1, fOutputList);
435    }
436    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
437
438    if(!primVtx){
439      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
440      fHistEvtSelection->Fill(3);
441      PostData(1, fOutputList);
442      return;
443    }
444
445    Int_t nTracksPrim = primVtx->GetNContributors();
446    if ((nTracksPrim < fMinContribVtx) ||
447          (primVtx->GetZ() < fVtxZMin) ||
448          (primVtx->GetZ() > fVtxZMax) ){
449       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
450       fHistEvtSelection->Fill(3);
451       PostData(1, fOutputList);
452       return;
453    }
454
455    // event class selection (from jet helper task)
456    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
457    if(fDebug) Printf("Event class %d", eventClass);
458    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
459       fHistEvtSelection->Fill(4);
460       PostData(1, fOutputList);
461       return;
462    }
463
464    // centrality selection
465    AliCentrality *cent = 0x0;
466    Double_t centValue = 0.; 
467    if(fESD) cent = fESD->GetCentrality();
468    if(cent) centValue = cent->GetCentralityPercentile("V0M");
469    if(fDebug) printf("centrality: %f\n", centValue);
470    if (centValue < fCentMin || centValue > fCentMax){
471       fHistEvtSelection->Fill(4);
472       PostData(1, fOutputList);
473       return;
474    }
475
476
477    fHistEvtSelection->Fill(0); 
478    // accepted events  
479    // -- end event selection --
480   
481    // get background
482    AliAODJetEventBackground* externalBackground = 0;
483    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
484       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
485       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
486    }
487    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
488      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
489       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
490    }
491    
492    Float_t rho = 0;
493    if(externalBackground)rho = externalBackground->GetBackground(0);
494
495
496    // fetch jets
497    TClonesArray *aodJets[2];
498    aodJets[0]=0;
499    if(fAOD&&!aodJets[0]){
500    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); 
501    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
502    if(fAODExtension && !aodJets[0]){ 
503    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); 
504    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
505
506    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
507    //Int_t inord[aodJets[0]->GetEntriesFast()];
508    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
509    //  ptsub[n]=0;
510    //  inord[n]=0;}   
511
512    TList ParticleList;
513    Int_t nT = GetListOfTracks(&ParticleList);
514      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
515       fListJets[iJetType]->Clear();
516       if (!aodJets[iJetType]) continue;
517
518       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
519       
520    
521       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
522          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
523          if (jet) fListJets[iJetType]->Add(jet);
524          // if(iJetType==0){
525          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
526       }}
527    
528    Double_t etabig=0;
529    Double_t ptbig=0;
530    Double_t areabig=0;
531    Double_t phibig=0.;
532    Double_t etasmall=0;
533    Double_t ptsmall=0;
534    Double_t areasmall=0;
535    //Double_t distr=0.;
536    Double_t phismall=0.;
537          
538   
539    // Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
540    // Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
541    // Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
542    // Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
543    // Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
544    // Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
545    // Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
546    // Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
547    Int_t iCount=0; 
548    Int_t trigJet=-1;
549    Int_t trigBBTrack=-1;
550    Int_t trigInTrack=-1;
551      
552    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
553            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
554            etabig  = jetbig->Eta();
555            phibig  = jetbig->Phi();
556            ptbig   = jetbig->Pt();
557            if(ptbig==0) continue; 
558            areabig = jetbig->EffectiveAreaCharged();
559            Double_t ptcorr=ptbig-rho*areabig;
560            if(ptcorr<=0) continue;
561            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
562                    Double_t dismin=100.;
563                    Double_t ptmax=-10.; 
564                    Int_t index1=-1;
565                    Int_t index2=-1;
566                   
567                Int_t point=GetHardestTrackBackToJet(jetbig);    
568                AliVParticle *partback = (AliVParticle*)ParticleList.At(point);                            
569                if(!partback) continue; 
570                    fh3spectriggered->Fill(centValue,ptcorr,partback->Pt());
571                        AliAODTrack* leadtrack; 
572                        Int_t ippt=0;
573                        Double_t ppt=-10;   
574                        TRefArray *genTrackList = jetbig->GetRefTracks();
575                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
576                        AliAODTrack* genTrack;
577                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
578                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
579                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
580                        ippt=ir;}}
581                         leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
582                         if(!leadtrack) continue;
583                         fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
584                         if(centValue<10)fh3spectot->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
585                         //store one trigger info                   
586                         if((partback->Pt()>10.)&&(iCount==0)){                        
587                           trigJet=i;
588                           trigBBTrack=point;
589                           trigInTrack=ippt;
590                           iCount=iCount+1;} 
591
592
593                   if(fCheckMethods){
594                   for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
595                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
596                   etasmall  = jetsmall->Eta();
597                   phismall = jetsmall->Phi();
598                   ptsmall   = jetsmall->Pt();
599                   areasmall = jetsmall->EffectiveAreaCharged();
600                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
601                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
602                      //Fraction in the jet core  
603                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
604                     index2=j;}  
605                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
606                       index1=j;}} //en of loop over R=0.2 jets
607                 //method1:most concentric jet=core 
608                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
609                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
610                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
611                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
612                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
613                 //method2:hardest contained jet=core   
614                 if(index2!=-1){ 
615                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
616                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
617                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
618                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
619                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}  
620
621           for(int it = 0;it<nT;++it){
622           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
623           Double_t deltaR = jetbig->DeltaR(part);
624           Double_t deltaEta = etabig-part->Eta();
625           Double_t deltaPhi=phibig-part->Phi();
626           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
627           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
628           Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()};                     fhnDeltaR->Fill(jetEntries);
629           }
630           //end of track loop
631           }
632           //end of jet loop
633
634
635
636
637           if(fDoEventMixing){
638             //check before if the trigger exists
639             // fTrigBuffer[i][0] = zvtx
640             // fTrigBuffer[i][1] = phi
641             // fTrigBuffer[i][2] = eta
642             // fTrigBuffer[i][3] = pt_jet
643             // fTrigBuffer[i][4] = pt_trig
644             // fTrigBuffer[i][5]= pt_track_in
645             // fTrigBuffer[i][6]= centrality
646             if(fTindex==11) fTindex=0;
647             if(fTrigBuffer[fTindex][3]>0){
648             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
649             if (TMath::Abs(fTrigBuffer[fTindex][6]-centValue<10)){  
650               
651                         for(int it = 0;it<nT;++it){
652                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);                        Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
653                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
654                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
655                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
656                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
657                         Double_t triggerEntries[8] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4],fTrigBuffer[fTindex][5]};                      
658                         fhnMixedEvents->Fill(triggerEntries);
659                         }
660                         fNevents=fNevents+1;  
661                         if(fNevents==9) {fTindex=fTindex+1;
662                         fNevents=0;} 
663             }}}
664         
665
666                // Copy the triggers from the current event into the buffer.
667                //again, only if the trigger exists:
668                 if(trigJet>-1){
669                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                 
670                 AliVParticle *partL = (AliVParticle*)ParticleList.At(trigInTrack);
671                 AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);         
672                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
673                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
674                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
675                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
676                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
677                 fTrigBuffer[fTrigBufferIndex][5] = partL->Pt();
678                 fTrigBuffer[fTrigBufferIndex][6] = centValue;
679                 fTrigBufferIndex++;
680                 if(fTrigBufferIndex==9) fTrigBufferIndex=0;
681                 }
682           }
683           
684
685
686
687
688      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
689
690      //tracks up to R=0.8 distant from the jet axis
691    //   if(fAngStructCloseTracks==1){
692    //    TList CloseTrackList;
693    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
694    //    Double_t difR=0.04;
695    //    for(Int_t l=0;l<15;l++){
696    //    Double_t rr=l*0.1+0.1;
697    //     for(int it = 0;it<nn;++it){
698    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
699    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
700    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
701    //         Double_t ptm=part1->Pt();
702    //         Double_t ptn=part2->Pt(); 
703    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
704    //                    Rnm=TMath::Sqrt(Rnm);
705    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
706    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
707    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
708    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
709    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
710    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
711    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
712    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
713    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
714    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
715    //   }
716     
717    //   //only jet constituents
718    //    if(fAngStructCloseTracks==2){
719
720    //    Double_t difR=0.04;
721    //    for(Int_t l=0;l<15;l++){
722    //    Double_t rr=l*0.1+0.1;
723
724     
725    //    AliAODTrack* part1;
726    //    AliAODTrack* part2;
727
728    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
729    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
730           
731              
732
733    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
734    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
735    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
736    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
737    //         Double_t ptm=part1->Pt();
738    //         Double_t ptn=part2->Pt(); 
739    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
740    //                    Rnm=TMath::Sqrt(Rnm);
741    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
742    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
743    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
744    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
745    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
746    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
747    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
748    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
749    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
750    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
751    // }
752      // //end loop over R=0.4 jets      
753      // if(fAngStructCloseTracks>0){
754      // for(Int_t l=0;l<15;l++){
755      // Double_t rr=l*0.1+0.1;
756      //    if(down1[l]!=0){  
757      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
758      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
759      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
760      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
761      //    if(down2[l]!=0){  
762      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
763      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
764      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
765      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
766      //    if(down3[l]!=0){  
767      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
768      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
769      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
770      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
771      //    if(down4[l]!=0){  
772      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
773      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
774      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
775      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
776
777     
778
779
780
781
782
783    PostData(1, fOutputList);
784 }
785
786 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
787 {
788    // Draw result to the screen
789    // Called once at the end of the query
790
791    if (!GetOutputData(1))
792    return;
793 }
794
795
796
797
798
799
800
801
802
803
804
805 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
806
807     Int_t iCount = 0;
808  
809     
810     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
811       AliAODTrack *tr = fAOD->GetTrack(it);
812       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
813       if(TMath::Abs(tr->Eta())>0.9)continue;
814       if(tr->Pt()<0.15)continue;
815       list->Add(tr);
816       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
817       iCount++;
818     }
819   
820    
821   return iCount;
822  
823 }
824
825    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
826
827    
828     Int_t index=-1;
829     Double_t ptmax=-10;
830     Double_t dphi=0;
831     Double_t dif=0;
832     Int_t iCount=0;
833     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
834       AliAODTrack *tr = fAOD->GetTrack(it);
835       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
836       if(TMath::Abs(tr->Eta())>0.9)continue;
837       if(tr->Pt()<0.15)continue;
838       iCount=iCount+1;
839       dphi=RelativePhi(tr->Phi(),jetbig->Phi());  
840       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
841       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
842       index=iCount-1;
843       dif=dphi;  }}
844   
845       return index;
846
847    }
848
849
850
851
852
853
854
855
856
857  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
858
859     Int_t iCount = 0;
860  
861   
862     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
863       AliAODTrack *tr = fAOD->GetTrack(it);
864       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
865       if(TMath::Abs(tr->Eta())>0.9)continue;
866       if(tr->Pt()<0.15)continue;
867       Double_t disR=jetbig->DeltaR(tr);
868       if(disR>0.8)  continue;
869       list->Add(tr);
870       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
871       iCount++;
872     }
873   
874    list->Sort();
875    return iCount;
876
877 }
878
879
880
881
882
883
884
885
886
887
888
889 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
890 {
891
892    Int_t nInputTracks = 0;
893
894    TString jbname(fJetBranchName[1]);
895    //needs complete event, use jets without background subtraction
896    for(Int_t i=1; i<=3; ++i){
897       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
898    }
899    // use only HI event
900    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
901    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
902
903    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
904    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
905    if(!tmpAODjets){
906       Printf("Jet branch %s not found", jbname.Data());
907       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
908       return -1;
909    }
910    
911    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
912       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
913       if(!jet) continue;
914       TRefArray *trackList = jet->GetRefTracks();
915       Int_t nTracks = trackList->GetEntriesFast();
916       nInputTracks += nTracks;
917       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
918    }
919    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
920
921    return nInputTracks;  
922 }
923
924
925
926 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
927
928   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
929   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
930   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
931   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
932   double dphi = mphi-vphi;
933   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
934   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
935   return dphi;//dphi in [-Pi, Pi]
936 }
937
938
939
940 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
941 {
942    // generate new THnSparseF, axes are defined in GetDimParams()
943
944    Int_t count = 0;
945    UInt_t tmp = entries;
946    while(tmp!=0){
947       count++;
948       tmp = tmp &~ -tmp;  // clear lowest bit
949    }
950
951    TString hnTitle(name);
952    const Int_t dim = count;
953    Int_t nbins[dim];
954    Double_t xmin[dim];
955    Double_t xmax[dim];
956
957    Int_t i=0;
958    Int_t c=0;
959    while(c<dim && i<32){
960       if(entries&(1<<i)){
961       
962          TString label("");
963          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
964          hnTitle += Form(";%s",label.Data());
965          c++;
966       }
967       
968       i++;
969    }
970    hnTitle += ";";
971
972    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
973 }
974
975 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
976 {
977    // stores label and binning of axis for THnSparse
978
979    const Double_t pi = TMath::Pi();
980    
981    switch(iEntry){
982       
983    case 0:
984       label = "V0 centrality (%)";
985      
986          nbins = 10;
987          xmin = 0.;
988          xmax = 100.;
989          break;
990       
991       
992    case 1:
993       label = "corrected jet pt";
994          nbins = 50;
995          xmin = 0.;
996          xmax = 200.;
997           break;
998       
999       
1000    case 2:
1001       label = "track pT";
1002      
1003          nbins = 1000;
1004          xmin = 0.;
1005          xmax = 50;
1006          break;
1007       
1008       
1009     case 3:
1010       label = "deltaR";
1011       nbins = 15;
1012       xmin = 0.;
1013       xmax = 1.5;
1014       break;
1015
1016
1017
1018    case 4:
1019       label = "deltaEta";
1020       nbins = 30;
1021       xmin = -1.5;
1022       xmax = 1.5;
1023       break;
1024
1025
1026   case 5:
1027       label = "deltaPhi";
1028       nbins = 90;
1029       xmin = -0.5*pi;
1030       xmax = 1.5*pi;
1031       break;   
1032    
1033       
1034         
1035     case 6:
1036       label = "leading track";
1037       nbins = 50;
1038       xmin = 0;
1039       xmax = 50;
1040       break;
1041            
1042      case 7:
1043     
1044       label = "trigger track";
1045       nbins =50;
1046       xmin = 0;
1047       xmax = 50;
1048       break;
1049
1050
1051    
1052   
1053
1054
1055
1056
1057    }
1058
1059 }
1060