PWG4/JetTasks/AliAnalysisTaskCorrectionsUE.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id:$ */
  17
  18 ////////////////////////////////////////////////////////////////////////
  19 //
  20 // Analysis class to Correct Underlying Event studies
  21 //
  22 // This class needs as input ESDs.\r
  23 // The output is an analysis-specific container.\r
  24 //
  25 // The class is used to get the contamination from secondaries\r
  26 // from tracks DCA distribution \r
  27 // as function of track pT and pseudo-rapidity.\r
  28 // It provides additional information for the corrections \r
  29 // that can not be retrieved by the AliAnalysisTaskLeadingTackUE\r
  30 // task, which is running on AODs.\r
  31 //
  32 ////////////////////////////////////////////////////////////////////////
  33
  34 #include <TROOT.h>
  35 #include <TChain.h>
  36 //#include <TCanvas.h>\r
  37 //#include <TFile.h>\r
  38 #include <TList.h>
  39 #include <TMath.h>
  40 //#include <TProfile.h>\r
  41 #include <TTree.h>
  42 //#include <TVector3.h>\r
  43 #include <TH3F.h>\r
  44
  45 #include "AliAnalyseLeadingTrackUE.h"\r
  46 #include "AliAnalysisFilter.h"\r
  47 #include "AliAnalysisHelperJetTasks.h"
  48 #include "AliAnalysisManager.h"
  49 #include "AliAnalysisTaskCorrectionsUE.h"
  50 #include "AliAODHandler.h"
  51 #include "AliESDHandler.h"
  52 #include "AliESDtrack.h"
  53 #include "AliESDtrackCuts.h"\r
  54 #include "AliESDEvent.h"
  55 #include "AliAODInputHandler.h"
  56 #include "AliESDInputHandler.h"
  57 #include "AliCFContainer.h"
  58 #include "AliCFManager.h"
  59 #include "AliGenDPMjetEventHeader.h"
  60 #include "AliGenPythiaEventHeader.h"
  61 #include "AliInputEventHandler.h"
  62 #include "AliLog.h"
  63 #include "AliMCEventHandler.h"
  64 #include "AliAnalysisHelperJetTasks.h"
  65
  66 class TCanvas;\r
  67 class TFile;\r
  68 class TProfile;\r
  69 class TVector3; \r
  70 \r
  71 ClassImp( AliAnalysisTaskCorrectionsUE)
  72
  73 // Define global pointer
  74 AliAnalysisTaskCorrectionsUE* AliAnalysisTaskCorrectionsUE::fgTaskCorrectionsUE=NULL;
  75
  76 //____________________________________________________________________
  77 AliAnalysisTaskCorrectionsUE:: AliAnalysisTaskCorrectionsUE(const char* name):
  78 AliAnalysisTask(name,""),
  79 fAnalyseUE(0x0),\r
  80 fDebug(0),\r
  81 fESDEvent(0x0),\r
  82 fESDHandler(0x0),
  83 fInputHandler(0x0),
  84 fListOfHistos(0x0),
  85 fMcEvent(0x0),\r
  86 fMcHandler(0x0),
  87 fMode(0),\r
  88 fOutCFcont(0x0),\r
  89 fhEntries(0x0),\r
  90 fhFakes(0x0),\r
  91 fhPtMCAll(0x0),\r
  92 fhPtMCPrim(0x0),\r
  93 fhPtMCSec(0x0),\r
  94 fhPtMCPrimFake(0x0),\r
  95 fhPtMCSecFake(0x0),\r
  96 fnTracksVertex(1),  // QA tracks pointing to principal vertex  \r
  97 fZVertex(10.),\r
  98 fhVertexContributors(0x0),\r
  99 fhVertexReso(0x0),\r
 100 fTrackEtaCut(0.9),
 101 fTrackPtCut(0.),
 102 fEsdTrackCuts(0x0),\r
 103 fEsdTrackCutsSPD(0x0),\r
 104 fEsdTrackCutsSDD(0x0),\r
 105 fEsdTrackCutsDCA(0x0)\r
 106 {
 107   // Default constructor
 108   // Define input and output slots here
 109   // Input slot #0 works with a TChain
 110   DefineInput(0, TChain::Class());
 111   // Output slot #0 writes into a TList container
 112   DefineOutput(0, TList::Class());
 113
 114 }
 115
 116 //______________________________________________________________
 117 AliAnalysisTaskCorrectionsUE & AliAnalysisTaskCorrectionsUE::operator = (const AliAnalysisTaskCorrectionsUE & /*source*/)
 118 {
 119   // assignment operator
 120   return *this;
 121 }
 122
 123 \r
 124 /************** INTERFACE METHODS *****************************/\r
 125 \r
 126 //______________________________________________________________
 127 Bool_t AliAnalysisTaskCorrectionsUE::Notify()
 128 {
 129
 130   return kTRUE;
 131
 132 }
 133
 134 //____________________________________________________________________
 135 void AliAnalysisTaskCorrectionsUE::ConnectInputData(Option_t* /*option*/)
 136 {
 137   // Connect the input data
 138
 139   if (fDebug > 1) AliInfo("ConnectInputData() ");
 140
 141   //Get the input handler
 142   fESDHandler = (AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();\r
 143   if ( !fESDHandler && fDebug > 0 ) {\r
 144         AliFatal(" No ESD event handler connected !!! ");\r
 145         return;\r
 146         }\r
 147
 148   //Get ESD event\r
 149    fESDEvent = (AliESDEvent*)fESDHandler->GetEvent();\r
 150    if (!fESDEvent && fDebug > 1) {\r
 151         AliFatal(" No ESD event retrieved !!! ");\r
 152         return;\r
 153         }\r
 154
 155   //Get MC handler
 156   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 157   \r
 158   // Define track cuts\r
 159   fEsdTrackCuts = new AliESDtrackCuts("ITSTPC", "ITS+TPC standard 2009 cuts w.o. SPD cluster requirement nor DCA cut");\r
 160   // TPC  \r
 161   fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin\r
 162   fEsdTrackCuts->SetMinNClustersTPC(70);\r
 163   //fEsdTrackCuts->SetMinNClustersTPC(90); // ***** TMP *****\r
 164   fEsdTrackCuts->SetMaxChi2PerClusterTPC(4);\r
 165   fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);\r
 166   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);\r
 167   // ITS\r
 168   fEsdTrackCuts->SetRequireITSRefit(kTRUE);\r
 169   fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
 170   fEsdTrackCuts->SetMaxDCAToVertexZ(2); // new for pile-up !!!\r
 171 \r
 172   // Add SPD requirement \r
 173   fEsdTrackCutsSPD = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");\r
 174   fEsdTrackCutsSPD->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
 175   \r
 176   // Add SDD requirement \r
 177   fEsdTrackCutsSDD = new AliESDtrackCuts("SDD", "Require 1 cluster in first layer SDD");\r
 178   fEsdTrackCutsSDD->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kFirst);\r
 179   \r
 180   // Add DCA cuts \r
 181   fEsdTrackCutsDCA = new AliESDtrackCuts("DCA", "pT dependent DCA cut");\r
 182   // 7*(0.0050+0.0060/pt^0.9)\r
 183   fEsdTrackCutsDCA->SetMaxDCAToVertexXYPtDep("0.0350+0.0420/pt^0.9");\r
 184   \r
 185   fEsdTrackCutsDCA->SetMaxDCAToVertexZ(1.e6);\r
 186   fEsdTrackCutsDCA->SetDCAToVertex2D(kFALSE);\r
 187 \r
 188   // emulates filterbit when getting leading-track from ESD\r
 189   fAnalyseUE->DefineESDCuts(16); // any number = standard ITS+TPC + (SPD or SDD)\r
 190 }
 191
 192 //____________________________________________________________________
 193 void  AliAnalysisTaskCorrectionsUE::CreateOutputObjects()
 194 {
 195   // Create the output container
 196
 197   if (fDebug > 1) AliInfo("CreateOutPutData()");
 198
 199   // Create pointer to AliAnalyseLeadingTrackUE, a class implementing the main analysis algorithms\r
 200   fAnalyseUE = new AliAnalyseLeadingTrackUE();\r
 201   if (!fAnalyseUE){\r
 202      AliError("UE analysis class not initialized!!!");
 203      return;
 204   }
 205   \r
 206   // Create list of output histograms
 207   if (fListOfHistos != NULL){
 208         delete fListOfHistos;
 209         fListOfHistos = NULL;
 210         }
 211   if (!fListOfHistos){
 212         fListOfHistos = new TList();
 213         fListOfHistos->SetOwner(kTRUE);
 214         }
 215
 216   // Create CF container\r
 217   CreateContainer();\r
 218   // number of events\r
 219   fhEntries = new TH1F("fhEntries","Entries",1,0.,2.); \r
 220   fListOfHistos->Add(fhEntries);\r
 221   // tracks contributing to the vertex\r
 222   fhVertexContributors = new TH1F("fhVertexContributors","Tracks in vertex",51, -0.5, 50.5);\r
 223   fhVertexContributors->GetXaxis()->SetTitle("# tracks in vertex");\r
 224   fhVertexContributors->GetYaxis()->SetTitle("Entries");\r
 225   fListOfHistos->Add(fhVertexContributors);\r
 226   // vertex resolution\r
 227   fhVertexReso = new TH3F("fhVertexReso","Vertex resolution",51,-0.5,50.5,100,0.,0.05,100.,0.,0.1); \r
 228   fhVertexContributors->GetXaxis()->SetTitle("# tracks in vertex");\r
 229   fhVertexContributors->GetYaxis()->SetTitle("Resolution XY (cm)");\r
 230   fhVertexContributors->GetZaxis()->SetTitle("Resolution Z (cm)");\r
 231   fListOfHistos->Add(fhVertexReso);\r
 232
 233   // number of fake tracks\r
 234   fhFakes = new TH1F("fhFakes","Fraction of fake tracks",5,-0.5,4.5);\r
 235   fhFakes->GetXaxis()->SetBinLabel(1,"No MC");\r
 236   fhFakes->GetXaxis()->SetBinLabel(2,"Unique MC primary");\r
 237   fhFakes->GetXaxis()->SetBinLabel(3,"Unique MC secondary");\r
 238   fhFakes->GetXaxis()->SetBinLabel(4,"Multiple MC");\r
 239   fListOfHistos->Add(fhFakes);\r
 240   \r
 241   //pT distributions\r
 242   fhPtMCAll = new TH1F("fhPtMCAll","All MC particles reconstructed",100,0., 20.);\r
 243   fListOfHistos->Add(fhPtMCAll);\r
 244   fhPtMCPrim = new TH1F("fhPtMCPrim","Primary MC particles reconstructed",100,0., 20.);\r
 245   fListOfHistos->Add(fhPtMCPrim);\r
 246   fhPtMCSec = new TH1F("fhPtMCSec","Secondary MC particles reconstructed",100,0., 20.);\r
 247   fListOfHistos->Add(fhPtMCSec);\r
 248   fhPtMCPrimFake = new TH1F("fhPtMCPrimFake","Fake primary MC particles reconstructed",100,0., 20.);\r
 249   fListOfHistos->Add(fhPtMCPrimFake);\r
 250   fhPtMCSecFake = new TH1F("fhPtMCSecFake","Fake secondary MC particles reconstructed",100,0., 20.);\r
 251   fListOfHistos->Add(fhPtMCSecFake);\r
 252 \r
 253 \r
 254   // Add task configuration to output list \r
 255   AddSettingsTree();
 256     \r
 257
 258   //Post outputs
 259   PostData(0,fListOfHistos);
 260  \r
 261 }
 262
 263 //____________________________________________________________________
 264 void  AliAnalysisTaskCorrectionsUE::Exec(Option_t */*option*/)
 265 {
 266   \r
 267   // Get MC event\r
 268   if (fMcHandler){
 269         fMcEvent = fMcHandler->MCEvent();
 270         if ( fDebug > 3 ) AliInfo( " Processing MC event..." );
 271         if (fMode && !fMcEvent) return;\r
 272         }
 273
 274   // Do the analysis\r
 275   AnalyseCorrectionMode();\r
 276
 277   \r
 278  PostData(0,fListOfHistos);\r
 279
 280 }\r
 281 \r
 282 //____________________________________________________________________\r
 283 void  AliAnalysisTaskCorrectionsUE::Terminate(Option_t */*option*/)\r
 284 {\r
 285   // Terminate analysis\r
 286
 287   if (fDebug >1) AliAnalysisHelperJetTasks::PrintDirectorySize("PWG4_JetTasksOutput.root");\r
 288   \r
 289   if (!gROOT->IsBatch()){\r
 290         fListOfHistos = dynamic_cast<TList*> (GetOutputData(0));\r
 291         if (!fListOfHistos){\r
 292                 AliError("Histogram List is not available");\r
 293                 return;\r
 294                 }
 295 \r
 296   } else {\r
 297         AliInfo(" Batch mode, not histograms will be shown...");\r
 298   }\r
 299 \r
 300   \r
 301 }
 302
 303 \r
 304 /******************** ANALYSIS METHODS *****************************/\r
 305 \r
 306 //____________________________________________________________________
 307 void  AliAnalysisTaskCorrectionsUE::AddSettingsTree()
 308 {
 309   //Write settings to output list
 310   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");\r
 311   settingsTree->Branch("fnTracksVertex", &fnTracksVertex, "TracksInVertex/D");\r
 312   settingsTree->Branch("fZVertex", &fZVertex, "VertexZCut/D");\r
 313   settingsTree->Branch("fTrackPtCut", &fTrackPtCut, "TrackPtCut/D");
 314   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
 315   settingsTree->Fill();
 316   fListOfHistos->Add(settingsTree);
 317 }  \r
 318
 319 //____________________________________________________________________
 320 void AliAnalysisTaskCorrectionsUE::AnalyseCorrectionMode()\r
 321 {
 322
 323   // Analyse the event\r
 324   Int_t labelMC  = -1;\r
 325   if (fMcHandler && fMcEvent){\r
 326         // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex\r
 327         if (!fAnalyseUE->VertexSelection(fMcEvent, 0, fZVertex))  \r
 328         return; \r
 329         // Get MC-true leading particle \r
 330         TObjArray *ltMC = (TObjArray*)fAnalyseUE->FindLeadingObjects(fMcEvent);\r
 331         AliVParticle* leadingMC = 0;\r
 332         if (ltMC){\r
 333                 leadingMC = (AliVParticle*) ltMC->At(0);\r
 334                 }
 335         if (!leadingMC)return; \r
 336         labelMC = TMath::Abs(leadingMC->GetLabel());\r
 337         }\r
 338
 339   // Trigger selection ************************************************\r
 340   if (!fAnalyseUE->TriggerSelection(fESDHandler)) return;\r
 341   \r
 342   // PILEUP-CUT ****** NEW !!!!!!!!! **************\r
 343   Bool_t select = kTRUE;\r
 344   //select = AliAnalysisHelperJetTasks::TestSelectInfo(AliAnalysisHelperJetTasks::kIsPileUp);\r
 345   if (! select) return;\r
 346   \r
 347   // Vertex selection *************************************************\r
 348   \r
 349   if(!fAnalyseUE->VertexSelection(fESDEvent, 0, fZVertex)) return;\r
 350   AliESDVertex* vertex = (AliESDVertex*)fESDEvent->GetPrimaryVertex();\r
 351   Int_t nvtx = vertex->GetNContributors();\r
 352   fhVertexContributors->Fill(nvtx);\r
 353   if (fMcHandler){\r
 354                 AliVVertex *vertexMC = (AliVVertex*)fMcEvent->GetPrimaryVertex();\r
 355                 if (vertexMC){\r
 356                         Double_t diffx = vertexMC->GetX()-vertex->GetX();\r
 357                         Double_t diffy = vertexMC->GetY()-vertex->GetY();\r
 358                         Double_t diffxy = TMath::Sqrt(TMath::Power(diffx,2)+TMath::Power(diffy,2));\r
 359                         //Double_t diffxy = TMath::Abs(diffx); // **** TMP ****\r
 360                         //Double_t diffxy = diffy;\r
 361                         Double_t diffz = TMath::Abs(vertexMC->GetZ()-vertex->GetZ());\r
 362                 \r
 363                         fhVertexReso->Fill(nvtx,diffxy,diffz);\r
 364                 }else if (fDebug>1)Printf("******* NO MC VERTEX ********");\r
 365         }
 366 \r
 367   if(!fAnalyseUE->VertexSelection(fESDEvent, fnTracksVertex, fZVertex)) return;\r
 368 \r
 369   // Get Reconstructed leading particle *******************************\r
 370   TObjArray *ltRECO = fAnalyseUE->FindLeadingObjects(fESDEvent);\r
 371   if (!ltRECO){\r
 372         delete[] ltRECO;\r
 373         return;\r
 374         }\r
 375   Int_t labelReco= TMath::Abs(((AliVParticle*)ltRECO->At(0))->GetLabel());\r
 376
 377
 378   // Loop on tracks\r
 379   Int_t nTracks = fESDEvent->GetNumberOfTracks();\r
 380   if (!nTracks)return;\r
 381   // count accepted events\r
 382   fhEntries->Fill(1.);\r
 383 \r
 384   Int_t npart=0;\r
 385   Bool_t *labelsArray = 0; \r
 386   if (fMcHandler){\r
 387         npart = fMcEvent->GetNumberOfTracks();\r
 388         labelsArray = new Bool_t[npart];\r
 389         for(Int_t j = 0; j<npart; j++){\r
 390                 labelsArray[j] = kFALSE;\r
 391                 }\r
 392         }\r
 393 \r
 394   for (Int_t i = 0; i < nTracks; i++){\r
 395         AliESDtrack *track = fESDEvent->GetTrack(i);\r
 396         // only charged\r
 397         if (!track || !track->Charge())continue;\r
 398         // apply cuts\r
 399         Bool_t cut = fEsdTrackCuts->AcceptTrack(track);\r
 400         Bool_t cutSPD = fEsdTrackCutsSPD->AcceptTrack(track);\r
 401         Bool_t cutSDD = fEsdTrackCutsSDD->AcceptTrack(track);\r
 402         Bool_t cutDCA = fEsdTrackCutsDCA->AcceptTrack(track);\r
 403         \r
 404         //Exclude the MC leading track\r
 405         Double_t matchLeading = 1.;//no match\r
 406         if (fMcHandler){\r
 407                 if (TMath::Abs(track->GetLabel())==labelMC) matchLeading=0.; //match MC leading\r
 408                 if (TMath::Abs(track->GetLabel())==labelReco) {\r
 409                         matchLeading=2.;//match RECO leading\r
 410                         if (labelMC == labelReco)matchLeading = 3.; // match both (mc = reco leading)\r
 411                         }\r
 412                 }\r
 413         // Fill step0 (all tracks) \r
 414         FillContainer(track, 0,kFALSE,matchLeading);         \r
 415         // Fill step1 (no SPD cluster requirement - no DCA cut )\r
 416         if ( cut ) FillContainer(track,1,kFALSE,matchLeading);\r
 417         // Fill step2\r
 418         if ( cut && cutDCA && (cutSPD || cutSDD)) FillContainer(track,2,kFALSE,matchLeading);\r
 419         // Fill step3-step4-step5 \r
 420         if ( cut && cutDCA && !cutSPD && !cutSDD) FillContainer(track,3,kTRUE,matchLeading);\r
 421         if ( cut && cutDCA && cutSPD) FillContainer(track,4,kTRUE,matchLeading);\r
 422         if ( cut && cutDCA && !cutSPD && cutSDD) FillContainer(track,5,kTRUE,matchLeading);\r
 423         // Fill step 6 - temporary just to define the standard track cuts  \r
 424         if ( cut && cutSPD ) FillContainer(track,6,kFALSE,matchLeading);\r
 425         //if ( cut && cutSPD ) FillContainer(track,6,kTRUE,matchLeading); // ***** TMP *****\r
 426         if ( cut && (cutSPD || cutSDD) ) FillContainer(track,7,kFALSE,matchLeading);\r
 427         // Study contamination form fakes\r
 428         if (fMcHandler){\r
 429                         \r
 430                 // consider standard ITS+TPC cuts without SPD cluster requirement\r
 431                 if (cut && cutDCA){\r
 432                         //check if it points back to any MC\r
 433                         Int_t label = TMath::Abs(track->GetLabel());\r
 434                         AliVParticle *part = (AliVParticle*)fMcEvent->GetTrack(label);  \r
 435                         if (!part){\r
 436                                 fhFakes->Fill(0.);\r
 437                                 //Printf("*************** NO MC PARTICLE ************************");\r
 438                                 continue;\r
 439                                 }\r
 440                                 \r
 441                         fhPtMCAll->Fill(part->Pt()); \r
 442                         // Check if label is not already in array\r
 443                         if (!labelsArray[label]){\r
 444                                 labelsArray[label]= kTRUE;\r
 445                                 if (fMcEvent->IsPhysicalPrimary(label)){\r
 446                                         fhFakes->Fill(1.);\r
 447                                         fhPtMCPrim->Fill(part->Pt());\r
 448                                 }else{\r
 449                                         fhFakes->Fill(2.);\r
 450                                         fhPtMCSec->Fill(part->Pt());\r
 451                                         }\r
 452                         }else{\r
 453                                 fhFakes->Fill(3.);\r
 454                                 if (fMcEvent->IsPhysicalPrimary(label))fhPtMCPrimFake->Fill(part->Pt());\r
 455                                 else fhPtMCSecFake->Fill(part->Pt());\r
 456                                 }\r
 457                    }    \r
 458                 }\r
 459         } // end loop on tracks\r
 460         if(labelsArray)\r
 461         delete[] labelsArray;\r
 462 \r
 463 }\r
 464 \r
 465 \r
 466 //____________________________________________________________________\r
 467 void AliAnalysisTaskCorrectionsUE::CreateContainer()\r
 468 {\r
 469
 470   // Create the output CF container\r
 471   // relevant variables \r
 472   UInt_t iprim  = 0; // 0: primaries, 1: secondaries from strangness, 2: secondaries form material \r
 473   UInt_t iptmc  = 1;\r
 474   UInt_t ipt    = 2;\r
 475   UInt_t ieta   = 3;\r
 476   UInt_t idcaxy = 4;\r
 477   UInt_t idcaz  = 5;\r
 478   UInt_t imatch = 6;\r
 479   UInt_t icharge = 7;\r
 480   // set-up the grid\r
 481   UInt_t nstep = 8;\r
 482   const Int_t nvar  = 8;\r
 483   const Int_t nbin0 = 5;  // prim\r
 484   const Int_t nbin1 = 20; // pt resolution\r
 485   const Int_t nbin2 = 39; // pt \r
 486   const Int_t nbin3 = 20; // eta\r
 487   const Int_t nbin4 = 100; // dca xy\r
 488   const Int_t nbin5 = 100; // dca z\r
 489   const Int_t nbin6 = 4; // matching with leading track\r
 490   const Int_t nbin7 = 2;\r
 491 \r
 492   // array for the number of bins in each dimension\r
 493   Int_t iBin[nvar];\r
 494   iBin[0] = nbin0;\r
 495   iBin[1] = nbin1;\r
 496   iBin[2] = nbin2;\r
 497   iBin[3] = nbin3;\r
 498   iBin[4] = nbin4;\r
 499   iBin[5] = nbin5;\r
 500   iBin[6] = nbin6;\r
 501   iBin[7] = nbin7;\r
 502
 503   // primaries\r
 504   Double_t primBins[7] = {-0.5,0.5,1.5,2.5,3.5,4.5,5.5};\r
 505   // matching with leading-track\r
 506   Double_t matchBins[5] = {-0.5,0.5,1.5,2.5,3.5};\r
 507
 508   // pT resolution\r
 509   Double_t resoBins[nbin1+1];\r
 510   for (Int_t i=0; i<=nbin1; i++)\r
 511     resoBins[i] = -1.0 + 0.1 * i;\r
 512
 513   // pT\r
 514   Double_t ptBins[nbin2+1] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};\r
 515
 516   // eta\r
 517   Double_t etaBins[nbin3+1];\r
 518   for (Int_t i=0; i<=nbin3; i++)\r
 519     etaBins[i] = -1.0 + 0.1 * i;\r
 520
 521   // dca xy\r
 522   Double_t dcaxyBins[nbin4+1];\r
 523   for (Int_t i=0; i<=nbin4; i++)\r
 524     dcaxyBins[i] = -1.+0.02 * i;\r
 525
 526   // dca z\r
 527   Double_t dcazBins[nbin5+1];\r
 528   for (Int_t i=0; i<=nbin5; i++)\r
 529     dcazBins[i] = -5.0 + 0.1 * i;\r
 530  \r
 531   // charge \r
 532   Double_t chargeBins[nbin7+1] = {-1.,0.,1.};\r
 533
 534   // create container\r
 535   // set variables\r
 536   fOutCFcont = new AliCFContainer("fOutCFcont","Output Container",nstep,nvar,iBin);\r
 537   fOutCFcont->SetBinLimits(iprim,primBins);\r
 538   fOutCFcont->SetVarTitle(iprim, "Particle type");\r
 539   fOutCFcont->SetBinLimits(iptmc,resoBins);\r
 540   fOutCFcont->SetVarTitle(iptmc, "#Delta p_{T} (DATA-MC) (GeV/c)");\r
 541   fOutCFcont->SetBinLimits(ipt,ptBins);\r
 542   fOutCFcont->SetVarTitle(ipt, "p_{T} (GeV/c)");\r
 543   fOutCFcont->SetBinLimits(ieta,etaBins);\r
 544   fOutCFcont->SetVarTitle(ieta, "#eta");\r
 545   fOutCFcont->SetBinLimits(idcaxy,dcaxyBins);\r
 546   fOutCFcont->SetVarTitle(idcaxy, " DCA_{XY} (cm)");\r
 547   fOutCFcont->SetBinLimits(idcaz,dcazBins);\r
 548   fOutCFcont->SetVarTitle(idcaz, " DCA_{Z} (cm)");\r
 549   fOutCFcont->SetBinLimits(imatch,matchBins);\r
 550   fOutCFcont->SetVarTitle(imatch, "Matching with leading-track");\r
 551   fOutCFcont->SetBinLimits(icharge,chargeBins);\r
 552   fOutCFcont->SetVarTitle(icharge, "Charge");\r
 553 \r
 554   // set steps\r
 555   fOutCFcont->SetStepTitle(0,"all tracks");\r
 556   fOutCFcont->SetStepTitle(1,"ITS+TPC cuts (no SPD cluster requirement and DCA cut)");\r
 557   fOutCFcont->SetStepTitle(2,"add DCA cut");\r
 558   fOutCFcont->SetStepTitle(3,"NO SPD cluster, NO SDD cluster in first layer");\r
 559   fOutCFcont->SetStepTitle(4,"YES SPD cluster, NO SDD cluster in first layer");\r
 560   fOutCFcont->SetStepTitle(5,"NO SPD cluster, YES SDD cluster in first layer");\r
 561   fOutCFcont->SetStepTitle(6,"ITS+TPC cuts - no DCA cut - SPD cut");\r
 562   fOutCFcont->SetStepTitle(7,"ITS+TPC cuts - no DCA cut - SPD or SDD cut");\r
 563   fListOfHistos->Add(fOutCFcont);\r
 564 \r
 565 }\r
 566 \r
 567 \r
 568 void  AliAnalysisTaskCorrectionsUE::FillContainer(AliESDtrack *track, Int_t step,Bool_t mcvertex, Double_t matchLeading)\r
 569 {\r
 570 \r
 571   // Fill the CF container\r
 572 \r
 573   Double_t vars[8];\r
 574   Double_t prim = -1.;\r
 575
 576   if (track->Charge() > 0.) vars[7] = 0.5;\r
 577   else vars[7] = -0.5;\r
 578
 579   if (fMcHandler){\r
 580         // determine if points back to a primary\r
 581         Int_t label = TMath::Abs(track->GetLabel());\r
 582         \r
 583         AliMCParticle *part = (AliMCParticle*)fMcEvent->GetTrack(label);  \r
 584         for (Int_t i=0; i<=6;i++) vars[i] = -999.;\r
 585         if (part) { //PRIMARY\r
 586                 if (fMcEvent->IsPhysicalPrimary(label)){\r
 587                         prim = 0.;\r
 588                 }else { //SECONDARY\r
 589                         // decide if strange\r
 590                         Int_t labelm = TMath::Abs(part->GetMother());\r
 591                         AliMCParticle *mother = (AliMCParticle*)fMcEvent->GetTrack(labelm);  \r
 592                         Int_t code = mother->PdgCode();\r
 593                         Int_t mfl = Int_t (code/ TMath::Power(10, Int_t(TMath::Log10(code))));\r
 594                          if  (mfl == 3) prim = 1.;\r
 595                                 else{   \r
 596                                         //Printf("***** PROCESS : %d",part->Particle()->GetUniqueID());  \r
 597                                         if (TMath::Abs(code) == 211) prim = 2.; // charged pion decay\r
 598                                         else if (part->Particle()->GetUniqueID() == 13 )prim = 3.; // hadronic interactions\r
 599                                         else if (part->Particle()->GetUniqueID() == 5 )prim = 4.; // photon conversions\r
 600                                         else prim = 5.; // other?\r
 601                                 }\r
 602                         }\r
 603                 vars[1]= part->Pt()-track->Pt();\r
 604                 // In step 2 fill MC pT for contamination study\r
 605                 if (step == 2)vars[2]=part->Pt();\r
 606                 }
 607         }\r
 608   vars[0]=prim;\r
 609
 610   if (step != 2)vars[2]=track->Pt();\r
 611   vars[3]=track->Eta();\r
 612
 613   Bool_t dcaControlFlag = kFALSE;\r
 614   if (mcvertex && fMcHandler){\r
 615           // we want DCA w.r.t. MC vertex\r
 616           AliVVertex *vtxMC = (AliVVertex*)fMcEvent->GetPrimaryVertex();\r
 617           dcaControlFlag = track->RelateToVertex((AliESDVertex*)vtxMC,(Double_t)fESDEvent->GetMagneticField(),10000.);\r
 618           }else{\r
 619           AliESDVertex* vertex = (AliESDVertex*)fESDEvent->GetPrimaryVertex();\r
 620           dcaControlFlag = track->RelateToVertex(vertex,(Double_t)fESDEvent->GetMagneticField(),10000.);\r
 621           }\r
 622   if (dcaControlFlag){\r
 623         Float_t dca[2];\r
 624         Float_t dcaCov[2];\r
 625         track->GetImpactParameters(dca,dcaCov);\r
 626         vars[4]=dca[0];\r
 627         vars[5]=dca[1];\r
 628         }         \r
 629
 630
 631   vars[6]= matchLeading;\r
 632   fOutCFcont->Fill(vars,step);\r
 633 \r
 634 }
 635
 636 \r
 637 //____________________________________________________________________
 638 AliAnalysisTaskCorrectionsUE* AliAnalysisTaskCorrectionsUE::Instance()
 639 {
 640
 641   //Create instance
 642   if (fgTaskCorrectionsUE) {
 643         return fgTaskCorrectionsUE;
 644   } else {
 645         fgTaskCorrectionsUE = new AliAnalysisTaskCorrectionsUE();
 646         return fgTaskCorrectionsUE;
 647         }
 648 }
 649