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(),nsigmaEleTPC,nEoverP,cluster->GetM02(),cluster->GetM20(),cluster->GetDispersion(),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(),jet->GetRefTracks()->GetEntriesFast(),0,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