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     Printf("ERROR: aod not available");
 284     return;
 285   }
 286   fhEventCounter->Fill(0);
 287
 288   if(!fCuts->IsEventSelected(aod)){
 289     PostData(1,fhEventCounter);
 290     return;
 291   }
 292   fhEventCounter->Fill(1);
 293
 294   TClonesArray *arrayJets;
 295   if(!aod && AODEvent() && IsStandardAOD()) {
 296     // In case there is an AOD handler writing a standard AOD, use the AOD
 297     // event in memory rather than the input (ESD) event.
 298     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 299     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 300     // have to taken from the AOD event hold by the AliAODExtension
 301     AliAODHandler* aodHandler = (AliAODHandler*)
 302       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 303
 304     if(fLoadJet>=1&&aodHandler->GetExtensions()) {
 305       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.Jets.root");
 306       AliAODEvent* aodFromExt = ext->GetAOD();
 307
 308       // load jet array
 309       arrayJets = (TClonesArray*)aodFromExt->GetList()->FindObject(fJetArrayString.Data());
 310       if(!arrayJets) {
 311         Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired branch %s not found!",fJetArrayString.Data());
 312         PostData(1,fhEventCounter);
 313         return;
 314
 315       }
 316
 317     }
 318   } else {
 319     // load jet array
 320     if(fLoadJet>=1){
 321       arrayJets = (TClonesArray*)aod->GetList()->FindObject(fJetArrayString.Data());
 322       if(!arrayJets) {
 323         Printf("AliAnalysisTaskSEHFCJqa::UserExec: desired Jet branch %s not found!",fJetArrayString.Data());
 324         PostData(1,fhEventCounter);
 325         return;
 326       }
 327
 328     }
 329   }
 330
 331   if(fLoadJet!=0 && !arrayJets) {
 332     printf("AliAnalysisTaskSEHFCJqa::UserExec: desired jet input branch not found!\n");
 333     PostData(1,fhEventCounter);
 334     return;
 335   }
 336
 337   fhEventCounter->Fill(2);
 338
 339   // AOD primary vertex
 340   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 341   TString primTitle = vtx1->GetTitle();
 342     if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
 343     fhEventCounter->Fill(3);
 344
 345   }
 346   else {
 347     PostData(1,fhEventCounter);
 348     return;
 349
 350   }
 351     // Convert primary vertex to esd vertex (needed for later usage of RelateToVertex)
 352     Double_t pos[3],cov[6];
 353     vtx1->GetXYZ(pos);
 354     vtx1->GetCovarianceMatrix(cov);
 355     const AliESDVertex vESD(pos,cov,100.,100);
 356     Double_t magfield=aod->GetMagneticField();
 357
 358   TClonesArray *arrayMC=0x0;
 359   AliAODMCHeader *aodmcHeader=0x0;
 360   Double_t vtxTrue[3];
 361
 362
 363   if(fReadMC){
 364     // load MC particles
 365     arrayMC =
 366       (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 367     if(!arrayMC) {
 368       Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC particles branch not found!\n");
 369       return;
 370     }
 371     // load MC header
 372     aodmcHeader =
 373       (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 374     if(!aodmcHeader) {
 375       Printf("AliAnalysisTaskSEHFCJqa::UserExec: MC header branch not found!\n");
 376       return;
 377     }
 378     // MC primary vertex
 379     aodmcHeader->GetVertex(vtxTrue);
 380
 381   }
 382
 383   // Starting the fun part
 384   SetupPIDresponse();// this sets the pid reponse to fPPIDRespons; could also get from the cut object (AliAODPidHF::GetPidResponse)
 385
 386   // Looping over aod tracks
 387   for(Int_t j=0;j<aod->GetNTracks();j++){
 388
 389     AliAODTrack *aodtrack=aod->GetTrack(j);
 390     // CHECK FILTER MAPS
 391     if(!FillTrackHistosAndSelectTrack(aodtrack,vESD,magfield))continue;
 392     //    if(j%100==0)
 393
 394     Double_t p=aodtrack->P();
 395     Double_t eta=aodtrack->Eta();
 396     // START PID: TPC
 397
 398     Double_t nsigmaEleTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kElectron);
 399     Double_t nsigmaPionTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kPion);
 400     Double_t nsigmaKaonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kKaon);
 401     Double_t nsigmaProtonTPC=fpidResp->NumberOfSigmasTPC(aodtrack, AliPID::kProton);
 402
 403
 404     // TOF
 405     Double_t nsigmaEleTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kElectron);
 406     Double_t nsigmaPionTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kPion);
 407     Double_t nsigmaKaonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kKaon);
 408     Double_t nsigmaProtonTOF=fpidResp->NumberOfSigmasTOF(aodtrack, AliPID::kProton);
 409
 410
 411     fhnSigmaTPCTOFEle->Fill(p,nsigmaEleTPC,nsigmaEleTOF);
 412     fhnSigmaTPCTOFPion->Fill(p,nsigmaPionTPC,nsigmaPionTOF);
 413     fhnSigmaTPCTOFKaon->Fill(p,nsigmaKaonTPC,nsigmaKaonTOF);
 414     fhnSigmaTPCTOFProton->Fill(p,nsigmaProtonTPC,nsigmaProtonTOF);
 415     //    if(j%100==0)
 416
 417
 418     // NOW EMCAL
 419     // CHECK WHETHER THERE IS A EMCAL CLUSTER
 420     Int_t nClsId = aodtrack->GetEMCALcluster();
 421     if(nClsId >0) {
 422       AliVCluster *cluster = aod->GetCaloCluster(nClsId);
 423       Double_t clsE = cluster->E();
 424       if(!cluster->IsEMCAL())       continue;
 425
 426       // Check whether it is an electron candidate: do not reject here hadrons to allow QA checks of (E/p,nsigmaTPC,pt)
 427       Double_t nEoverP = clsE/p;
 428       Double_t eOverPpidResp;
 429       Double_t showerShape[4];
 430       Double_t nsigmaEleEMCal=fpidResp->NumberOfSigmasEMCAL(aodtrack,AliPID::kElectron,eOverPpidResp,showerShape);
 431       Double_t poix[4]={p, nsigmaEleTPC, nsigmaEleEMCal, nEoverP};
 432       fhSparseEoverPeleTPC->Fill(poix);
 433
 434
 435       Double_t emcTrackDx=cluster->GetTrackDx();
 436       Double_t emcTrackDz=cluster->GetTrackDz();
 437       Double_t pointET[6]={p,eta,clsE,nEoverP,emcTrackDx,emcTrackDz};
 438       fhTrackEMCal->Fill(pointET);
 439
 440       Double_t pointEmShSh[7]={aodtrack->Pt(),nsigmaEleTPC,nEoverP,cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),cluster->GetNCells()};
 441
 442       fhSparseShowShapeEleTPC->Fill(pointEmShSh);
 443
 444     }
 445   }
 446
 447
 448   // NOW LOOP OVER JETS
 449
 450   Int_t nJets=arrayJets->GetEntries();//calcolo numero di jet nell'event
 451   AliAODJet *jet;
 452   for(Int_t jetcand=0;jetcand<nJets;jetcand++){
 453     //    if(jetcand%100==0)
 454
 455     jet=(AliAODJet*)arrayJets->UncheckedAt(jetcand);
 456
 457     Double_t contribution=0,ptpart=-1;
 458     Int_t partonnat=0;
 459     if(fReadMC){
 460
 461       AliAODMCParticle *parton=IsMCJet(arrayMC,jet,contribution);
 462       if(parton){
 463         Int_t pdg=TMath::Abs(parton->PdgCode());
 464         //printf("pdg parton: %d \n",pdg);
 465         if(pdg==21)partonnat=1;
 466         else if(pdg<4)partonnat=2;
 467         else if(pdg==4)partonnat=3;
 468         else if(pdg==5)partonnat=4;
 469         ptpart=parton->Pt();
 470       }
 471
 472     }
 473
 474
 475     FillJetRecoHisto(jet,partonnat,contribution,ptpart);
 476   }
 477
 478
 479
 480   PostData(1,fhEventCounter);
 481   PostData(2,fListTrackAndPID);
 482   PostData(3,fListJets);
 483
 484   return;
 485 }
 486
 487
 488
 489 void AliAnalysisTaskSEHFCJqa::SetupPIDresponse(){
 490
 491   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 492   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 493   fpidResp=inputHandler->GetPIDResponse();
 494   if(!fpidResp)AliFatal("No PID response could be set");
 495 }
 496
 497 //_______________________________________________________________
 498 Bool_t AliAnalysisTaskSEHFCJqa::FillTrackHistosAndSelectTrack(AliAODTrack *aodtr, const AliESDVertex primary, const Double_t magfield){
 499
 500   Bool_t retval=kTRUE;
 501   // THnSparse for filter bits
 502   Double_t point[16]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.,-999.};
 503   Double_t pointAcc[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,999.,-999.,-999.,-999.,-999.};
 504   Double_t pointImpPar[6]={999.,-999.,-999.,-999.,-999.,-999.};
 505
 506   for(Int_t j=0;j<10;j++){
 507     if(aodtr->TestFilterBit(TMath::Power(2,j))){
 508       point[j]=1;
 509       pointAcc[j]=1;
 510     }
 511   }
 512   // check ID
 513   Int_t trid=aodtr->GetID();
 514   if(aodtr->GetID()>0)point[10]=1.;
 515   else if(aodtr->GetID()>0)point[10]=0.;
 516   if(aodtr->GetID()==0)point[10]=-1.;
 517
 518   Float_t iparxy,iparz;
 519   Float_t pt=aodtr->Pt();
 520   Int_t clustITS=aodtr->GetITSNcls();// for histo
 521   AliESDtrack esdtrack(aodtr);
 522   // needed to calculate impact parameters
 523
 524
 525   // check refit status
 526   Int_t refit=-1;
 527   UInt_t status = esdtrack.GetStatus();
 528   if(status&AliESDtrack::kTPCrefit)refit+=1;
 529   if(status&AliESDtrack::kITSrefit)refit+=2;
 530   point[11]=refit;
 531   // CHECK SPD
 532   point[12]=-1;
 533   if(aodtr->HasPointOnITSLayer(0)){
 534     clustITS+=10;
 535     point[12]+=1;
 536   }
 537   if(aodtr->HasPointOnITSLayer(1)){
 538     point[12]+=2;
 539     clustITS+=20;
 540   }
 541   point[13]=aodtr->GetTPCNcls();
 542   point[14]=aodtr->GetTPCNCrossedRows();
 543   esdtrack.RelateToVertex(&primary,magfield,4.);// CHECK THIS : I put 4.. usually we set it to 3
 544   esdtrack.GetImpactParameters(iparxy,iparz);
 545   point[15]=iparxy;
 546
 547
 548
 549
 550   fhSparseFilterMask->Fill(point);
 551   if(!aodtr->TestBit(ffilterbit))retval =kFALSE;
 552
 553   if(aodtr->GetID()<0&&!fKeepTrackNegID)retval = kFALSE;
 554   if(retval)  fhImpParResolITSsel->Fill(pt,iparxy*10000.,clustITS);
 555
 556   AliESDtrackCuts *cuts=fCuts->GetTrackCuts();
 557   if(!cuts->IsSelected(&esdtrack)) retval = kFALSE;
 558
 559   if(fCuts->GetUseKinkRejection()){
 560     AliAODVertex *maybeKink=aodtr->GetProdVertex();
 561     if(maybeKink->GetType()==AliAODVertex::kKink) retval=kFALSE;
 562   }
 563
 564   if(retval){
 565     pointAcc[10]=1*trid;
 566   }
 567   else {
 568     if(trid!=0)
 569     pointAcc[10]=2*trid;
 570     else pointAcc[10]=-3;
 571   }
 572
 573   pointAcc[11]=point[12];
 574   pointAcc[12]=aodtr->Phi();
 575   if(pointAcc[12]<0.)pointAcc[12]+=2.*TMath::Pi();
 576   pointAcc[13]=aodtr->Eta();
 577   pointAcc[14]=pt;
 578   fhSparseFilterMaskTrackAcc->Fill(pointAcc);
 579
 580   pointImpPar[0]=pointAcc[10];
 581   pointImpPar[1]=pointAcc[11];
 582   pointImpPar[2]=pointAcc[12];
 583   pointImpPar[3]=pointAcc[13];
 584   pointImpPar[4]=pointAcc[14];
 585   pointImpPar[5]=iparxy;
 586   fhSparseFilterMaskImpPar->Fill(pointImpPar);
 587
 588   if(retval)  {
 589     fhImpParResolITSselGoodTracks->Fill(pt,iparxy*10000.,clustITS);
 590   }
 591
 592   return retval;
 593
 594 }
 595
 596
 597 //---------------------------------------------------------------
 598 AliAODMCParticle* AliAnalysisTaskSEHFCJqa::IsMCJet(TClonesArray *arrayMC,const AliAODJet *jet, Double_t &contribution){// assignment of parton ID to jet
 599   // method by L. Feldkamp
 600   std::vector< int >           idx;
 601   std::vector< int >           idx2;
 602   std::vector< double >     weight;
 603
 604   int counter =0;
 605   int num = jet->GetRefTracks()->GetEntries();
 606
 607   for(int k=0;k<num;++k){
 608
 609     AliAODTrack * track = (AliAODTrack*)(jet->GetRefTracks()->At(k));
 610     if (track->GetLabel() >=0){
 611       AliAODMCParticle* part =  (AliAODMCParticle*)  arrayMC->At(track->GetLabel());
 612       if(!part)continue;
 613
 614       int label =0 ;
 615       AliAODMCParticle* motherParton=GetMCPartonOrigin(arrayMC,part, label);
 616       if (!motherParton) Printf("no mother");
 617       else {
 618         counter++;
 619         idx.push_back(label);                       //! Label  of Mother
 620         idx2.push_back(label);
 621         weight.push_back(track->Pt());  //! Weight : P_t trak /  P_t jet ... the latter used at the end
 622       }
 623     }///END LOOP OVER REFERENCED TRACKS
 624   }
 625   //! Remove duplicate indices for counting
 626   sort( idx2.begin(), idx2.end() );
 627   idx2.erase( unique( idx2.begin(), idx2.end() ), idx2.end() );
 628   if (idx2.size() == 0) return 0x0;
 629   Double_t* arrayOfWeights = new Double_t [(int)idx2.size()];
 630   for(Int_t ii=0;ii<(Int_t)idx2.size();ii++)arrayOfWeights[ii]=0;
 631
 632   for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
 633     for (unsigned int z=0; z< idx.size() ; ++z){
 634       int     a = idx.at(z);
 635       double w = weight.at(z);
 636       if(a == idx2.at(idxloop))    arrayOfWeights[idxloop] += w;
 637     }
 638   }
 639
 640   int winner = -1;
 641   double c=-1.;
 642   for (unsigned int idxloop =0 ;idxloop<idx2.size();idxloop++){
 643     if(c < arrayOfWeights[idxloop]){
 644       winner =idxloop;
 645       c=arrayOfWeights[idxloop];
 646     }
 647   }
 648
 649   AliAODMCParticle *parton = 0x0;
 650   parton=(AliAODMCParticle*)arrayMC->At(idx.at(winner));
 651   contribution = arrayOfWeights[winner]/jet->Pt();
 652
 653   if(arrayOfWeights)    delete arrayOfWeights;
 654
 655
 656   return parton;
 657
 658
 659 }
 660
 661 //---------------------------------------------------------------
 662 AliAODMCParticle *AliAnalysisTaskSEHFCJqa::GetMCPartonOrigin(TClonesArray *arrayMC,AliAODMCParticle *p, Int_t &idx)
 663 {  //Follows chain of track mothers until q/g or idx = -1
 664   AliAODMCParticle *p2=0x0;
 665   Int_t mlabel = TMath::Abs(p->GetMother()) ;
 666   Double_t pz=0.;
 667   while(mlabel > 1){
 668     p2 = (AliAODMCParticle*)arrayMC->At(mlabel);
 669     pz=TMath::Abs(p2->Pz());
 670     //printf("Mother label %d, pdg %d, pz %f\n",mlabel,p2->PdgCode(),pz);
 671     if( p2->PdgCode() == 21 || (p2->PdgCode() != 0 && abs(p2->PdgCode()) <6) )
 672       {
 673         idx = mlabel;
 674         return p2;
 675       }
 676     mlabel = TMath::Abs(p2->GetMother());
 677   }
 678   idx=-1;
 679   return 0x0;
 680
 681 }
 682
 683 //_____________________________________________________________
 684 void AliAnalysisTaskSEHFCJqa::FillJetRecoHisto(const AliAODJet *jet,Int_t partonnat,Double_t contribution,Double_t ptpart){
 685 //FIll sparse with reco jets properties
 686   Double_t point[8]={jet->Pt(),jet->Eta(),jet->Phi()-TMath::Pi(),jet->GetRefTracks()->GetEntriesFast(),0,partonnat,contribution,ptpart};
 687   fSparseRecoJets->Fill(point);
 688 }
 689
 690
 691 //_______________________________________________________________
 692 //void AliAnalysisTaskSEHFCJqa::FillTrackHistosPID(AliAODTrack *aodtr){
 693 //
 694 //}
 695
 696 //_______________________________________________________________
 697 void AliAnalysisTaskSEHFCJqa::Terminate(const Option_t*){
 698   //TERMINATE METHOD: NOTHING TO DO
 699
 700
 701
 702 }
 703