Jet and Particle identification tasks moved to different directories
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskUE.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id:$ */
17  
18 #include <TROOT.h>
19 #include <TSystem.h>
20 #include <TChain.h>
21 #include <TFile.h>
22 #include <TList.h>
23 #include <TH1I.h>
24 #include <TH1F.h>
25 #include <TH2F.h>
26 #include <TCanvas.h>
27 #include <TVector3.h>
28 #include <TLorentzVector.h>
29 #include <TMath.h>
30
31 #include "AliAnalysisTaskUE.h"
32 #include "AliAnalysisManager.h"
33 #include "AliMCEventHandler.h"
34 #include "AliAODEvent.h"
35 #include "AliAODInputHandler.h"
36 #include "AliAODHandler.h"
37 #include "AliStack.h"
38 #include "AliAODJet.h"
39 #include "AliAODTrack.h"
40
41 #include "AliLog.h"
42
43 //
44 // Analysis class for Underlying Event studies
45 //
46 // Look for correlations on the tranverse regions to 
47 // the leading charged jet
48 //
49 // This class needs as input AOD with track and Jets
50 // the output is a list of histograms
51 //
52 // AOD can be either connected to the InputEventHandler  
53 // for a chain of AOD files 
54 // or 
55 // to the OutputEventHandler
56 // for a chain of ESD files, so this case class should be 
57 // in the train after the Jet finder
58 //
59 //    Arian.Abrahantes.Quintana@cern.ch 
60 //    Ernesto.Lopez.Torres@cern.ch
61 // 
62
63
64 ClassImp( AliAnalysisTaskUE)
65
66 ////////////////////////////////////////////////////////////////////////
67
68
69 //____________________________________________________________________
70 AliAnalysisTaskUE:: AliAnalysisTaskUE(const char* name):
71   AliAnalysisTask(name, ""),
72   fDebug(kFALSE),          
73   fAOD(0x0),            
74   fListOfHistos(0x0),  
75   fBinsPtInHist(30),     
76   fMinJetPtInHist(0.),
77   fMaxJetPtInHist(300.),  
78   fAnaType(1),         
79   fRegionType(1),
80   fConeRadius(0.7),
81   fOrdering(1),     
82   fJet1EtaCut(0.2),
83   fJet2DeltaPhiCut(2.616),    // 150 degrees
84   fJet2RatioPtCut(0.8),
85   fJet3PtCut(15.),
86   fTrackPtCut(0.),
87   fTrackEtaCut(0.9),
88   fhNJets(0x0),
89   fhEleadingPt(0x0),
90   fhMinRegPtDist(0x0),
91   fhRegionMultMin(0x0),
92   fhMinRegAvePt(0x0), 
93   fhMinRegSumPt(0x0),            
94   fhMinRegMaxPtPart(0x0),
95   fhMinRegSumPtvsMult(0x0),
96   fhdNdEta_PhiDist(0x0),        
97   fhFullRegPartPtDistVsEt(0x0), 
98   fhTransRegPartPtDistVsEt(0x0),
99   fhRegionSumPtMaxVsEt(0x0),
100   fhRegionMultMax(0x0),         
101   fhRegionMultMaxVsEt(0x0),     
102   fhRegionSumPtMinVsEt(0x0), //fhRegionMultMin(0x0),         
103   fhRegionMultMinVsEt(0x0),     
104   fhRegionAveSumPtVsEt(0x0),    
105   fhRegionDiffSumPtVsEt(0x0),   
106   fhRegionAvePartPtMaxVsEt(0x0),
107   fhRegionAvePartPtMinVsEt(0x0),
108   fhRegionMaxPartPtMaxVsEt(0x0)//,   fhValidRegion(0x0)
109 {
110   // Default constructor
111   // Define input and output slots here
112   // Input slot #0 works with a TChain
113   DefineInput(0, TChain::Class());
114   // Output slot #0 writes into a TList container
115   DefineOutput(0, TList::Class());
116 }
117
118
119 //____________________________________________________________________
120 void AliAnalysisTaskUE::ConnectInputData(Option_t* /*option*/)
121 {
122 // Connect the input data  
123
124 // We need AOD with tracks and jets.
125 // Since AOD can be either connected to the InputEventHandler or to the OutputEventHandler
126 // we need to check where it is and get the pointer to AODEvent in the right way
127   
128   if (fDebug > 1) AliInfo("ConnectInputData() \n");
129   
130   
131   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
132   
133   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
134      fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
135   } else {
136       handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
137       if( handler && handler->InheritsFrom("AliAODHandler") ) {
138          fAOD  = ((AliAODHandler*)handler)->GetAOD();
139       } else {
140         AliFatal("I can't get any AOD Event Handler");
141       }
142    }
143 }
144
145 //____________________________________________________________________
146 void  AliAnalysisTaskUE::CreateOutputObjects()
147 {
148 // Create the output container
149 //
150   if (fDebug > 1) AliInfo("CreateOutPutData()");
151   //
152   //  Histograms
153   CreateHistos();
154   fListOfHistos->SetOwner(kTRUE);  
155 //  OpenFile(0);
156 }
157
158
159 //____________________________________________________________________
160 void  AliAnalysisTaskUE::Exec(Option_t */*option*/)
161 {
162 // Execute analysis for current event
163 //
164   AnalyseUE();   
165   // Post the data
166   PostData(0, fListOfHistos);
167 }
168   
169   
170 //____________________________________________________________________
171 void  AliAnalysisTaskUE::Terminate(Option_t */*option*/)
172 {
173 // Terminate analysis
174 //
175
176   // Normalize histos to region area TODO: 
177   Double_t areafactor = 1.;
178   if( fRegionType == 1 ) {
179     areafactor = 1./ (fTrackEtaCut *2. * TMath::Pi()*2.);
180   } else {
181     areafactor = 1./ ( fConeRadius * fConeRadius * TMath::Pi() );
182   }
183   
184   // Draw some Test plot to the screen
185   if (!gROOT->IsBatch()) {
186     
187     // Update pointers reading them from the output slot
188     fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));
189     if( !fListOfHistos ) {
190       AliError("Histogram List is not available");
191       return;
192     }
193     fhEleadingPt         = (TH1F*)fListOfHistos->At(1);
194     fhdNdEta_PhiDist     = (TH1F*)fListOfHistos->At(8);
195     fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11);
196     fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12);
197     fhRegionMultMaxVsEt  = (TH1F*)fListOfHistos->At(14);
198     fhRegionMultMinVsEt  = (TH1F*)fListOfHistos->At(15);
199     fhRegionAveSumPtVsEt = (TH1F*)fListOfHistos->At(16);
200
201     fhNJets              = (TH1F*)fListOfHistos->At(0);
202
203     // Canvas
204     TCanvas* c1 = new TCanvas("c1",Form("sumPt dist (%s)", GetTitle()),60,60,1200,800);
205     c1->Divide(2,2);
206     c1->cd(1);
207     TH1F *h1r = new TH1F("hRegionEtvsSumPtMax" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
208     h1r->Divide(fhRegionSumPtMaxVsEt,fhEleadingPt,1,1);
209     h1r->Scale( areafactor );
210     h1r->SetMarkerStyle(20);
211     h1r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
212     h1r->SetYTitle("P_{T}^{90, max}");
213     h1r->DrawCopy("p");
214     
215     c1->cd(2);
216     
217     TH1F *h2r = new TH1F("hRegionEtvsSumPtMin" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
218     h2r->Divide(fhRegionSumPtMinVsEt,fhEleadingPt,1,1);
219     h2r->Scale( areafactor );
220     h2r->SetMarkerStyle(20);
221     h2r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
222     h2r->SetYTitle("P_{T}^{90, min}");
223     h2r->DrawCopy("p");
224
225     c1->cd(3);
226     TH1F *h4r = new TH1F("hRegionEtvsDiffPt" , "", fBinsPtInHist,  fMinJetPtInHist, fMaxJetPtInHist);
227     h4r->Divide(fhRegionAveSumPtVsEt,fhEleadingPt,1,1);
228     h4r->Scale(2.); // make average
229     h4r->Scale( areafactor );
230     h4r->SetYTitle("#DeltaP_{T}^{90}");
231     h4r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
232     h4r->SetMarkerStyle(20);
233     h4r->DrawCopy("p");
234     
235     c1->cd(4);
236     TH1F *h5r = new TH1F("hRegionMultMaxVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
237     TH1F *h6r = new TH1F("hRegionMultMinVsEtleading",   "",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
238     
239     h5r->Divide(fhRegionMultMaxVsEt,fhEleadingPt,1,1);
240     h6r->Divide(fhRegionMultMinVsEt,fhEleadingPt,1,1);
241     h5r->Scale( areafactor );
242     h6r->Scale( areafactor );
243     h5r->SetYTitle("N_{Tracks}^{90}");
244     h5r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
245     h5r->SetMarkerStyle(20);
246     h5r->DrawCopy("p");
247     h6r->SetMarkerStyle(21);
248     h6r->SetMarkerColor(2);
249     h6r->SetYTitle("N_{Tracks}^{90}");
250     h6r->SetXTitle("P_{T} of Leading Jet (GeV/c)");
251     h6r->DrawCopy("p same");
252     c1->Update();
253     
254     TCanvas* c2 = new TCanvas("c2","Jet Pt dist",160,160,1200,800);
255     c2->Divide(2,2);
256     c2->cd(1);
257     fhEleadingPt->SetMarkerStyle(20);
258     fhEleadingPt->SetMarkerColor(2);
259     fhEleadingPt->DrawCopy("p");
260     gPad->SetLogy();
261     
262     c2->cd(2);
263     fhdNdEta_PhiDist->SetMarkerStyle(20);
264     fhdNdEta_PhiDist->SetMarkerColor(2);
265     fhdNdEta_PhiDist->DrawCopy("p");
266     gPad->SetLogy();
267     
268     c2->cd(3);      
269     fhNJets->DrawCopy();
270     
271     //fhTransRegPartPtDist  = (TH1F*)fListOfHistos->At(2);
272     fhRegionMultMin          = (TH1F*)fListOfHistos->At(3);
273     fhMinRegAvePt     = (TH1F*)fListOfHistos->At(4);
274     fhMinRegSumPt     = (TH1F*)fListOfHistos->At(5);   
275     //fhMinRegMaxPtPart     = (TH1F*)fListOfHistos->At(6);
276     fhMinRegSumPtvsMult   = (TH1F*)fListOfHistos->At(7);
277
278     // Canvas
279     TCanvas* c3 = new TCanvas("c3"," p_{T} dist",160,160,1200,800);
280     c3->Divide(2,2);
281     c3->cd(1);
282     /*fhTransRegPartPtDist->SetMarkerStyle(20);
283     fhTransRegPartPtDist->SetMarkerColor(2); 
284     fhTransRegPartPtDist->Scale(areafactor/fhTransRegPartPtDist->GetEntries());
285     fhTransRegPartPtDist->DrawCopy("p");
286     gPad->SetLogy();
287     */
288     c3->cd(2); 
289     fhMinRegSumPt->SetMarkerStyle(20);
290     fhMinRegSumPt->SetMarkerColor(2);  
291     fhMinRegSumPt->Scale(areafactor);
292     fhMinRegSumPt->DrawCopy("p");
293     gPad->SetLogy();
294     
295     c3->cd(3);
296     fhMinRegAvePt->SetMarkerStyle(20);
297     fhMinRegAvePt->SetMarkerColor(2);  
298     fhMinRegAvePt->Scale(areafactor);
299     fhMinRegAvePt->DrawCopy("p");
300     gPad->SetLogy();
301     
302     c3->cd(4);
303     
304     TH1F *h7r = new TH1F("hRegionMultMinVsMult",   "",  21, -0.5,   20.5);
305     
306     h7r->Divide(fhMinRegSumPtvsMult,fhRegionMultMin,1,1);
307     
308     h7r->SetMarkerStyle(20);
309     h7r->SetMarkerColor(2);   
310     h7r->DrawCopy("p");
311
312     c3->Update();
313
314     
315 /*    c2->cd(3);      
316     fhValidRegion->SetMarkerStyle(20);
317     fhValidRegion->SetMarkerColor(2);
318     fhValidRegion->DrawCopy("p");
319 */    
320     c2->Update();
321   } else {
322     AliInfo(" Batch mode, not histograms will shown...");
323   }
324   
325   if( fDebug > 1 ) 
326      AliInfo("End analysis");
327   
328 }
329
330
331 //____________________________________________________________________
332 void  AliAnalysisTaskUE::AnalyseUE()
333 {
334   //
335   // Look for correlations on the tranverse regions to 
336   // the leading charged jet
337   // 
338   // ------------------------------------------------
339   // Find Leading Jets 1,2,3 
340   // (could be skipped if Jets are sort by Pt...)
341   Double_t maxPtJet1 = 0.; 
342   Int_t    index1 = -1;
343   Double_t maxPtJet2 = 0.;   // jet 2 need for back to back inclusive
344   Int_t    index2 = -1;
345   Double_t maxPtJet3 = 0.;   // jet 3 need for back to back exclusive
346   Int_t    index3 = -1;
347   Int_t nJets = fAOD->GetNJets();
348   for( Int_t i=0; i<nJets; ++i ) {
349     AliAODJet* jet = fAOD->GetJet(i);
350     if( jet->Pt() > maxPtJet1 ) {
351       maxPtJet3 = maxPtJet2; index3 = index2;
352       maxPtJet2 = maxPtJet1; index2 = index1;
353       maxPtJet1 = jet->Pt(); index1 = i;
354     } else if( jet->Pt() > maxPtJet2 ) {
355       maxPtJet3 = maxPtJet2; index3 = index2;
356       maxPtJet2 = jet->Pt(); index2 = i;
357     } else if( jet->Pt() > maxPtJet3 ) {
358       maxPtJet3 = jet->Pt(); index3 = i;
359     }
360   }
361
362         if( index1 < 0 ) {
363       AliInfo("\nSkipping Event, not jet found...");
364       return;
365         } else
366       AliInfo(Form("\nPt Leading Jet = %6.1f eta=%5.1f ",  maxPtJet1, fAOD->GetJet(index1)->Eta() ));
367   
368   fhNJets->Fill(nJets);
369   
370   TVector3 jetVect[2];
371   
372   // ----------------------------------------------
373   // Cut events by jets topology 
374   // fAnaType:
375   //     1 = inclusive, 
376   //         - Jet1 |eta| < fJet1EtaCut
377   //     2 = back to back inclusive
378   //         - fulfill case 1
379   //         - |Jet1.Phi - Jet2.Phi| > fJet2DeltaPhiCut
380   //         - Jet2.Pt/Jet1Pt > fJet2RatioPtCut
381   //     3 = back to back exclusive   
382   //         - fulfill case 2
383   //         - Jet3.Pt < fJet3PtCut
384   
385   Bool_t isInclusive = kFALSE;
386   
387   if(  TMath::Abs(fAOD->GetJet(index1)->Eta()) > fJet1EtaCut) {
388     if( fDebug > 1 ) AliInfo("Skipping Event...Jet1 |eta| > fJet1EtaCut");
389     return;
390   }
391   isInclusive = kTRUE;
392   jetVect[0].SetXYZ(fAOD->GetJet(index1)->Px(),
393                     fAOD->GetJet(index1)->Py(),
394                     fAOD->GetJet(index1)->Pz());
395   // back to back inclusive
396   Bool_t isB2Binclusive = kFALSE;                  
397   if( fAnaType > 1 && index2 > 0 && isInclusive) {  
398     jetVect[1].SetXYZ(fAOD->GetJet(index2)->Px(),               
399                       fAOD->GetJet(index2)->Py(),
400                       fAOD->GetJet(index2)->Pz());
401     if( TMath::Abs(jetVect[0].DeltaPhi(jetVect[1])) < fJet2DeltaPhiCut ||
402         maxPtJet2/maxPtJet1 < fJet2RatioPtCut ) {
403       if( fDebug > 1 ) AliInfo("Skipping Event... |Jet1.Phi - Jet2.Phi| < fJet2DeltaPhiCut");
404       return;                 
405     }
406     isB2Binclusive = kTRUE;
407   }
408   if (isInclusive && !isB2Binclusive && fAnaType>1) return;
409   // back to back exclusive
410   Bool_t isB2Bexclusive = kFALSE;
411   if( fAnaType > 2 && index3 > 0 && isB2Binclusive) {  
412     if( fAOD->GetJet(index3)->Pt() > fJet3PtCut ) {
413       if( fDebug > 1 ) AliInfo("Skipping Event... Jet3.Pt > fJet3PtCut ");
414       return;
415     }
416     isB2Bexclusive = kTRUE;
417   }
418   if (isB2Binclusive && !isB2Bexclusive && fAnaType>2) return;
419   
420   AliInfo(Form("njets = %d",nJets));
421   fhEleadingPt->Fill( maxPtJet1 );
422
423   // ----------------------------------------------
424   // Find max and min regions
425   Double_t sumPtRegionPosit = 0.;
426   Double_t sumPtRegionNegat = 0.;
427   Double_t maxPartPtRegion  = 0.;
428   Int_t    nTrackRegionPosit = 0;
429   Int_t    nTrackRegionNegat = 0;
430   static const Double_t k270rad = 270.*TMath::Pi()/180.;
431   
432   Int_t nTracks = fAOD->GetNTracks();
433   for (Int_t ipart=0; ipart<nTracks; ++ipart) {
434     AliAODTrack* part = fAOD->GetTrack( ipart );
435       
436     if ( !part->Charge() ) continue; 
437     if ( part->Pt() < fTrackPtCut ) continue;
438
439     TVector3 partVect(part->Px(), part->Py(), part->Pz());
440     
441     Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad;
442     if( deltaPhi > 2.*TMath::Pi() )  deltaPhi-= 2.*TMath::Pi();
443     fhdNdEta_PhiDist->Fill( deltaPhi );
444     fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
445      
446     Int_t region = IsTrackInsideRegion( jetVect, &partVect );  
447     
448     if (region > 0) {
449        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
450        sumPtRegionPosit += part->Pt();
451        nTrackRegionPosit++;
452        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
453     }
454     if (region < 0) {
455        if( maxPartPtRegion < part->Pt() ) maxPartPtRegion = part->Pt();
456        sumPtRegionNegat += part->Pt();
457        nTrackRegionNegat++;
458        fhTransRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 );
459     }
460   }
461
462   //How quantities will be sorted before Fill Min and Max Histogram
463   //  1=Plots will be CDF-like
464   //  2=Plots will be Marchesini-like
465   if (fOrdering==1){
466     if (sumPtRegionPosit > sumPtRegionNegat) {
467       FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat );
468     } else {
469       FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit );
470     }
471     if (nTrackRegionPosit > nTrackRegionNegat ) {
472       FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat );
473     } else {
474       FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit );
475     }
476   } else if (fOrdering==2) {
477     if (sumPtRegionPosit > sumPtRegionNegat) {
478       FillSumPtRegion( maxPtJet1, sumPtRegionPosit, sumPtRegionNegat );
479       FillMultRegion( maxPtJet1, nTrackRegionPosit, nTrackRegionNegat, sumPtRegionNegat );
480     } else {
481       FillSumPtRegion( maxPtJet1, sumPtRegionNegat, sumPtRegionPosit );
482       FillMultRegion( maxPtJet1, nTrackRegionNegat, nTrackRegionPosit, sumPtRegionPosit );
483     }
484   }
485   
486   Double_t avePosRegion = (nTrackRegionPosit) ? sumPtRegionPosit/nTrackRegionPosit : 0.;
487   Double_t aveNegRegion = (nTrackRegionNegat) ? sumPtRegionNegat/nTrackRegionNegat : 0.;
488   if (avePosRegion > aveNegRegion) {
489     FillAvePartPtRegion( maxPtJet1, avePosRegion, aveNegRegion );
490   } else {
491     FillAvePartPtRegion( maxPtJet1, aveNegRegion, avePosRegion );
492   }
493
494   fhRegionMaxPartPtMaxVsEt->Fill(maxPtJet1, maxPartPtRegion );
495        
496   // Compute pedestal like magnitude
497   fhRegionDiffSumPtVsEt->Fill(maxPtJet1, TMath::Abs(sumPtRegionPosit-sumPtRegionNegat)/2.0);
498   fhRegionAveSumPtVsEt->Fill(maxPtJet1, (sumPtRegionPosit+sumPtRegionNegat)/2.0);
499
500 }
501
502
503 //____________________________________________________________________
504 void AliAnalysisTaskUE::FillSumPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
505 {
506    // Fill sumPt of control regions
507    
508    // Max cone
509    fhRegionSumPtMaxVsEt->Fill( leadingE, ptMax );
510    // Min cone
511    fhRegionSumPtMinVsEt->Fill( leadingE, ptMin );
512    // MAke distributions for UE comparison with MB data
513    fhMinRegSumPt->Fill(ptMin);
514 }
515
516 //____________________________________________________________________
517 void AliAnalysisTaskUE::FillAvePartPtRegion( Double_t leadingE, Double_t ptMax, Double_t ptMin  )
518 {
519    // Fill average particle Pt of control regions
520    
521    // Max cone
522    fhRegionAvePartPtMaxVsEt->Fill( leadingE, ptMax );
523    // Min cone
524    fhRegionAvePartPtMinVsEt->Fill( leadingE, ptMin );
525    // MAke distributions for UE comparison with MB data
526    fhMinRegAvePt->Fill(ptMin);
527 }
528
529 //____________________________________________________________________
530 void AliAnalysisTaskUE::FillMultRegion( Double_t leadingE, Double_t nTrackPtmax, Double_t nTrackPtmin, Double_t ptMin  )
531 {
532    // Fill Nch multiplicity of control regions
533    
534    // Max cone
535    fhRegionMultMaxVsEt->Fill( leadingE, nTrackPtmax );
536    fhRegionMultMax->Fill( nTrackPtmax );
537    // Min cone
538    fhRegionMultMinVsEt->Fill( leadingE, nTrackPtmin );
539    fhRegionMultMin->Fill( nTrackPtmin );
540    // MAke distributions for UE comparison with MB data
541    fhMinRegSumPtvsMult->Fill(nTrackPtmin,ptMin);
542  
543 }
544
545 //____________________________________________________________________
546 Int_t AliAnalysisTaskUE::IsTrackInsideRegion(TVector3 *jetVect, TVector3 *partVect) 
547 {  
548   // return de region in delta phi
549   // -1 negative delta phi 
550   //  1 positive delta phi
551   //  0 outside region
552   static const Double_t k60rad  = 60.*TMath::Pi()/180.;
553   static const Double_t k120rad = 120.*TMath::Pi()/180.;
554   
555   Int_t region = 0;
556   if( fRegionType == 1 ) {
557     if( TMath::Abs(partVect->Eta()) > fTrackEtaCut ) return 0;
558     // transverse regions
559     if (jetVect[0].DeltaPhi(*partVect) < -k60rad && jetVect[0].DeltaPhi(*partVect) > -k120rad ) region = -1;
560     if (jetVect[0].DeltaPhi(*partVect) > k60rad && jetVect[0].DeltaPhi(*partVect) < k120rad ) region = 1;
561   
562   } else if( fRegionType == 2 ) {
563     // Cone regions
564     Double_t deltaR = 0.;
565     
566     TVector3 positVect,negatVect;
567     positVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()+TMath::PiOver2());
568     negatVect.SetMagThetaPhi(1, 2.*atan(exp(-jetVect[0].Eta())), jetVect[0].Phi()-TMath::PiOver2());
569     
570     if (TMath::Abs(positVect.DeltaPhi(*partVect)) < fConeRadius ) { 
571        region = 1;  
572        deltaR = positVect.DrEtaPhi(*partVect);
573     } else if (TMath::Abs(negatVect.DeltaPhi(*partVect)) < fConeRadius) { 
574        region = -1;  
575        deltaR = negatVect.DrEtaPhi(*partVect);
576     }
577     
578     if (deltaR > fConeRadius) region = 0;
579   
580   } else 
581     AliError("Unknow region type");
582
583   // For debug (to be removed)
584   //if( region != 0 )  fhValidRegion->Fill( partVect->Eta()-jetVect[0].Eta(), jetVect[0].DeltaPhi(*partVect) );
585   
586   return region;
587 }
588
589 //____________________________________________________________________
590 void  AliAnalysisTaskUE::CreateHistos()
591 {
592   fListOfHistos = new TList();
593   
594   
595   fhNJets = new TH1F("fhNJets", "n Jet",  10, 0, 10);
596   fhNJets->SetXTitle("# of jets");
597   fhNJets->Sumw2();
598   fhNJets->Sumw2();
599   fListOfHistos->Add( fhNJets );                 // At(0) 
600   
601   fhEleadingPt  = new TH1F("hEleadingPt",   "Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
602   fhEleadingPt->SetXTitle("P_{T} (GeV/c)");
603   fhEleadingPt->SetYTitle("dN/dP_{T}");
604   fhEleadingPt->Sumw2();
605   fListOfHistos->Add( fhEleadingPt );            // At(1)
606   
607   fhMinRegPtDist = new TH1F("hMinRegPtDist",   "P_{T} distribution in Min zone",  50,0.,20.);
608   fhMinRegPtDist->SetXTitle("P_{T} (GeV/c)");
609   fhMinRegPtDist->SetYTitle("dN/dP_{T}");
610   fhMinRegPtDist->Sumw2();
611   fListOfHistos->Add( fhMinRegPtDist );          // At(2)
612   
613   fhRegionMultMin = new TH1F("hRegionMultMin",      "N_{ch}^{90, min}",  21, -0.5,   20.5);
614   fhRegionMultMin->SetXTitle("N_{ch tracks}");
615   fhRegionMultMin->Sumw2();
616   fListOfHistos->Add( fhRegionMultMin );         // At(3)            
617   
618   fhMinRegAvePt = new TH1F("hMinRegAvePt", "#LTp_{T}#GT",  50, 0.,   20.);
619   fhMinRegAvePt->SetXTitle("P_{T} (GeV/c)");
620   fhMinRegAvePt->Sumw2();
621   fListOfHistos->Add( fhMinRegAvePt );           // At(4)
622   
623   fhMinRegSumPt = new TH1F("hMinRegSumPt", "#Sigma p_{T} ",  50, 0.,   20.);
624   fhMinRegSumPt->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
625   fhMinRegSumPt->SetXTitle("#Sigma p_{T} (GeV/c)");
626   fhMinRegSumPt->Sumw2();
627   fListOfHistos->Add( fhMinRegSumPt );           // At(5)
628               
629   fhMinRegMaxPtPart = new TH1F("hMinRegMaxPtPart", "max(p_{T})|_{event} ",  50, 0.,   20.);
630   fhMinRegMaxPtPart->SetYTitle("Ed^{3}N_{tracks}/dp^{3} (c^{3}/GeV^{2})");  
631   fhMinRegMaxPtPart->SetXTitle("p_{T} (GeV/c)");
632   fhMinRegMaxPtPart->Sumw2();
633   fListOfHistos->Add( fhMinRegMaxPtPart );       // At(6)
634   
635   fhMinRegSumPtvsMult = new TH1F("hMinRegSumPtvsMult", "#Sigma p_{T} vs. Multiplicity ",  21, -0.5,   20.5);
636   fhMinRegSumPtvsMult->SetYTitle("#Sigma p_{T} (GeV/c)");  
637   fhMinRegSumPtvsMult->SetXTitle("N_{charge}");
638   fhMinRegSumPtvsMult->Sumw2();
639   fListOfHistos->Add( fhMinRegSumPtvsMult );     // At(7);
640
641   fhdNdEta_PhiDist  = new TH1F("hdNdEta_PhiDist",   "Charge particle density |#eta|< 1 vs #Delta#phi",  120, 0.,   2.*TMath::Pi());
642   fhdNdEta_PhiDist->SetXTitle("#Delta#phi");
643   fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi");
644   fhdNdEta_PhiDist->Sumw2();
645   fListOfHistos->Add( fhdNdEta_PhiDist );        // At(8)
646    
647   // Can be use to get part pt distribution for differente Jet Pt bins
648   fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}",  
649                50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
650   fhFullRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
651   fhFullRegPartPtDistVsEt->SetXTitle("p_{T}");
652   fhFullRegPartPtDistVsEt->Sumw2();
653   fListOfHistos->Add( fhFullRegPartPtDistVsEt );  // At(9) 
654    
655    // Can be use to get part pt distribution for differente Jet Pt bins
656   fhTransRegPartPtDistVsEt = new TH2F("hTransRegPartPtDistVsEt", "dN/dP_{T} in tranvese regions |#eta|<1 vs Leading Jet P_{T}",  
657                50,0.,50., fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
658   fhTransRegPartPtDistVsEt->SetYTitle("Leading Jet P_{T}");
659   fhTransRegPartPtDistVsEt->SetXTitle("p_{T}");
660   fhTransRegPartPtDistVsEt->Sumw2();
661   fListOfHistos->Add( fhTransRegPartPtDistVsEt );  // At(10)
662   
663   
664   fhRegionSumPtMaxVsEt = new TH1F("hRegionSumPtMaxVsEt",  "P_{T}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
665   fhRegionSumPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
666   fhRegionSumPtMaxVsEt->Sumw2();
667   fListOfHistos->Add( fhRegionSumPtMaxVsEt );      // At(11)
668
669   fhRegionSumPtMinVsEt = new TH1F("hRegionSumPtMinVsEt",   "P_{T}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
670   fhRegionSumPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
671   fhRegionSumPtMinVsEt->Sumw2();
672   fListOfHistos->Add( fhRegionSumPtMinVsEt );      // At(12)
673   
674   fhRegionMultMax = new TH1I("hRegionMultMax",      "N_{ch}^{90, max}",  21, -0.5,   20.5);
675   fhRegionMultMax->SetXTitle("N_{ch tracks}");
676   fhRegionMultMax->Sumw2();
677   fListOfHistos->Add( fhRegionMultMax );           // At(13)
678
679   fhRegionMultMaxVsEt = new TH1F("hRegionMultMaxVsEt",  "N_{ch}^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
680   fhRegionMultMaxVsEt->SetXTitle("E (GeV hRegionAveSumPtVsEt/c)");
681   fhRegionMultMaxVsEt->Sumw2();
682   fListOfHistos->Add( fhRegionMultMaxVsEt );       // At(14)
683
684   fhRegionMultMinVsEt = new TH1F("hRegionMultMinVsEt",  "N_{ch}^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
685   fhRegionMultMinVsEt->SetXTitle("E (GeV/c)");
686   fhRegionMultMinVsEt->Sumw2();
687   fListOfHistos->Add( fhRegionMultMinVsEt );      // At(15)
688          
689   fhRegionAveSumPtVsEt = new TH1F("hRegionAveSumPtVsEt", "(P_{T}^{90, max} + P_{T}^{90, min})/2 vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
690   fhRegionAveSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
691   fhRegionAveSumPtVsEt->Sumw2();
692   fListOfHistos->Add( fhRegionAveSumPtVsEt );     // At(16)
693
694   fhRegionDiffSumPtVsEt= new TH1F("hRegionPtDiffVsEt", "(P_{T}^{90, max} - P_{T}^{90, min}) vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
695   fhRegionDiffSumPtVsEt->SetXTitle("P_{T} (GeV/c)");
696   fhRegionDiffSumPtVsEt->Sumw2();
697   fListOfHistos->Add( fhRegionDiffSumPtVsEt );    // At(17)
698   
699   fhRegionAvePartPtMaxVsEt = new TH1F("hRegionAvePartPtMaxVsEt", "#LTp_{T}#GT^{90, max} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
700   fhRegionAvePartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
701   fhRegionAvePartPtMaxVsEt->Sumw2();
702   fListOfHistos->Add( fhRegionAvePartPtMaxVsEt );  // At(18)
703
704   fhRegionAvePartPtMinVsEt = new TH1F("hRegionAvePartPtMinVsEt", "#LTp_{T}#GT^{90, min} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
705   fhRegionAvePartPtMinVsEt->SetXTitle("P_{T} (GeV/c)");
706   fhRegionAvePartPtMinVsEt->Sumw2();
707   fListOfHistos->Add( fhRegionAvePartPtMinVsEt );   // At(19)
708
709   fhRegionMaxPartPtMaxVsEt = new TH1F("hRegionMaxPartPtMaxVsEt", "max(p_{T})^{90} vs Leading Jet P_{T}",  fBinsPtInHist, fMinJetPtInHist,   fMaxJetPtInHist);
710   fhRegionMaxPartPtMaxVsEt->SetXTitle("P_{T} (GeV/c)");
711   fhRegionMaxPartPtMaxVsEt->Sumw2();
712   fListOfHistos->Add( fhRegionMaxPartPtMaxVsEt );    // At(20)
713
714 /*   
715   // For debug region selection
716   fhValidRegion = new TH2F("hValidRegion", "dN/d#eta/d#phi",      
717                20, -2.,2., 62, -TMath::Pi(),   TMath::Pi());
718   fhValidRegion->SetYTitle("#Delta#phi");
719   fhValidRegion->Sumw2();
720   fListOfHistos->Add( fhValidRegion );  // At(15)
721 */
722 }
723  
724
725