Compilation with Root6
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetHF.cxx
1 //
2 // Author: A Castro
3
4 #include "AliAnalysisTaskEmcalJetHF.h"
5
6 // general ROOT includes                                                                                                                                                  
7 #include <TCanvas.h>
8 #include <TChain.h>
9 #include <TClonesArray.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TH3F.h>
13 #include <THnSparse.h>
14 #include <TList.h>
15 #include <TLorentzVector.h>
16 #include <TParameter.h>
17 #include <TParticle.h>
18 #include <TTree.h>
19 #include <TVector3.h>
20 #include <TObjArray.h>
21
22 // AliROOT includes                                                                                                                         
23 #include "AliAODEvent.h"
24 #include "AliESDEvent.h"
25 #include "AliAnalysisManager.h"
26 #include "AliAnalysisTask.h"
27 #include "AliCentrality.h"
28 #include "AliEmcalJet.h"
29 #include "AliAODJet.h"
30 #include "AliVCluster.h"
31 #include "AliVTrack.h"
32 #include <AliVEvent.h>
33 #include <AliVParticle.h>
34 #include "AliRhoParameter.h"
35 #include "AliEmcalParticle.h"
36
37 // event handler (and pico's) includes                                                                                                      
38 #include <AliInputEventHandler.h>
39 #include <AliVEventHandler.h>
40 #include "AliESDInputHandler.h"
41 #include "AliPicoTrack.h"
42 #include "AliEventPoolManager.h"
43
44 // PID includes                                                                                                                             
45 #include "AliPIDResponse.h"
46 #include "AliTPCPIDResponse.h"
47 #include "AliESDpid.h"
48
49 #include <AliInputEventHandler.h>
50 #include <AliVEventHandler.h>
51
52 // magnetic field includes
53 #include "TGeoGlobalMagField.h"
54 #include "AliMagF.h"
55
56 using std::cout;
57 using std::endl;
58
59 ClassImp(AliAnalysisTaskEmcalJetHF)
60
61 //________________________________________________________________________
62 AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF() : 
63   AliAnalysisTaskEmcalJet("heavyF",kFALSE), 
64   event(0),
65   fillHist(0),
66   //isESD(0),
67   doGlobalPID(0),
68   fCuts(0),
69   fPhimin(-10), fPhimax(10),
70   fEtamin(-0.9), fEtamax(0.9),
71   fAreacut(0.0),
72   fJetHIpt(50.0),
73   fTrackPtCut(2.0),
74   fTrackEta(0.9),
75   fTrkQAcut(0),
76   fPIDResponse(0x0), fTPCResponse(),
77   fEsdtrackCutsITSTPC(),
78   fEsdtrackCutsTPC(),
79   fEsdtrackCutsITS(),
80   fESD(0), fAOD(0),
81   fHistRhovsCent(0),
82   fHistJetPhi(0),
83   fHistCorJetPt(0), fHistJetPt(0),
84   fHistdEdx(0),
85   fHistdEdxvPt(0),
86   fHistClusE(0),
87   fHistEovPTracks(0),
88   fHistEovPvsPtTracks(0),
89   fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
90   fHistnsigelectron(0),
91   fHistnSigElecPt(0),
92   fHistTrackPhivEta(0),
93   fHistClusterPhivEta(0),
94   fHistnJetTrackvnJetClusters(0),
95   fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
96 {
97   // Default constructor.
98   for (Int_t i = 0;i<6;++i){
99     fHistJetPtvsTrackPt[i]      = 0;
100     fHistTrackPt[i]             = 0;
101     fHistEP0[i]                 = 0;
102     fHistEP0A[i]                = 0;
103     fHistEP0C[i]                = 0;
104     fHistEPAvsC[i]              = 0;
105   }
106
107   SetMakeGeneralHistograms(kTRUE);
108
109 }
110
111 //________________________________________________________________________
112 AliAnalysisTaskEmcalJetHF::AliAnalysisTaskEmcalJetHF(const char *name) :
113   AliAnalysisTaskEmcalJet(name,kTRUE),
114   event(0),
115   fillHist(0),
116   //isESD(0),
117   doGlobalPID(0),
118   fCuts(0),
119   fPhimin(-10), fPhimax(10),
120   fEtamin(-0.9), fEtamax(0.9),
121   fAreacut(0.0),
122   fJetHIpt(50.0),
123   fTrackPtCut(2.0),
124   fTrackEta(0.9),
125   fTrkQAcut(0),
126   fPIDResponse(0x0), fTPCResponse(),
127   fEsdtrackCutsITSTPC(),
128   fEsdtrackCutsTPC(),
129   fEsdtrackCutsITS(),
130   fESD(0), fAOD(0),
131   fHistRhovsCent(0),
132   fHistJetPhi(0),
133   fHistCorJetPt(0), fHistJetPt(0),
134   fHistdEdx(0),
135   fHistdEdxvPt(0),
136   fHistClusE(0),
137   fHistEovPTracks(0),
138   fHistEovPvsPtTracks(0),
139   fHistPID(0), fHistPIDtpc(0), fHistPIDits(0), fHistPIDtof(0),
140   fHistnsigelectron(0),
141   fHistnSigElecPt(0),
142   fHistTrackPhivEta(0),
143   fHistClusterPhivEta(0),
144   fHistnJetTrackvnJetClusters(0),
145   fhnPIDHF(0x0), fhnQA(0x0), fhnJetQA(0x0), fhnClusQA(0x0), fhnTrackQA(0x0), fhnGlobalPID(0x0)
146
147   for (Int_t i = 0;i<6;++i){
148     fHistJetPtvsTrackPt[i]      = 0;
149     fHistTrackPt[i]             = 0;
150     fHistEP0[i]                 = 0;
151     fHistEP0A[i]                = 0;
152     fHistEP0C[i]                = 0;
153     fHistEPAvsC[i]              = 0;
154    }
155    SetMakeGeneralHistograms(kTRUE);
156  
157    DefineInput(0,TChain::Class());
158    DefineOutput(1, TList::Class());
159 }
160
161 //_______________________________________________________________________
162 AliAnalysisTaskEmcalJetHF::~AliAnalysisTaskEmcalJetHF()
163 {
164   // destructor
165   //
166   if (fOutput) {
167     delete fOutput;
168     fOutput = 0;
169   }
170 }
171
172 //________________________________________________________________________
173 void AliAnalysisTaskEmcalJetHF::UserCreateOutputObjects()
174 {
175   if (! fCreateHisto)
176     return;
177   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
178
179   fOutput = new TList();
180   fOutput->SetOwner(kTRUE);
181
182   fHistJetPhi                = new TH1F("NjetvsPhi", "NjetvsPhi", 288,-2*TMath::Pi(),2*TMath::Pi());
183   fHistJetPt                 = new TH1F("NjetvsJetPt", "NjetvsJetPt", 300, 0, 300);
184   fOutput->Add(fHistJetPhi);
185   fOutput->Add(fHistJetPt);
186
187   fillHist = 1;
188
189   fHistClusE = new TH1F("NumberClustersvsEnergy","NumberClustersvsEnergy", 500, 0, 10);
190
191     
192   if(fillHist>0){
193   fHistRhovsCent              = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 400, 0, 400);
194   fHistCorJetPt                     = new TH1F("NjetvsCorrJetPt", "NjetvsCorrJetPt", 300, -100, 200);
195   fHistdEdx                             = new TH1F("dEdxSignal", "dEdxSignal", 500, 0, 500);
196   fHistdEdxvPt                = new TH2F("dEdxvPt", "dEdxvPt", 200, 0, 100,500, 0 ,500);
197   fHistEovPTracks             = new TH1F("EovPTracks","EovPTracks",200,0.0,2.0);
198   fHistEovPvsPtTracks         = new TH2F("E/p_vs_Pt","E/p_vs_Pt",200, 0 ,100 ,50, 0, 2);
199   fHistPID                    = new TH1F("fHistPID", "Detector PID", 8, 0, 8);
200   fHistPIDtpc                 = new TH1F("fHistPIDtpc", "TPC pid", 4, 0, 4);
201   fHistPIDits                 = new TH1F("fHistPIDits", "ITS pid", 4, 0, 4);
202   fHistPIDtof                 = new TH1F("fHistPIDtof", "TOF pid", 4, 0, 4);
203   fHistnsigelectron           = new TH1F("nsig_elec_TPC","nsig_elec_TPC",500,-10,10);
204   fHistnSigElecPt             = new TH2F("nsig_v_pt(TPC)","nsig_v_pt(TPC)",200,0,100,100,-10,10);
205   fHistTrackPhivEta           = new TH2F("TrackPhi_v_Eta","TrackPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
206   fHistClusterPhivEta         = new TH2F("ClusterPhi_v_Eta","ClusterPhiEta",64,-1.6,1.6,72,1,2*TMath::Pi());
207   fHistnJetTrackvnJetClusters = new TH2F("NumbJetTracksvJetClusters","NumbJetTracksvJetClusters",21,0,20,21,0,20);
208
209   // PT bins used to be (2000, -100, 300) 
210   TString name;
211   TString title;
212
213   // creating centrality dependent histos that don't involve Global Rho
214   for (Int_t i = 0;i<6;++i){
215     name = TString(Form("JetPtvsTrackPt_%i",i));
216     title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
217     fHistJetPtvsTrackPt[i] = new TH2F(name,title, 500, -100, 400, 100,0,100);
218     fOutput->Add(fHistJetPtvsTrackPt[i]);
219
220     name = TString(Form("TrackPt_%i",i));
221     title = TString(Form("Track pT cent bin %i",i));
222     fHistTrackPt[i] = new TH1F(name,title,400,0,200);
223     fOutput->Add(fHistTrackPt[i]);   
224
225     name = TString(Form("EP0_%i",i));
226     title = TString(Form("EP VZero cent bin %i",i));
227     fHistEP0[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
228     fOutput->Add(fHistEP0[i]);
229
230     name = TString(Form("EP0A_%i",i));
231     title = TString(Form("EP VZero cent bin %i",i));
232     fHistEP0A[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
233     fOutput->Add(fHistEP0A[i]);
234
235     name = TString(Form("EP0C_%i",i));
236     title = TString(Form("EP VZero cent bin %i",i));
237     fHistEP0C[i] = new TH1F(name,title,144,-TMath::Pi(),TMath::Pi());
238     fOutput->Add(fHistEP0C[i]);
239
240     name = TString(Form("EPAvsC_%i",i));
241     title = TString(Form("EP VZero cent bin %i",i));
242     fHistEPAvsC[i] = new TH2F(name,title,144,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
243     fOutput->Add(fHistEPAvsC[i]);
244
245   }
246
247   fOutput->Add(fHistRhovsCent);
248   fOutput->Add(fHistCorJetPt);
249   fOutput->Add(fHistdEdx);
250   fOutput->Add(fHistdEdxvPt);
251   fOutput->Add(fHistEovPTracks);
252   fOutput->Add(fHistEovPvsPtTracks);
253   fOutput->Add(fHistPID);
254   fOutput->Add(fHistPIDtpc);
255   fOutput->Add(fHistPIDits);
256   fOutput->Add(fHistPIDtof);
257   fOutput->Add(fHistnsigelectron);
258   fOutput->Add(fHistnSigElecPt);
259   fOutput->Add(fHistTrackPhivEta);
260   fOutput->Add(fHistClusterPhivEta);
261   fOutput->Add(fHistnJetTrackvnJetClusters);
262   
263   }//Fill Histograms
264
265     fOutput->Add(fHistClusE);
266
267
268   // ****************************** PID *****************************************************                                               
269   // set up PID handler                                                                                                                     
270   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
271   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
272   if(!inputHandler) {
273     AliFatal("Input handler needed");
274   }
275
276   // PID response object                                                                                                                    
277   //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();                                                                         
278   //  inputHandler->CreatePIDResponse(fIsMC);         // needed to create object, why though?                                                 
279   fPIDResponse = inputHandler->GetPIDResponse();
280   if (!fPIDResponse) {
281     AliError("PIDResponse object was not created");
282   }
283   // ****************************************************************************************
284   UInt_t bitcoded = 0;  // bit coded, see GetDimParamsPID() below
285   bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
286   fhnPIDHF = NewTHnSparseDHF("fhnPIDHF", bitcoded);
287   
288   UInt_t bitcoded1 = 0;  // bit coded, see GetDimParamsPID() below
289   bitcoded1 = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4;
290   fhnJetQA = NewTHnSparseDJetQA("fhnJetQA", bitcoded1);
291   
292   UInt_t bitcoded2 = 0;  // bit coded, see GetDimParamsPID() below
293   bitcoded2 = 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15;
294   fhnClusQA = NewTHnSparseDJetQA("fhnClusQA", bitcoded2);
295   
296   UInt_t bitcoded3 = 0;  // bit coded, see GetDimParamsPID() below
297   bitcoded3 = 1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14;
298   fhnTrackQA = NewTHnSparseDJetQA("fhnTrackQA", bitcoded3);
299   
300   cout << "_______________Created Sparse__________________" << endl;
301   
302   fOutput->Add(fhnPIDHF);
303   fOutput->Add(fhnJetQA);
304   fOutput->Add(fhnClusQA);
305   fOutput->Add(fhnTrackQA);
306   fOutput->Add(fhnGlobalPID);
307
308
309   PostData(1, fOutput);
310 }
311
312 //________________________________________________________
313 void AliAnalysisTaskEmcalJetHF::ExecOnce()
314 {
315   //  Initialize the analysis 
316   AliAnalysisTaskEmcalJet::ExecOnce();
317
318 } // end of ExecOnce
319
320 //________________________________________________________________________
321 Bool_t AliAnalysisTaskEmcalJetHF::Run()
322 {
323   // check to see if we have any tracks
324   if (!fTracks)  return kTRUE;
325   if (!fJets)  return kTRUE;
326   
327   // what kind of event do we have: AOD or ESD?
328   Bool_t useAOD;
329   if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
330   else useAOD = kFALSE;
331   
332   // if we have ESD event, set up ESD object
333   if(!useAOD){
334     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
335     if (!fESD) {
336       AliError(Form("ERROR: fESD not available\n"));
337       return kTRUE;
338     }
339   }
340   
341   // if we have AOD event, set up AOD object
342   if(useAOD){
343     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
344     if(!fAOD) {
345       AliError(Form("ERROR: fAOD not available\n"));
346       return kTRUE;
347     }
348   }
349   
350   // get magnetic field info for DCA
351   Double_t  MagF = fESD->GetMagneticField();
352   Double_t MagSign = 1.0;
353   if(MagF<0)MagSign = -1.0;
354   // set magnetic field
355   if (!TGeoGlobalMagField::Instance()->GetField()) {
356     AliMagF* field = new AliMagF("Maps","Maps", MagSign, MagSign, AliMagF::k5kG);
357     TGeoGlobalMagField::Instance()->SetField(field);
358   }
359
360   // get centrality bin
361   Int_t centbin = GetCentBin(fCent);
362   //for pp analyses we will just use the first centrality bin
363   if (centbin == -1)  centbin = 0;
364
365   // get vertex information
366   Double_t fvertex[3]={0,0,0};
367   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
368   //Double_t zVtx=fvertex[2];
369
370   // create pointer to list of input event                                                                                                  
371   TList *list = InputEvent()->GetList();
372   if(!list) {
373     AliError(Form("ERROR: list not attached\n"));
374     return kTRUE;
375   }
376
377   // background density                                                                                                                                                                                                                               
378   fRhoVal = fRho->GetVal();
379
380   // initialize TClonesArray pointers to jets and tracks                                                                                    
381   TClonesArray *jets = 0;
382   TClonesArray *tracks = 0;
383   TClonesArray *clusters = 0;
384     
385   // get Jets object                                                                                                                        
386   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
387   if(!jets){
388     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
389     return kTRUE;
390   } // verify existence of jets                                                                                                             
391
392   // get Tracks object                                                                                                                      
393   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
394   if (!tracks) {
395     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
396     return kTRUE;
397   } // verify existence of tracks  (fTracksName.Data())
398     
399   //get Clusters object
400   clusters = dynamic_cast<TClonesArray*>(list->FindObject(fCaloClusters));
401   if (!clusters){
402     AliError(Form("Pointer to tracks %s == 0",fCaloClusters->GetName()));
403     return kTRUE;
404   }  //verify cluster existence
405   
406   
407   //const Int_t nclus = clusters->GetEntries();
408   /* Globalcluster
409   for(int icluster = 0; icluster < nclus; icluster++) {
410     AliVCluster* andyCluster = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
411     if (!icluster) continue;
412   }
413   */
414   // get all emcal clusters
415   TRefArray* caloClusters = new TRefArray();
416   fESD->GetEMCALClusters( caloClusters );
417         
418   //TObjArray* listcuts = fEsdtrackCutsTPC->GetAcceptedTracks(fESD);
419   //Int_t nGoodTracks = list->GetEntries();
420   Int_t nCluster = caloClusters->GetEntries();
421
422   // event plane histograms filled
423   if(fillHist>0) fHistEP0[centbin]->Fill(fEPV0);
424   if(fillHist>0) fHistEP0A[centbin]->Fill(fEPV0A);
425   if(fillHist>0) fHistEP0C[centbin]->Fill(fEPV0C);
426   if(fillHist>0) fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
427   
428   event++;
429   cout<<"Event #: "<<event<<endl;
430   
431   // get number of jets and tracks                                                                                                          
432   const Int_t Njets = jets->GetEntries();
433   const Int_t Ntracks = tracks->GetEntries();
434   //const Int_t Nclusters = clusters->GetEntries();
435   if(Ntracks<1)   return kTRUE;
436   if(Njets<1)     return kTRUE;
437   //if(Nclusters<1)  return kTRUE;
438   
439   if(doGlobalPID){
440   // loop over Globap tracks in the event for PID: Runs ON ESD's
441   for (int i = 0;i<Ntracks;i++){
442     //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
443     AliVTrack *trackGlobal = static_cast<AliVTrack*>(fTracks->At(i));
444     //AliESDtrack* track = fESD->GetTrack(i);
445     if (! trackGlobal) {
446       AliError(Form("Couldn't get ESD track %d\n", i));
447       continue;
448     }
449     
450
451     // apply track cuts                                                                                                        
452     if(TMath::Abs(trackGlobal->Eta())>fTrackEta) continue;
453     if (trackGlobal->Pt()<0.15) continue;
454
455     // initialize track info 
456     Double_t dEdxGlobal = -99;
457     Double_t trackphiGlobal = -99;
458     Double_t trackptGlobal = 0;
459     Double_t p = trackGlobal->P();
460     Double_t fClsEGlobal = -99;
461     Double_t EovPGlobal = -99;
462
463     // distance of closest approach
464     Float_t DCAxyGlobal = -999;
465     Float_t DCAzGlobal = -999;
466     //Double_t DCAXYGlobal = -999;
467     //Double_t DCAZGlobal = -999;
468     
469     // track info    
470     trackGlobal->GetTPCsignal();
471     //trackGlobal->GetImpactParameters(DCAxyGlobal, DCAzGlobal);
472     dEdxGlobal = trackGlobal->GetTPCsignal();
473     trackptGlobal = trackGlobal->Pt();
474     trackphiGlobal = trackGlobal->Phi();
475
476     // fill track histo's
477     //if(fillHist>0) fHistdEdx->Fill(dEdx);
478     //if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
479     //if(fillHist>0) fHistTrackPt[centbin]->Fill(track->Pt());
480
481     // clusters                                                                                                        
482     Int_t nMatchClusGlobal = -1;
483     AliESDCaloCluster *matchedClusGlobal = NULL;
484                                                                                                                                    
485     //////////////////////////////////////////////////////////////////////////                                                                            
486  
487     // cut on 1.5 GeV for EMCal Cluster                                                                                                                        //if(track->GetEMCALcluster()<0 || pt<1.5) continue;                                                                                            
488     //////////////////////////////////////////////////////////////////////////                                                                                   
489     // particles in TOF                                                                                                                                    
490     //Double_t nSigmaPion_TOFGlobal, nSigmaProton_TOFGlobal, nSigmaKaon_TOFGlobal = -1.;
491
492     // misc quantities                                                                                                                                     
493     //Double_t TOFsigGlobal  = -1.;
494     //Double_t nClustersTPCGlobal = -1;
495     //Int_t chargeGlobal     = 0;
496     //Int_t trackCutsGlobal  = 0;     // set to 0 to get to run
497     //Int_t myPIDGlboal      = 0;     // set to 0 because myPID is to compare with MC unless I hard code in cuts
498     Double_t nSigmaElectron_TPCGlobal = -999;
499     Double_t nSigmaElectron_TOFGlobal = -999;
500     //Double_t nSigmaElectron_EMCALGlobal = -999;
501
502     // get clusters                                                                                                                                       
503     //Int_t clusiGlobal = trackGlobal->GetEMCALcluster();
504     //nClustersTPCGlobal = trackGlobal->GetTPCclusters(0);
505     
506     AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
507  
508     nSigmaElectron_TPCGlobal = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
509     nSigmaElectron_TOFGlobal = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
510     
511     //for EMCAL
512     nMatchClusGlobal = trackGlobal->GetEMCALcluster();
513     if(nMatchClusGlobal > 0){
514       matchedClusGlobal = (AliESDCaloCluster*)fESD->GetCaloCluster(nMatchClusGlobal);
515       //AliESDCaloCluster* matchedClus = fESD->GetCaloCluster(clusi);                                                                                                                                                                                    
516       //double eop2Global = -1;
517       //double ssGl[4]={0.,0.,0.,0.};
518       //Double_t nSigmaEop = fPID->GetPIDResponse()-m >NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss);
519
520       fClsEGlobal       = matchedClusGlobal->E();
521       EovPGlobal        = fClsEGlobal/p;
522       //nSigmaElectron_EMCALGlobal = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
523       //if(fillHist>0) fHistEovPTracks->Fill(EovP);
524       //if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
525     }
526
527     //if(fillHist>0) fHistnsigelectron->Fill(nSigmaElectron_TPC);
528     //if(fillHist>0) fHistnSigElecPt->Fill(trackpt,nSigmaElectron_TPC);
529
530     // extra attempt                                                                                                                                      
531     AliVEvent *eventQA=InputEvent();
532     if (!eventQA||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
533     //////////////////////////////////////////////////////////////////////////                                                                             
534         // cut on 1.5 GeV for EMCal Cluster
535     if(trackGlobal->Pt() < fTrackPtCut) continue;
536
537     Double_t HF_tracks[12] = {fCent, trackGlobal->Pt(), trackGlobal->P() ,trackGlobal->Eta(), trackGlobal->Phi(), EovPGlobal, DCAxyGlobal, DCAzGlobal, dEdxGlobal, nSigmaElectron_TPCGlobal, nSigmaElectron_TOFGlobal, 0};//nSigmaElectron_EMCAL};
538     fhnGlobalPID->Fill(HF_tracks);    // fill Sparse Histo with trigger entries
539
540   } // Loop over tracks
541   }//PID Switch
542   
543     //  Start Jet Analysis
544     // initialize jet parameters
545     Int_t ijethi=-1;
546     Double_t highestjetpt=0.0;
547     
548     // loop over jets in an event - to find highest jet pT and apply some cuts && JetQA Sparse
549     for (Int_t ijet = 0; ijet < Njets; ijet++){
550         // get our jets
551         AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
552         if (!jet) continue;
553       
554         // apply jet cuts
555         if(!AcceptMyJet(jet)) continue;
556  
557         Double_t JetQA[5] = {0, static_cast<Double_t>(jet->GetNumberOfTracks()), static_cast<Double_t>(jet->GetNumberOfClusters()),jet->Eta(), jet->Phi()};
558         fhnJetQA->Fill(JetQA);
559       
560       
561         // Loop over clusters for JetQA
562         for(int iCluster = 0; iCluster < nCluster; iCluster++) {
563           AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
564           //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(jet->GetNumberOfClusters());
565           //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
566           //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
567           if (!caloCluster){
568             AliError(Form("ERROR: Could not get cluster %d", iCluster));
569             continue;
570           }
571           if(!IsJetCluster(jet,iCluster,kFALSE)) continue ;
572           Float_t pos[3];
573           caloCluster->GetPosition(pos);  // Get cluster position
574           TVector3 cp(pos);
575           Double_t NtrMatched = -999.0;
576           NtrMatched = caloCluster->GetNTracksMatched();
577           
578           //loop over tracks for Jet QA
579           for(int itrack = 0; itrack < NtrMatched; itrack++){
580             //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
581             AliVTrack *trackcluster = static_cast<AliVTrack*>(fTracks->At(itrack));
582             //AliESDtrack* track = fESD->GetTrack(i);
583             if (! trackcluster) {
584               AliError(Form("Couldn't get ESD track %d\n", itrack));
585               continue;
586             }
587             if(!IsJetTrack(jet,itrack,kFALSE)) continue;
588             Double_t ClusQA[6] = {caloCluster->E(),cp.PseudoRapidity() ,cp.Phi(), static_cast<Double_t>(caloCluster->IsEMCAL()), NtrMatched, trackcluster->Phi()};
589             fhnClusQA->Fill(ClusQA); //,1./Njets);    // fill Sparse Histo with trigger entries
590           }//loop over tracks for JetQA
591
592         }//loop over clusters for JetQA
593       
594       // loop over tracks in the event
595       for (int iTrack = 0; iTrack<jet->GetNumberOfTracks(); iTrack++){  //loop over tracks in jet for TrackQA
596         //AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
597         AliVTrack *jetTrack = static_cast<AliVTrack*>(fTracks->At(iTrack));
598         //AliESDtrack* track = fESD->GetTrack(i);
599         if (! jetTrack) {
600           AliError(Form("Couldn't get ESD track %d\n", iTrack));
601           continue;
602         }
603         if(!IsJetTrack(jet,iTrack,kFALSE)) continue;
604         Double_t trackQA[5] = {jetTrack->Pt(), jetTrack->Eta(), jetTrack->Phi(), static_cast<Double_t>(jetTrack->IsEMCAL()), static_cast<Double_t>(jetTrack->GetEMCALcluster())};
605         fhnTrackQA->Fill(trackQA); //,1./Njets);
606         
607       }//track loop for TrackQA
608
609         if(highestjetpt<jet->Pt()){
610             ijethi=ijet;
611             highestjetpt=jet->Pt();
612         }
613     } // end of looping over jets
614
615   // loop over jets in the event and make appropriate cuts
616   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
617      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
618      if (!jet)  // see if we have a jet
619        continue;
620       
621      // phi of jet, constrained to 1.6 < Phi < 2.94
622      float jetphi = jet->Phi();      // phi of jet 
623   
624     // apply jet cuts
625     if(!AcceptMyJet(jet)) continue;
626  
627     Int_t JetClusters = jet->GetNumberOfClusters();
628     Int_t JetTracks = jet -> GetNumberOfTracks();
629     fHistnJetTrackvnJetClusters->Fill(JetClusters,JetTracks);
630    
631      // Initializations and Calculations
632      //Double_t jeteta = jet->Eta();                                    // ETA of jet
633      Double_t jetptraw = jet->Pt();                             // raw pT of jet
634      Double_t jetPt = -500;                                     // initialize corr jet pt LOCAL
635      Double_t jetarea = -500;                                   // initialize jet area
636      jetarea = jet->Area();                                     // jet area
637      jetPt = jet->Pt() - jetarea*fRhoVal;                 // semi-corrected pT of jet from GLOBAL rho value
638
639      // make histo's
640      if(fillHist>0) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
641      fHistJetPhi->Fill(jetphi);
642      if(fillHist>0) fHistCorJetPt->Fill(jetPt);
643      fHistJetPt->Fill(jetptraw);
644     
645       if(jet->Pt() > fJetHIpt) {
646         
647       //loop over clusters
648       //for (int i = 0; i < njetclusters; i++){
649       for(int iCluster = 0; iCluster < nCluster; iCluster++) {
650           AliVCluster* caloCluster = static_cast<AliVCluster*>(fCaloClusters->At(iCluster));
651           //AliVCluster* caloCluster = (AliVCluster* )caloClusters->At(iCluster);
652           //AliESDCaloCluster* caloCluster = (AliESDCaloCluster* )caloClusters->At(iCluster);
653           //AliESDCaloCluster* clus = fESD->GetCaloCluster(iclus);
654           if (!caloCluster){
655               AliError(Form("ERROR: Could not get cluster %d", iCluster));
656               continue;
657           }
658         
659           //if (!caloCluster -> IsEMCAL()) continue; //Check that Cluster is EMCal Cluster
660           if(! IsJetCluster(jet, iCluster, kFALSE)) continue;
661         
662           AliESDtrack *track = 0;
663           if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
664           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
665           Double_t fclusE = -999;
666           fclusE = caloCluster->E();
667         
668           if (fclusE<0.) continue;  //Check that cluster has positive energy
669         
670           fHistClusE->Fill(fclusE);
671           //Int_t  cluslabel = caloCluster->GetID();
672           //Int_t nmatched = caloCluster->GetNTracksMatched();
673         
674           Float_t pos[3];
675           caloCluster->GetPosition(pos);  // Get cluster position
676           TVector3 cp(pos);
677           if(fillHist>0) fHistClusterPhivEta->Fill(cp.PseudoRapidity(),cp.Phi());
678         
679           Int_t trackMatchedindex=caloCluster->GetTrackMatchedIndex();
680           if(trackMatchedindex<0)continue;  // Make sure we don't have a bad index
681  
682           if (caloCluster->GetTrackMatchedIndex() > 0) // tender's matching
683           track = fESD->GetTrack(caloCluster->GetTrackMatchedIndex());
684           Double_t NtrMatched = -999;
685           NtrMatched = caloCluster->GetNTracksMatched();
686           for(int itrack = 0; itrack < NtrMatched; itrack++){   // Loop over tracks matched to clusters from jets
687             //AliVParticle *trackCluster = static_cast<AliVParticle*>(fTracks->At(i));
688             AliVTrack *trackCluster = static_cast<AliVTrack*>(fTracks->At(itrack));
689             //AliESDtrack* trackCluster = fESD->GetTrack(itrack);
690             if (! trackCluster) {
691               AliError(Form("Couldn't get ESD track %d\n", itrack));
692               continue;
693             }
694             if(!IsJetTrack(jet,itrack,kFALSE)) continue;  // Check that track is still part of ith jet
695             
696               //if (trackCluster->Phi()>3.2 || trackCluster->Phi()<1.4)cout <<"Out of Range Track!    Track Phi:  " << trackCluster->Phi() << endl;
697   
698               //if(track->Pt() < 4.0) continue;
699               // initialize track info
700               Double_t dEdx = -99;
701               //Double_t trackphi = -99;
702               Double_t trackpt = 0;
703               Double_t p = trackCluster->P();
704               Double_t EovP = -99;
705               
706               // distance of closest approach
707               Float_t DCAxy = -999;
708               Float_t DCAz = -999;
709               //Double_t DCAXY = -999;
710               //Double_t DCAZ = -999;
711               
712               // track info
713               //trackCluster->GetImpactParameters(DCAxy, DCAz);
714               trackpt = trackCluster->Pt();
715               //trackphi = track->Phi();
716               
717               // fill track histo's
718               if(fillHist>0) fHistTrackPhivEta->Fill(trackCluster->Eta(),trackCluster->Phi());
719               if(fillHist>0) fHistdEdx->Fill(dEdx);
720               if(fillHist>0) fHistdEdxvPt->Fill(trackpt,dEdx);
721               if(fillHist>0) fHistTrackPt[centbin]->Fill(trackCluster->Pt());
722             
723               Double_t nSigmaElectron_TPC = -999;
724               Double_t nSigmaElectron_TOF = -999;
725               Double_t nSigmaElectron_EMCAL = -999;
726   
727             if(!useAOD){
728               
729               AliESDtrack *trackESD = fESD->GetTrack(Ntracks);
730               
731               dEdx = trackESD->GetTPCsignal();
732               trackESD->GetImpactParameters(DCAxy, DCAz);
733               nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kElectron);
734               nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kElectron);
735             }
736            if (useAOD) {
737               AliAODTrack *trackAOD = fAOD->GetTrack(Ntracks);
738               
739               // get detector signals
740               dEdx = trackAOD->GetTPCsignal();
741              
742              //trackAOD->GetImpactParameters(DCAxy,DCAz);
743               nSigmaElectron_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kElectron);
744               nSigmaElectron_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kElectron);
745               
746             } // end of AOD pid
747       
748               //double eop2 = -1;
749               //double ss[4]={0.,0.,0.,0.};
750               EovP        = fclusE/p;
751               //nSigmaElectron_EMCAL = fPIDResponse->NumberOfSigmasEMCAL(trackESD,AliPID::kElectron,eop2,ss);
752               if(fillHist>0) fHistEovPTracks->Fill(EovP);
753               if(fillHist>0) fHistEovPvsPtTracks->Fill(trackpt,EovP);
754               Double_t HF_tracks[12] = {fCent, trackCluster->Pt(), trackCluster->P() ,trackCluster->Eta(), trackCluster->Phi(), EovP, 0/*DCAxy*/, 0/*DCAz*/, dEdx, nSigmaElectron_TPC, nSigmaElectron_TOF, nSigmaElectron_EMCAL};
755               fhnPIDHF->Fill(HF_tracks);    // fill Sparse Histo with trigger entries
756           
757           //Int_t trackMatchedIndex = caloCluster->GetTrackMatchedIndex();//find the index of the matched track. This by default returns the best match
758           //AliESDtrack *track = event->GetTrack(trackMatchedIndex);
759           //if this is a good track, accept track will return true. The track matched is good, so not track matched is false
760           //Int_t matched = fESD->AcceptTrack(track);//If the track is bad, don't count it.  By default even bad tracks are accepted
761
762         }  //loop over tracks matched to clusters in jet
763         
764       } // loop over jet clusters
765     
766           
767       } // highest pt jet cut
768   } // LOOP over JETS in event
769
770   return kTRUE;
771    
772 }
773 //________________________________________________________________________
774 void AliAnalysisTaskEmcalJetHF::Terminate(Option_t *)
775 {
776   cout<<"###########################"<<endl;
777   cout<<"####   Task Finished   ####"<<endl;
778   cout<<"###########################"<<endl;
779   cout<<"###########################"<<endl;
780 } // end of terminate
781
782 //________________________________________________________________________
783 Int_t AliAnalysisTaskEmcalJetHF::AcceptMyJet(AliEmcalJet *jet) {
784   //applies all jet cuts except pt
785   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
786   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
787   if (jet->Area()<fAreacut) return 0;
788   // prevents 0 area jets from sneaking by when area cut == 0
789   if (jet->Area()==0) return 0;
790   //exclude jets with extremely high pt tracks which are likely misreconstructed
791   if(jet->MaxTrackPt()>100) return 0;
792   
793   //passed all above cuts
794   return 1;
795 }
796
797 //________________________________________________________________________
798 Int_t AliAnalysisTaskEmcalJetHF::GetCentBin(Double_t cent) const
799 {  // Get centrality bin.
800
801   Int_t centbin = -1;
802   if (cent>=0 && cent<10)
803     centbin = 0; 
804   else if (cent>=10 && cent<20)
805     centbin = 1;
806   else if (cent>=20 && cent<30)
807     centbin = 2;
808   else if (cent>=30 && cent<40)
809     centbin = 3;
810   else if (cent>=40 && cent<50)
811     centbin = 4;
812   else if (cent>=50 && cent<90)
813     centbin = 5;
814   return centbin;
815
816 //____________________________________________________________________________________________                                                                            
817 THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDHF(const char* name, UInt_t entries)
818 {
819   // generate new THnSparseD PID, axes are defined in GetDimParams()                                                                                                     
820   Int_t count = 0;
821   UInt_t tmp = entries;
822   while(tmp!=0){
823     count++;
824     tmp = tmp &~ -tmp;  // clear lowest bit                                                                                                                             
825   }
826
827   TString hnTitle(name);
828   const Int_t dim = count;
829   Int_t nbins[dim];
830   Double_t xmin[dim];
831   Double_t xmax[dim];
832
833   Int_t i=0;
834   Int_t c=0;
835   while(c<dim && i<32){
836     if(entries&(1<<i)){
837
838       TString label("");
839       GetDimParamsHF(i, label, nbins[c], xmin[c], xmax[c]);
840       hnTitle += Form(";%s",label.Data());
841       c++;
842     }
843
844     i++;
845   }
846   hnTitle += ";";
847
848   return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
849 } // end of NewTHnSparseF PID
850
851 THnSparse* AliAnalysisTaskEmcalJetHF::NewTHnSparseDJetQA(const char* name, UInt_t entries)
852 {
853   // generate new THnSparseD JetQA, axes are defined in GetDimParamsJetQA()
854   Int_t count = 0;
855   UInt_t tmp = entries;
856   while(tmp!=0){
857     count++;
858     tmp = tmp &~ -tmp;  // clear lowest bit
859   }
860   
861   TString hnTitle(name);
862   const Int_t dim = count;
863   Int_t nbins[dim];
864   Double_t xmin[dim];
865   Double_t xmax[dim];
866   
867   Int_t i=0;
868   Int_t c=0;
869   while(c<dim && i<32){
870     if(entries&(1<<i)){
871       
872       TString label("");
873       GetDimParamsJetQA(i, label, nbins[c], xmin[c], xmax[c]);
874       hnTitle += Form(";%s",label.Data());
875       c++;
876     }
877     
878     i++;
879   }
880   hnTitle += ";";
881   
882   return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
883 } // end of NewTHnSparseF JetQA
884
885
886 void AliAnalysisTaskEmcalJetHF::GetDimParamsHF(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
887 {
888   // stores label and binning of axis for THnSparse                                                                                                                      
889   const Double_t pi = TMath::Pi();
890
891   switch(iEntry){
892
893   case 0:
894     label = "V0 centrality (%)";
895     nbins = 10;
896     xmin = 0.;
897     xmax = 100.;
898     break;
899
900   case 1:
901     label = "Track p_{T}";
902     nbins = 750;
903     xmin = 0.;
904     xmax = 75.;
905     break;
906
907   case 2:
908     label = "Track p";
909     nbins = 750;
910     xmin = 0.;
911     xmax = 75.;
912     break;
913
914   case 3:
915     label = "Track Eta";
916     nbins = 64;
917     xmin = -1.6;
918     xmax = 1.6;
919     break;
920
921   case 4:
922     label = "Track Phi";
923     nbins = 72;
924     xmin = 0;
925     xmax = 2*pi;
926     break;
927
928   case 5:
929     label = "E/p of track";
930     nbins = 400;
931     xmin = 0;
932     xmax = 4.0;
933     break;
934
935  case 6:
936     label = "DCA xy";
937     nbins = 20;
938     xmin = -10;
939     xmax =  10;
940     break;
941
942   case 7:
943     label = "DCA z";
944     nbins = 20;
945     xmin = -10;
946     xmax = 10;
947     break;
948
949   case 8:                                                                                                                                               
950     label = "dEdX of track - TPC";
951     nbins = 300;
952     xmin = 0;
953     xmax = 300;
954     break;
955
956   case 9:                                                                                                                                                
957     label = "nSigma electron TPC";
958     nbins = 50;
959     xmin = -5;
960     xmax = 5;
961     break;
962
963    case 10:
964     label = "nSigma electron TOF";
965     nbins = 50;
966     xmin = -5;
967     xmax = 5;
968     break;
969
970    case 11:
971     label = "nSigma electron Emcal";
972     nbins = 50;
973     xmin = -5;
974     xmax = 5;
975     break;
976
977
978   } // end of switch                                                                                                                                                     
979 } // end of getting dim-params
980
981 void AliAnalysisTaskEmcalJetHF::GetDimParamsJetQA(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
982 {
983   // stores label and binning of axis for THnSparse
984   const Double_t pi = TMath::Pi();
985   
986   switch(iEntry){
987       
988     case 0:
989       label = "number of Jets in Event";
990       nbins = 2;
991       xmin = 0.;
992       xmax = 1.;
993       break;
994       
995     case 1:
996       label = "number of Clusters in a Jet";
997       nbins = 100;
998       xmin = 0.;
999       xmax = 100.;
1000       break;
1001       
1002     case 2:
1003       label = "number of Tracks in a Jet";
1004       nbins = 100;
1005       xmin = 0.;
1006       xmax = 100.;
1007       break;
1008       
1009     case 3:
1010       label = "Jet Eta";
1011       nbins = 24;
1012       xmin = -1.2;
1013       xmax = 1.2;
1014       break;
1015       
1016     case 4:
1017       label = "Jet Phi";
1018       nbins = 72;
1019       xmin = 0;
1020       xmax = 2*pi;
1021       break;
1022       
1023     case 5:
1024       label = "Cluster Energy";
1025       nbins = 1000;
1026       xmin = 0;
1027       xmax = 10;
1028       break;
1029       
1030     case 6:
1031       label = "Cluster Eta";
1032       nbins = 24;
1033       xmin = -1.2;
1034       xmax =  1.2;
1035       break;
1036       
1037     case 7:
1038       label = "Cluster Phi";
1039       nbins = 72;
1040       xmin = 0;
1041       xmax = 2*pi;
1042       break;
1043       
1044     case 8:
1045       label = "Is EMCalCluster";
1046       nbins = 2;
1047       xmin = 0;
1048       xmax = 2;
1049       break;
1050       
1051     case 9:
1052       label = "Number of Matched Tracks to Cluster";
1053       nbins = 60;
1054       xmin = 0;
1055       xmax = 60;
1056       break;
1057       
1058     case 10:
1059       label = "Track Pt";
1060       nbins = 750;
1061       xmin = 0;
1062       xmax = 75;
1063       break;
1064       
1065     case 11:
1066       label = "Track Eta";
1067       nbins = 24;
1068       xmin = -1.2;
1069       xmax = 1.2;
1070       break;
1071       
1072     case 12:
1073       label= "Track Phi";
1074       nbins = 72;
1075       xmin = 0;
1076       xmax = 2*pi;
1077       break;
1078       
1079     case 13:
1080       label="Is Track EMCal";
1081       nbins = 2;
1082       xmin = 0;
1083       xmax = 2;
1084       break;
1085       
1086     case 14:
1087       label = "Get Track EMCal Cluster";
1088       nbins = 100;
1089       xmin = 0;
1090       xmax = 100;
1091       break;
1092       
1093     case 15:
1094       label = "Track Matched Phi";
1095       nbins = 72;
1096       xmin = 0;
1097       xmax = 2*pi;
1098       
1099       
1100       
1101       
1102   } // end of switch
1103 } // end of getting dim-params
1104
1105
1106