Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliAnalysisTaskSEHFCJqa.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 /////////////////////////////////////////////////////////////
17 //
18 // Class AliAnalysisTaskSEHFCJqa
19 // AliAnalysisTaskSE for the extraction of the fraction of prompt charm
20 // using the charm hadron impact parameter to the primary vertex
21 //
22 // Author: Andrea Rossi, andrea.rossi@pd.infn.it
23 /////////////////////////////////////////////////////////////
24
25
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TH3F.h>
29 #include <TAxis.h>
30 #include <THnSparse.h>
31 #include <TDatabasePDG.h>
32 #include <TMath.h>
33 #include <TROOT.h>
34 #include "AliAODEvent.h"
35 #include "AliAODRecoDecayHF2Prong.h"
36 #include "AliAODRecoDecayHF.h"
37 #include "AliAODRecoDecay.h"
38 #include "AliAnalysisDataSlot.h"
39 #include "AliAnalysisDataContainer.h"
40 #include "AliAODTrack.h"
41 #include "AliAODHandler.h"
42 #include "AliESDtrack.h"
43 #include "AliAODVertex.h"
44 #include "AliESDVertex.h"
45 #include "AliVertexerTracks.h"
46 #include "AliAODMCParticle.h"
47 #include "AliAODPid.h"
48 #include "AliTPCPIDResponse.h"
49 #include "AliAODMCHeader.h"
50 #include "AliAnalysisVertexingHF.h"
51 #include "AliAnalysisTaskSEHFCJqa.h"
52 #include "AliRDHFCutsD0toKpi.h"
53 #include "AliAODInputHandler.h"
54 #include "AliAnalysisManager.h"
55 #include "AliNormalizationCounter.h"
56
57 class TCanvas;
58 class TTree;
59 class TChain;
60 class AliAnalysisTaskSE;
61
62
63 ClassImp(AliAnalysisTaskSEHFCJqa)
64
65   AliAnalysisTaskSEHFCJqa::AliAnalysisTaskSEHFCJqa() 
66 : AliAnalysisTaskSE(),
67   fReadMC(),
68   ffilterbit(),
69   fKeepTrackNegID(),
70   fpidResp(),
71   fCuts(),
72   fhEventCounter(),
73   fhImpParResolITSsel(),
74   fhImpParResolITSselGoodTracks(),
75   fhSparseFilterMask(),
76   fhSparseFilterMaskTrackAcc(),
77   fhSparseFilterMaskImpPar(),
78   fhSparseEoverPeleTPC(),
79   fhSparseShowShapeEleTPC(),
80   fhnSigmaTPCTOFEle(),
81   fhnSigmaTPCTOFPion(),
82   fhnSigmaTPCTOFKaon(),
83   fhnSigmaTPCTOFProton(),
84   fhTrackEMCal(),
85   fSparseRecoJets(),
86   fLoadJet(),
87   fJetArrayString(),
88   fListTrackAndPID(),
89   fListJets()
90 {// standard constructor
91   
92 }
93
94
95
96 //________________________________________________________________________
97 AliAnalysisTaskSEHFCJqa::AliAnalysisTaskSEHFCJqa(const char *name) 
98   : AliAnalysisTaskSE(name),
99     fReadMC(),
100     ffilterbit(AliAODTrack::kTrkGlobalNoDCA),
101     fKeepTrackNegID(),
102     fpidResp(),
103     fCuts(),
104     fhEventCounter(),
105     fhImpParResolITSsel(),
106     fhImpParResolITSselGoodTracks(),
107     fhSparseFilterMask(),
108     fhSparseFilterMaskTrackAcc(),
109     fhSparseFilterMaskImpPar(),
110     fhSparseEoverPeleTPC(),
111     fhSparseShowShapeEleTPC(),
112     fhnSigmaTPCTOFEle(),
113     fhnSigmaTPCTOFPion(),
114     fhnSigmaTPCTOFKaon(),
115     fhnSigmaTPCTOFProton(),
116     fhTrackEMCal(),
117     fSparseRecoJets(),
118     fLoadJet(),
119     fJetArrayString(),
120     fListTrackAndPID(),
121     fListJets()
122 { // default constructor
123   
124   DefineOutput(1, TH1F::Class()); // counter
125   DefineOutput(2, TList::Class()); // single track properties list and PID
126   DefineOutput(3, TList::Class()); // jet properties list
127
128 }
129
130 //________________________________________________________________________
131 AliAnalysisTaskSEHFCJqa::~AliAnalysisTaskSEHFCJqa(){
132   // destructor
133   delete fCuts;
134   delete fpidResp;
135   delete fhEventCounter;
136   delete fhSparseFilterMask;
137   delete   fhSparseFilterMaskTrackAcc;
138   delete fhSparseFilterMaskImpPar;
139   delete fhSparseEoverPeleTPC;
140   delete fhSparseShowShapeEleTPC;
141   delete fhnSigmaTPCTOFEle;
142   delete fhnSigmaTPCTOFPion;
143   delete fhnSigmaTPCTOFKaon;
144   delete fhnSigmaTPCTOFProton;
145   delete fhTrackEMCal;
146   delete fSparseRecoJets;
147   delete fhImpParResolITSsel;
148   delete fhImpParResolITSselGoodTracks;
149   delete fListTrackAndPID;
150   delete fListJets;
151 }
152
153 //________________________________________________________________________
154 void AliAnalysisTaskSEHFCJqa::Init()
155 {
156   // Initialization
157
158 }
159
160
161 //________________________________________________________________________
162 void AliAnalysisTaskSEHFCJqa::UserCreateOutputObjects(){
163
164
165   //##########  DEFINE THE TLISTS ##################
166   fListTrackAndPID=new TList();
167   fListTrackAndPID->SetOwner();
168   fListTrackAndPID->SetName("fListTrackAndPID");
169
170   fListJets=new TList();
171   fListJets->SetOwner();
172   fListJets->SetName("fListJets");
173
174   
175   //  ########### DEFINE THE EVENT COUNTER ############
176   fhEventCounter=new TH1F("fhEventCounter","Counter of event selected",20,-0.5,19.5);
177   fhEventCounter->GetXaxis()->SetBinLabel(1,"Events analyzed");
178   fhEventCounter->GetXaxis()->SetBinLabel(2,"Event selected");
179   fhEventCounter->GetXaxis()->SetBinLabel(3,"Jet array present");
180   fhEventCounter->GetXaxis()->SetBinLabel(4,"Vtx Track Ncont");
181
182   
183   
184   Int_t nbinsRecoJets[8]={50,20,20,20,5,5,10,60};
185   Double_t binlowRecoJets[8]={5.,-1.,-TMath::Pi(),0.99,0.,-0.5,0,0.};
186   Double_t binupRecoJets[8]={55.,1.,TMath::Pi(),20.99,4.99,4.5,2.,60.};
187   fSparseRecoJets=new THnSparseF("fSparseRecoJets","fSparseRecoJets;jetpt;eta;phi;ntrks;nEle;parton;partContr;ptPart;",8,nbinsRecoJets,binlowRecoJets,binupRecoJets);
188
189   fListJets->Add(fSparseRecoJets);
190
191   // Num axes: 10 filter bits + ID + TPCrefit,ITSrefit,bothTPCandITSrefit + kAny,kFirst,kBoth+ 20Nclust TPC+ 20 NTPC crossed padRows + 20 DCA
192   Int_t nbinsFilterMask[16]={2,2,2,2,2,2,2,2,2,2,2,3,3,20,20,70};
193   Double_t binlowFilterMask[16]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0,0,-3.5};
194   Double_t binupFilterMask[16]={1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,2.5,2.5,160,160,3.5};
195   fhSparseFilterMask=new THnSparseF("fhSparseFilterMask","fhSparseFilterMask;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackID;refitting;SPD;NTPCclust;NTPCcrossRows;DCA",16,nbinsFilterMask,binlowFilterMask,binupFilterMask);
196   fListTrackAndPID->Add(fhSparseFilterMask);
197
198
199 // Num axes: 10 filter bits + ID*isSelected 5 + kAny,kFirst,kBoth + phi+ eta +pt 
200   Int_t nbinsFilterMaskTrackAcc[15]={2,2,2,2,2,2,2,2,2,2,5,3,36,30,30};
201   Double_t binlowFilterMaskTrackAcc[15]={-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-2.5,-0.5,0.,-1.5,0.};
202   Double_t binupFilterMaskTrackAcc[15]={1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,1.5,2.5,2.5,TMath::Pi()*2.,1.5,15.};
203   fhSparseFilterMaskTrackAcc=new THnSparseF("fhSparseFilterMaskTrackAcc","fhSparseFilterMaskTrackAcc;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackIDandsel;SPD;#phi (rad);#eta;p_{T} (GeV/c)",15,nbinsFilterMaskTrackAcc,binlowFilterMaskTrackAcc,binupFilterMaskTrackAcc);
204   fListTrackAndPID->Add(fhSparseFilterMaskTrackAcc);  
205
206
207 // Num axes: ID*isSelected 5 + kAny,kFirst,kBoth + phi+ eta +pt + imp par
208   Int_t nbinsFilterMaskImpPar[6]={5,3,36,30,30,200};
209   Double_t binlowFilterMaskImpPar[6]={-2.5,-0.5,0.,-1.5,0.,-300.};
210   Double_t binupFilterMaskImpPar[6]={2.5,2.5,TMath::Pi()*2.,1.5,15.,300.};
211   fhSparseFilterMaskImpPar=new THnSparseF("fhSparseFilterMaskImpPar","fhSparseFilterMaskImpPar;fb0;fb1;fb2;fb3;fb4;fb5;fb6;fb7;fb8;fb9;trackIDandsel;SPD;#phi (rad);#eta;p_{T} (GeV/c; imp par (#mum)",6,nbinsFilterMaskImpPar,binlowFilterMaskImpPar,binupFilterMaskImpPar);
212   fListTrackAndPID->Add(fhSparseFilterMaskImpPar);  
213
214   
215 fhImpParResolITSsel=new TH3F("fhImpParResolITSsel","fhImpParResolITSsel;p_{T} (GeV/c);imp. par (#mum);ITS clust",50.,0.,10.,200,-800.,800.,38,-0.5,37.5);
216   // the convention for ITS clust axis: number between 0 and 48, first digit (units) = number of clust in ITS, second digits (x10): SPD status: 0 = none, 1=kFirst, 2=kSecond; 3= kBoth 
217   fListTrackAndPID->Add(fhImpParResolITSsel);
218
219   fhImpParResolITSselGoodTracks=new TH3F("fhImpParResolITSselGoodTracks","fhImpParResolITSselGoodTracks;p_{T} (GeV/c);imp. par (#mum);ITS clust",50.,0.,10.,200,-800.,800.,38,-0.5,37.5);
220   // the convention for ITS clust axis: number between 0 and 48, first digit (units) = number of clust in ITS, second digits (x10): SPD status: 0 = none, 1=kFirst, 2=kSecond; 3= kBoth 
221   fListTrackAndPID->Add(fhImpParResolITSselGoodTracks);
222
223
224   // PID PLOTS
225   fhnSigmaTPCTOFEle=new TH3F("fhnSigmaTPCTOFEle","fhnSigmaTPCTOFEle;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
226   fhnSigmaTPCTOFPion=new TH3F("fhnSigmaTPCTOFPion","fhnSigmaTPCTOFPion;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
227   fhnSigmaTPCTOFKaon=new TH3F("fhnSigmaTPCTOFKaon","fhnSigmaTPCTOFKaon;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
228   fhnSigmaTPCTOFProton=new TH3F("fhnSigmaTPCTOFProton","fhnSigmaTPCTOFProton;p (GeV/c);n#sigma_{TPC}^{ele};n#sigma_{TOF}^{ele}",120,0.,30,64,-8.,8.,64,-8.,8.);
229
230   fListTrackAndPID->Add(fhnSigmaTPCTOFEle);
231   fListTrackAndPID->Add(fhnSigmaTPCTOFPion);
232   fListTrackAndPID->Add(fhnSigmaTPCTOFKaon);
233   fListTrackAndPID->Add(fhnSigmaTPCTOFProton);
234
235
236   /// EMCAL PLOTS
237   //study of NsigmaTPC, e/p, p
238   Int_t nbinsEoP[4]={202,400,100,400};
239   Double_t binlowEoP[4]= {-1., -20.,-20., -1};
240   Double_t binupEoP[4]= {100., 20, 20.,9};
241   fhSparseEoverPeleTPC = new THnSparseF("fhSparseEoP", "fhSparseEoP; p;nsigmatpc;nsigmaElePIDresp;E/p;",4, nbinsEoP, binlowEoP, binupEoP);
242   fListTrackAndPID->Add(fhSparseEoverPeleTPC);
243   //fSparseEoverPallHadr = new THnSparseF("fSparseEoP", "fSparseEoP; p;nsigmatpc;E/p;",3, nbinsEoP, binlowEoP, binupEoP);
244
245   Int_t nbinsEmShSh[7]={35,120,100,50,50,50,50};
246   Double_t binlowEmShSh[7]= {5., -20., -1,0,0,0,0};
247   Double_t binupEmShSh[7]= {40., 10, 9,1,1,2,50};
248   fhSparseShowShapeEleTPC = new THnSparseF("fhSparseShowShapeEleTPC", "fhSparseShowShapeEleTPC; pt;nsigmatpc;E/p;M02;M20;disp;Nclust",7, nbinsEmShSh, binlowEmShSh, binupEmShSh);
249   fListTrackAndPID->Add(fhSparseShowShapeEleTPC);
250
251   
252   
253   Int_t nbinsTrEM[6]={124,20,124,25,20,20};
254   Double_t binlowTrEM[6]={-1,-1.,-1.0,0.,0.};
255   Double_t binupTrEM[6]={31,1.,31.,5.,1.,1.};
256
257   fhTrackEMCal=new THnSparseF("fhTrackEMCal","fhTrackEMCal;p(GeV/c);eta;clusterE(GeV/c);E/p;trkDistX;trkDistZ",6,nbinsTrEM,binlowTrEM,binupTrEM);
258   
259   
260
261   //fhPhotonicEle=new TH3F("fhPhotonicEle","fhPhotonicEle; pele;pgamma;mass",20,0.,20.,20,0.,20.,50,0.,5);
262   //fhPhotonicEleLS=new TH3F("fhPhotonicEleLS","fhPhotonicEleLS; pele;pgamma;mass",20,0.,20.,20,0.,20.,50,0.,5);
263
264
265   PostData(1,fhEventCounter);
266   PostData(2,fListTrackAndPID);
267   PostData(3,fListJets);
268   
269 }
270
271
272
273 //________________________________________________________________________
274 void AliAnalysisTaskSEHFCJqa::UserExec(Option_t */*option*/){
275
276
277
278 // Execute analysis for current event:
279   // heavy flavor candidates association to MC truth
280   
281   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
282   if(aod){
283     fhEventCounter->Fill(0);
284
285     if(!fCuts->IsEventSelected(aod)){
286       PostData(1,fhEventCounter);
287       return;
288     }
289     fhEventCounter->Fill(1);
290   }
291
292   TClonesArray *arrayJets=0x0;
293   if(!aod && AODEvent() && IsStandardAOD()) {
294     // In case there is an AOD handler writing a standard AOD, use the AOD 
295     // event in memory rather than the input (ESD) event.    
296     aod = dynamic_cast<AliAODEvent*> (AODEvent());
297     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
298     // have to taken from the AOD event hold by the AliAODExtension
299     if (!aod) {
300       Printf("ERROR: aod not available");
301       return;
302     }
303     else {
304       fhEventCounter->Fill(0);
305
306       if(!fCuts->IsEventSelected(aod)){
307         PostData(1,fhEventCounter);
308         return;
309       }
310       fhEventCounter->Fill(1);
311     }
312
313     AliAODHandler* aodHandler = (AliAODHandler*) 
314       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
315
316     if(fLoadJet>=1&&aodHandler->GetExtensions()) {
317       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.Jets.root");
318       AliAODEvent* aodFromExt = ext->GetAOD();
319       
320       // load jet array
321       arrayJets = (TClonesArray*)aodFromExt->GetList()->FindObject(fJetArrayString.Data());
322       if(!arrayJets) {
323         Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired branch %s not found!",fJetArrayString.Data());
324         PostData(1,fhEventCounter);
325         return;
326
327       }
328         
329     }
330   } else {
331     // load jet array
332     if(fLoadJet>=1){
333       arrayJets = (TClonesArray*)aod->GetList()->FindObject(fJetArrayString.Data());
334       if(!arrayJets) {
335         Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired Jet branch %s not found!",fJetArrayString.Data());
336         PostData(1,fhEventCounter);
337         return;
338       }
339       
340     }
341   }
342
343   if(fLoadJet!=0 && !arrayJets) {
344     printf("AliAnalysisTaskSEHFCJqa::UserExec: desired jet input branch not found!\n");
345     PostData(1,fhEventCounter);
346     return;
347   }
348
349   fhEventCounter->Fill(2);
350   
351   // AOD primary vertex
352   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
353   TString primTitle = vtx1->GetTitle();
354     if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) { 
355     fhEventCounter->Fill(3);
356
357   }
358   else {
359     PostData(1,fhEventCounter);
360     return;
361     
362   }
363     // Convert primary vertex to esd vertex (needed for later usage of RelateToVertex)
364     Double_t pos[3],cov[6];
365     vtx1->GetXYZ(pos);
366     vtx1->GetCovarianceMatrix(cov);
367     const AliESDVertex vESD(pos,cov,100.,100);
368     Double_t magfield=aod->GetMagneticField();
369
370   TClonesArray *arrayMC=0x0;
371   AliAODMCHeader *aodmcHeader=0x0;
372   Double_t vtxTrue[3];
373  
374   
375   if(fReadMC){
376     // load MC particles
377     arrayMC = 
378       (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
379     if(!arrayMC) {
380       Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC particles branch not found!\n");
381       return;
382     }
383     // load MC header
384     aodmcHeader = 
385       (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
386     if(!aodmcHeader) {
387       Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC header branch not found!\n");
388       return;
389     }
390     // MC primary vertex
391     aodmcHeader->GetVertex(vtxTrue);
392
393   }
394
395   // Starting the fun part
396   SetupPIDresponse();// this sets the pid reponse to fPPIDRespons; could also get from the cut object (AliAODPidHF::GetPidResponse)
397
398   // Looping over aod tracks
399   for(Int_t j=0;j<aod->GetNTracks();j++){
400
401     AliAODTrack *aodtrack=aod->GetTrack(j);
402     // CHECK FILTER MAPS
403     if(!FillTrackHistosAndSelectTrack(aodtrack,&vESD,magfield))continue;
404     //    if(j%100==0)  
405   
406     Double_t p=aodtrack->P();
407     Double_t eta=aodtrack->Eta();
408     // START PID: TPC
409
410     Double_t nsigmaEleTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kElectron);
411     Double_t nsigmaPionTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kPion);
412     Double_t nsigmaKaonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kKaon);
413     Double_t nsigmaProtonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kProton);
414
415
416     // TOF
417     Double_t nsigmaEleTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kElectron);
418     Double_t nsigmaPionTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kPion);
419     Double_t nsigmaKaonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kKaon);
420     Double_t nsigmaProtonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kProton);
421
422
423     fhnSigmaTPCTOFEle->Fill(p,nsigmaEleTPC,nsigmaEleTOF);
424     fhnSigmaTPCTOFPion->Fill(p,nsigmaPionTPC,nsigmaPionTOF);
425     fhnSigmaTPCTOFKaon->Fill(p,nsigmaKaonTPC,nsigmaKaonTOF);
426     fhnSigmaTPCTOFProton->Fill(p,nsigmaProtonTPC,nsigmaProtonTOF);
427     //    if(j%100==0)
428
429
430     // NOW EMCAL
431     // CHECK WHETHER THERE IS A EMCAL CLUSTER
432     Int_t nClsId = aodtrack->GetEMCALcluster();
433     if(nClsId >0) {
434       AliVCluster *cluster = aod->GetCaloCluster(nClsId);
435       Double_t clsE = cluster->E();
436       if(!cluster->IsEMCAL())       continue;
437       
438       // Check whether it is an electron candidate: do not reject here hadrons to allow QA checks of (E/p,nsigmaTPC,pt)
439       Double_t nEoverP = clsE/p;
440       Double_t eOverPpidResp;
441       Double_t showerShape[4];
442       Double_t nsigmaEleEMCal=fpidResp->NumberOfSigmasEMCAL(aodtrack,AliPID::kElectron,eOverPpidResp,showerShape);
443       Double_t poix[4]={p, nsigmaEleTPC, nsigmaEleEMCal, nEoverP};
444       fhSparseEoverPeleTPC->Fill(poix);
445       
446       
447       Double_t emcTrackDx=cluster->GetTrackDx();
448       Double_t emcTrackDz=cluster->GetTrackDz();
449       Double_t pointET[6]={p,eta,clsE,nEoverP,emcTrackDx,emcTrackDz};
450       fhTrackEMCal->Fill(pointET);    
451       
452       Double_t pointEmShSh[7]={aodtrack->Pt(), static_cast<Double_t>(nsigmaEleTPC), static_cast<Double_t>(nEoverP),cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(), static_cast<Double_t>(cluster->GetNCells())};
453
454       fhSparseShowShapeEleTPC->Fill(pointEmShSh);
455
456     }
457   }
458   
459   
460   // NOW LOOP OVER JETS
461   
462   Int_t nJets=arrayJets->GetEntries();//calcolo numero di jet nell'event
463   AliAODJet *jet;
464   for(Int_t jetcand=0;jetcand<nJets;jetcand++){
465     //    if(jetcand%100==0)
466
467     jet=(AliAODJet*)arrayJets->UncheckedAt(jetcand);
468        
469     Double_t contribution=0,ptpart=-1;
470     Int_t partonnat=0;
471     if(fReadMC){
472       
473       AliAODMCParticle *parton=IsMCJet(arrayMC,jet,contribution);
474       if(parton){
475         Int_t pdg=TMath::Abs(parton->PdgCode());
476         //printf("pdg parton: %d \n",pdg);
477         if(pdg==21)partonnat=1;
478         else if(pdg<4)partonnat=2;
479         else if(pdg==4)partonnat=3;
480         else if(pdg==5)partonnat=4;
481         ptpart=parton->Pt();
482       }
483       
484     }
485     
486     
487     FillJetRecoHisto(jet,partonnat,contribution,ptpart);
488   }
489
490   
491
492   PostData(1,fhEventCounter);
493   PostData(2,fListTrackAndPID);
494   PostData(3,fListJets);
495
496   return;
497 }
498
499
500
501 void AliAnalysisTaskSEHFCJqa::SetupPIDresponse(){
502
503   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
504   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler(); 
505   fpidResp=inputHandler->GetPIDResponse();
506   if(!fpidResp)AliFatal("No PID response could be set");
507 }
508
509 //_______________________________________________________________
510 Bool_t AliAnalysisTaskSEHFCJqa::FillTrackHistosAndSelectTrack(AliAODTrack *aodtr, const AliESDVertex *primary, const Double_t magfield){
511   
512   Bool_t retval=kTRUE;
513   // THnSparse for filter bits
514   Double_t point[16]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.,-999.};
515   Double_t pointAcc[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.};
516   Double_t pointImpPar[6]={999.,-999.,-999.,-999.,-999.,-999.};
517
518   for(Int_t j=0;j<10;j++){
519     if(aodtr->TestFilterBit(TMath::Power(2,j))){
520       point[j]=1;      
521       pointAcc[j]=1;
522     }
523   }
524   // check ID
525   Int_t trid=aodtr->GetID();
526   if(aodtr->GetID()>0)point[10]=1.;
527   else if(aodtr->GetID()>0)point[10]=0.;
528   if(aodtr->GetID()==0)point[10]=-1.;
529
530   Float_t iparxy,iparz;
531   Float_t pt=aodtr->Pt();
532   Int_t clustITS=aodtr->GetITSNcls();// for histo
533   AliESDtrack esdtrack(aodtr);
534   // needed to calculate impact parameters
535   
536
537   // check refit status
538   Int_t refit=-1;
539   UInt_t status = esdtrack.GetStatus();
540   if(status&AliESDtrack::kTPCrefit)refit+=1;
541   if(status&AliESDtrack::kITSrefit)refit+=2;
542   point[11]=refit;
543   // CHECK SPD
544   point[12]=-1;
545   if(aodtr->HasPointOnITSLayer(0)){
546     clustITS+=10;
547     point[12]+=1;
548   }
549   if(aodtr->HasPointOnITSLayer(1)){
550     point[12]+=2;
551     clustITS+=20;
552   }
553   point[13]=aodtr->GetTPCNcls();
554   point[14]=aodtr->GetTPCNCrossedRows();
555   esdtrack.RelateToVertex(primary,magfield,4.);// CHECK THIS : I put 4.. usually we set it to 3 
556   esdtrack.GetImpactParameters(iparxy,iparz);
557   point[15]=iparxy;
558
559
560
561   
562   fhSparseFilterMask->Fill(point);
563   if(!aodtr->TestBit(ffilterbit))retval =kFALSE;
564
565   if(aodtr->GetID()<0&&!fKeepTrackNegID)retval = kFALSE;
566   if(retval)  fhImpParResolITSsel->Fill(pt,iparxy*10000.,clustITS);
567
568   AliESDtrackCuts *cuts=fCuts->GetTrackCuts();
569   if(!cuts->IsSelected(&esdtrack)) retval = kFALSE;
570
571   if(fCuts->GetUseKinkRejection()){
572     AliAODVertex *maybeKink=aodtr->GetProdVertex();
573     if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
574   }
575
576   if(retval){
577     pointAcc[10]=1*trid;
578   }
579   else {
580     if(trid!=0)
581     pointAcc[10]=2*trid;
582     else pointAcc[10]=-3;
583   }
584
585   pointAcc[11]=point[12];
586   pointAcc[12]=aodtr->Phi();
587   if(pointAcc[12]<0.)pointAcc[12]+=2.*TMath::Pi();    
588   pointAcc[13]=aodtr->Eta();
589   pointAcc[14]=pt;
590   fhSparseFilterMaskTrackAcc->Fill(pointAcc);
591
592   pointImpPar[0]=pointAcc[10];
593   pointImpPar[1]=pointAcc[11];
594   pointImpPar[2]=pointAcc[12];
595   pointImpPar[3]=pointAcc[13];
596   pointImpPar[4]=pointAcc[14];
597   pointImpPar[5]=iparxy;
598   fhSparseFilterMaskImpPar->Fill(pointImpPar);
599   
600   if(retval)  {
601     fhImpParResolITSselGoodTracks->Fill(pt,iparxy*10000.,clustITS);  
602   }
603   
604   return retval;
605   
606 }
607
608
609 //---------------------------------------------------------------
610 AliAODMCParticle* AliAnalysisTaskSEHFCJqa::IsMCJet(TClonesArray *arrayMC,const AliAODJet *jet, Double_t &contribution){// assignment of parton ID to jet
611   // method by L. Feldkamp
612   std::vector< int >           idx;
613   std::vector< int >           idx2;
614   std::vector< double >     weight;
615
616   int counter =0;
617   int num = jet->GetRefTracks()->GetEntries();
618   
619   for(int k=0;k<num;++k){
620     
621     AliAODTrack * track = (AliAODTrack*)(jet->GetRefTracks()->At(k));
622     if (track->GetLabel() >=0){
623       AliAODMCParticle* part =  (AliAODMCParticle*)  arrayMC->At(track->GetLabel());
624       if(!part)continue;
625
626       int label =0 ;
627       AliAODMCParticle* motherParton=GetMCPartonOrigin(arrayMC,part, label);
628       if (!motherParton) Printf("no mother");
629       else {
630         counter++;
631         idx.push_back(label);                       //! Label  of Mother
632         idx2.push_back(label);        
633         weight.push_back(track->Pt());  //! Weight : P_t trak /  P_t jet ... the latter used at the end
634       }
635     }///END LOOP OVER REFERENCED TRACKS   
636   }
637   //! Remove duplicate indices for counting
638   sort( idx2.begin(), idx2.end() );
639   idx2.erase( unique( idx2.begin(), idx2.end() ), idx2.end() );
640   if (idx2.size() == 0) return 0x0;
641   Double_t* arrayOfWeights = new Double_t[(UInt_t)idx2.size()];
642   if(!arrayOfWeights){
643     return 0x0;
644   }
645   for(UInt_t ii=0;ii<(UInt_t)idx2.size();ii++)arrayOfWeights[ii]=0;
646
647   for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
648     for (unsigned int z=0; z< idx.size() ; ++z){
649       int     a = idx.at(z);
650       double w = weight.at(z);
651       if(a == idx2.at(idxloop))    arrayOfWeights[idxloop] += w;
652     }
653   }
654   
655   int winner = -1;
656   double c=-1.;
657   for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
658     if(c < arrayOfWeights[idxloop]){
659       winner =idxloop; 
660       c=arrayOfWeights[idxloop];
661     }
662   }
663   
664   AliAODMCParticle *parton = 0x0;
665   if(winner>0){
666     parton=(AliAODMCParticle*)arrayMC->At(idx.at(winner));
667     contribution = arrayOfWeights[winner]/jet->Pt();
668   }
669   else {
670   
671     if(arrayOfWeights)    delete [] arrayOfWeights;
672     return 0x0;
673   }
674   if(arrayOfWeights)    delete [] arrayOfWeights;
675
676   return parton;
677   
678   
679 }
680
681 //---------------------------------------------------------------
682 AliAODMCParticle *AliAnalysisTaskSEHFCJqa::GetMCPartonOrigin(TClonesArray *arrayMC,AliAODMCParticle *p, Int_t &idx)
683 {  //Follows chain of track mothers until q/g or idx = -1       
684   AliAODMCParticle *p2=0x0;
685   Int_t mlabel = TMath::Abs(p->GetMother()) ; 
686   Double_t pz=0.;
687   while(mlabel > 1){
688     p2 = (AliAODMCParticle*)arrayMC->At(mlabel);
689     pz=TMath::Abs(p2->Pz());
690     //printf("Mother label %d, pdg %d, pz %f\n",mlabel,p2->PdgCode(),pz);
691     if( p2->PdgCode() == 21 || (p2->PdgCode() != 0 && abs(p2->PdgCode()) <6) )
692       {
693         idx = mlabel; 
694         return p2;
695       }
696     mlabel = TMath::Abs(p2->GetMother()); 
697   }
698   idx=-1;
699   return 0x0;
700
701
702
703 //_____________________________________________________________
704 void AliAnalysisTaskSEHFCJqa::FillJetRecoHisto(const AliAODJet *jet,Int_t partonnat,Double_t contribution,Double_t ptpart){
705 //FIll sparse with reco jets properties
706   Double_t point[8]={jet->Pt(),jet->Eta(),jet->Phi()-TMath::Pi(), static_cast<Double_t>(jet->GetRefTracks()->GetEntriesFast()),0, static_cast<Double_t>(partonnat),contribution,ptpart};
707   fSparseRecoJets->Fill(point);
708 }
709
710
711 //_______________________________________________________________
712 //void AliAnalysisTaskSEHFCJqa::FillTrackHistosPID(AliAODTrack *aodtr){
713 //
714 //}
715
716 //_______________________________________________________________
717 void AliAnalysisTaskSEHFCJqa::Terminate(const Option_t*){
718   //TERMINATE METHOD: NOTHING TO DO
719
720
721
722 }
723