ANALYSIS/AliAnalysisTaskESDfilter.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 #include <TChain.h>
  19 #include <TTree.h>
  20 #include <TList.h>
  21 #include <TArrayI.h>
  22 #include <TParameter.h>
  23 #include <TRandom.h>
  24 #include <TParticle.h>
  25 #include <TFile.h>
  26
  27 #include "AliAnalysisTaskESDfilter.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliESDEvent.h"
  30 #include "AliESDRun.h"
  31 #include "AliStack.h"
  32 #include "AliAODEvent.h"
  33 #include "AliMCEvent.h"
  34 #include "AliMCEventHandler.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAODMCParticle.h"
  38 #include "AliAnalysisFilter.h"
  39 #include "AliESDMuonTrack.h"
  40 #include "AliESDVertex.h"
  41 #include "AliCentrality.h"
  42 #include "AliEventplane.h"
  43 #include "AliESDv0.h"
  44 #include "AliESDkink.h"
  45 #include "AliESDcascade.h"
  46 #include "AliESDPmdTrack.h"
  47 #include "AliESDCaloCluster.h"
  48 #include "AliESDCaloCells.h"
  49 #include "AliMultiplicity.h"
  50 #include "AliLog.h"
  51 #include "AliCodeTimer.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "AliESDpid.h"
  54 #include "AliAODHMPIDrings.h"
  55 #include "AliV0vertexer.h"
  56 #include "AliCascadeVertexer.h"
  57 #include "Riostream.h"
  58 #include "AliExternalTrackParam.h"
  59 #include "AliTrackerBase.h"
  60 #include "TVector3.h"
  61 #include "AliTPCdEdxInfo.h"
  62
  63 using std::cout;
  64 using std::endl;
  65 ClassImp(AliAnalysisTaskESDfilter)
  66
  67 ////////////////////////////////////////////////////////////////////////
  68
  69 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():
  70   AliAnalysisTaskSE(),
  71   fTrackFilter(0x0),
  72   fKinkFilter(0x0),
  73   fV0Filter(0x0),
  74   fCascadeFilter(0x0),
  75   fHighPthreshold(0),
  76   fPtshape(0x0),
  77   fEnableFillAOD(kTRUE),
  78   fUsedTrack(0x0),
  79   fUsedKink(0x0),
  80   fUsedV0(0x0),
  81   fAODTrackRefs(0x0),
  82   fAODV0VtxRefs(0x0),
  83   fAODV0Refs(0x0),
  84   fMChandler(0x0),
  85   fNumberOfTracks(0),
  86   fNumberOfPositiveTracks(0),
  87   fNumberOfV0s(0),
  88   fNumberOfVertices(0),
  89   fNumberOfCascades(0),
  90   fNumberOfKinks(0),
  91   fOldESDformat(kFALSE),
  92   fPrimaryVertex(0x0),
  93   fTPCConstrainedFilterMask(0),
  94   fHybridFilterMaskTPCCG(0),
  95   fWriteHybridTPCCOnly(kFALSE),
  96   fGlobalConstrainedFilterMask(0),
  97   fHybridFilterMaskGCG(0),
  98   fWriteHybridGCOnly(kFALSE),
  99   fIsVZEROEnabled(kTRUE),
 100   fIsTZEROEnabled(kTRUE),
 101   fIsZDCEnabled(kTRUE),
 102   fIsHMPIDEnabled(kTRUE),
 103   fIsV0CascadeRecoEnabled(kFALSE),
 104   fAreCascadesEnabled(kTRUE),
 105   fAreV0sEnabled(kTRUE),
 106   fAreKinksEnabled(kTRUE),
 107   fAreTracksEnabled(kTRUE),
 108   fArePmdClustersEnabled(kTRUE),
 109   fAreCaloClustersEnabled(kTRUE),
 110   fAreEMCALCellsEnabled(kTRUE),
 111   fArePHOSCellsEnabled(kTRUE),
 112   fAreEMCALTriggerEnabled(kTRUE),
 113   fArePHOSTriggerEnabled(kTRUE),
 114   fAreTrackletsEnabled(kTRUE),
 115   fESDpid(0x0),
 116   fIsPidOwner(kFALSE),
 117   fTPCaloneTrackCuts(0),
 118   fDoPropagateTrackToEMCal(kTRUE)
 119 {
 120   // Default constructor
 121     fV0Cuts[0] =  33.   ;   // max allowed chi2
 122     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 123     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 124     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 125     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 126     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 127     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 128
 129     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 130     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 131     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 132     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 133     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 134     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 135     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 136     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 137 }
 138
 139 //______________________________________________________________________________
 140 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):
 141     AliAnalysisTaskSE(name),
 142     fTrackFilter(0x0),
 143     fKinkFilter(0x0),
 144     fV0Filter(0x0),
 145     fCascadeFilter(0x0),
 146     fHighPthreshold(0),
 147     fPtshape(0x0),
 148     fEnableFillAOD(kTRUE),
 149     fUsedTrack(0x0),
 150     fUsedKink(0x0),
 151     fUsedV0(0x0),
 152     fAODTrackRefs(0x0),
 153     fAODV0VtxRefs(0x0),
 154     fAODV0Refs(0x0),
 155     fMChandler(0x0),
 156     fNumberOfTracks(0),
 157     fNumberOfPositiveTracks(0),
 158     fNumberOfV0s(0),
 159     fNumberOfVertices(0),
 160     fNumberOfCascades(0),
 161     fNumberOfKinks(0),
 162     fOldESDformat(kFALSE),
 163     fPrimaryVertex(0x0),
 164   fTPCConstrainedFilterMask(0),
 165   fHybridFilterMaskTPCCG(0),
 166   fWriteHybridTPCCOnly(kFALSE),
 167   fGlobalConstrainedFilterMask(0),
 168   fHybridFilterMaskGCG(0),
 169   fWriteHybridGCOnly(kFALSE),
 170     fIsVZEROEnabled(kTRUE),
 171     fIsTZEROEnabled(kTRUE),
 172     fIsZDCEnabled(kTRUE),
 173     fIsHMPIDEnabled(kTRUE),
 174     fIsV0CascadeRecoEnabled(kFALSE),
 175     fAreCascadesEnabled(kTRUE),
 176     fAreV0sEnabled(kTRUE),
 177     fAreKinksEnabled(kTRUE),
 178     fAreTracksEnabled(kTRUE),
 179     fArePmdClustersEnabled(kTRUE),
 180     fAreCaloClustersEnabled(kTRUE),
 181     fAreEMCALCellsEnabled(kTRUE),
 182     fArePHOSCellsEnabled(kTRUE),
 183                 fAreEMCALTriggerEnabled(kTRUE),
 184                 fArePHOSTriggerEnabled(kTRUE),
 185                 fAreTrackletsEnabled(kTRUE),
 186     fESDpid(0x0),
 187     fIsPidOwner(kFALSE),
 188     fTPCaloneTrackCuts(0),
 189   fDoPropagateTrackToEMCal(kTRUE)
 190 {
 191   // Constructor
 192
 193     fV0Cuts[0] =  33.   ;   // max allowed chi2
 194     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
 195     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
 196     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
 197     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
 198     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
 199     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
 200
 201     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
 202     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
 203     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
 204     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
 205     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
 206     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
 207     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
 208     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
 209
 210
 211
 212 }
 213 AliAnalysisTaskESDfilter::~AliAnalysisTaskESDfilter(){
 214     if(fIsPidOwner) delete fESDpid;
 215 }
 216 //______________________________________________________________________________
 217 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()
 218 {
 219   //
 220   // Create Output Objects conenct filter to outputtree
 221   //
 222   if(OutputTree())
 223   {
 224     OutputTree()->GetUserInfo()->Add(fTrackFilter);
 225   }
 226   else
 227   {
 228     AliError("No OutputTree() for adding the track filter");
 229   }
 230   fTPCaloneTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 231 }
 232
 233 //______________________________________________________________________________
 234 void AliAnalysisTaskESDfilter::Init()
 235 {
 236   // Initialization
 237   if (fDebug > 1) AliInfo("Init() \n");
 238   // Call configuration file
 239 }
 240
 241 //______________________________________________________________________________
 242 Bool_t AliAnalysisTaskESDfilter::Notify()
 243 {
 244 // Notify method.
 245    AddMetadataToUserInfo();
 246    return kTRUE;
 247 }
 248
 249 //______________________________________________________________________________
 250 Bool_t AliAnalysisTaskESDfilter::AddMetadataToUserInfo()
 251 {
 252 // Copy metadata to AOD user info.
 253   static Bool_t copyFirst = kFALSE;
 254   if (!copyFirst) {
 255     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 256     if (!mgr) {
 257        AliError("AliAnalysisTaskESDfilter::AddMetadataToUserInfo() : No analysis manager !");
 258        return kFALSE;
 259     }
 260     TTree *esdTree = mgr->GetTree()->GetTree();
 261     if (!esdTree) return kFALSE;
 262     TNamed *alirootVersion = (TNamed*)esdTree->GetUserInfo()->FindObject("alirootVersion");
 263     if (!alirootVersion) return kFALSE;
 264     AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
 265          if (!aodHandler) return kFALSE;
 266     TTree *aodTree = aodHandler->GetTree();
 267     if (!aodTree) return kFALSE;
 268     aodTree->GetUserInfo()->Add(new TNamed(*alirootVersion));
 269     copyFirst = kTRUE;
 270   }
 271   return kTRUE;
 272 }
 273
 274 //______________________________________________________________________________
 275 void AliAnalysisTaskESDfilter::PrintTask(Option_t *option, Int_t indent) const
 276 {
 277 // Print selection task information
 278   AliInfo("");
 279
 280   AliAnalysisTaskSE::PrintTask(option,indent);
 281
 282   TString spaces(' ',indent+3);
 283
 284         cout << spaces.Data() << Form("Cascades       are %s",fAreCascadesEnabled ? "ENABLED":"DISABLED") << endl;
 285         cout << spaces.Data() << Form("V0s            are %s",fAreV0sEnabled ? "ENABLED":"DISABLED") << endl;
 286         cout << spaces.Data() << Form("Kinks          are %s",fAreKinksEnabled ? "ENABLED":"DISABLED") << endl;
 287         cout << spaces.Data() << Form("Tracks         are %s",fAreTracksEnabled ? "ENABLED":"DISABLED") << endl;
 288         cout << spaces.Data() << Form("PmdClusters    are %s",fArePmdClustersEnabled ? "ENABLED":"DISABLED") << endl;
 289         cout << spaces.Data() << Form("CaloClusters   are %s",fAreCaloClustersEnabled ? "ENABLED":"DISABLED") << endl;
 290   cout << spaces.Data() << Form("EMCAL cells    are %s",fAreEMCALCellsEnabled ? "ENABLED":"DISABLED") << endl;
 291         cout << spaces.Data() << Form("EMCAL triggers are %s",fAreEMCALTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 292         cout << spaces.Data() << Form("PHOS triggers  are %s",fArePHOSTriggerEnabled ? "ENABLED":"DISABLED") << endl;
 293         cout << spaces.Data() << Form("Tracklets      are %s",fAreTrackletsEnabled ? "ENABLED":"DISABLED") << endl;
 294         cout << spaces.Data() << Form("PropagateTrackToEMCal  is %s", fDoPropagateTrackToEMCal ? "ENABLED":"DISABLED") << endl;
 295 }
 296
 297 //______________________________________________________________________________
 298 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)
 299 {
 300 // Execute analysis for current event
 301 //
 302
 303   Long64_t ientry = Entry();
 304
 305   if (fDebug > 0) {
 306       printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
 307       if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
 308       if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
 309   }
 310   // Filters must explicitely enable AOD filling in their UserExec (AG)
 311   if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
 312   if(fEnableFillAOD) {
 313      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 314      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
 315   }
 316   ConvertESDtoAOD();
 317 }
 318
 319 //______________________________________________________________________________
 320 TClonesArray& AliAnalysisTaskESDfilter::Cascades()
 321 {
 322   return *(AODEvent()->GetCascades());
 323 }
 324
 325 //______________________________________________________________________________
 326 TClonesArray& AliAnalysisTaskESDfilter::Tracks()
 327 {
 328   return *(AODEvent()->GetTracks());
 329 }
 330
 331 //______________________________________________________________________________
 332 TClonesArray& AliAnalysisTaskESDfilter::V0s()
 333 {
 334   return *(AODEvent()->GetV0s());
 335 }
 336
 337 //______________________________________________________________________________
 338 TClonesArray& AliAnalysisTaskESDfilter::Vertices()
 339 {
 340   return *(AODEvent()->GetVertices());
 341 }
 342
 343 //______________________________________________________________________________
 344 AliAODHeader* AliAnalysisTaskESDfilter::ConvertHeader(const AliESDEvent& esd)
 345 {
 346 // Convert header information
 347
 348   AliCodeTimerAuto("",0);
 349
 350   AliAODHeader* header = AODEvent()->GetHeader();
 351
 352   header->SetRunNumber(esd.GetRunNumber());
 353   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 354
 355   TTree* tree = fInputHandler->GetTree();
 356   if (tree) {
 357     TFile* file = tree->GetCurrentFile();
 358     if (file) header->SetESDFileName(file->GetName());
 359   }
 360
 361   if (fOldESDformat) {
 362     header->SetBunchCrossNumber(0);
 363     header->SetOrbitNumber(0);
 364     header->SetPeriodNumber(0);
 365     header->SetEventType(0);
 366     header->SetMuonMagFieldScale(-999.);
 367     header->SetCentrality(0);
 368     header->SetEventplane(0);
 369   } else {
 370     header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
 371     header->SetOrbitNumber(esd.GetOrbitNumber());
 372     header->SetPeriodNumber(esd.GetPeriodNumber());
 373     header->SetEventType(esd.GetEventType());
 374
 375     header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
 376     if(const_cast<AliESDEvent&>(esd).GetCentrality()){
 377       header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
 378     }
 379     else{
 380       header->SetCentrality(0);
 381     }
 382     if(const_cast<AliESDEvent&>(esd).GetEventplane()){
 383       header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
 384     }
 385     else{
 386       header->SetEventplane(0);
 387     }
 388   }
 389
 390   // Trigger
 391   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
 392   header->SetTriggerMask(esd.GetTriggerMask());
 393   header->SetTriggerCluster(esd.GetTriggerCluster());
 394   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());
 395   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());
 396   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());
 397
 398   header->SetMagneticField(esd.GetMagneticField());
 399   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
 400   header->SetZDCN1Energy(esd.GetZDCN1Energy());
 401   header->SetZDCP1Energy(esd.GetZDCP1Energy());
 402   header->SetZDCN2Energy(esd.GetZDCN2Energy());
 403   header->SetZDCP2Energy(esd.GetZDCP2Energy());
 404   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
 405
 406   // ITS Cluster Multiplicty
 407   const AliMultiplicity *mult = esd.GetMultiplicity();
 408   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
 409
 410   // TPC only Reference Multiplicty
 411   Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
 412   header->SetTPConlyRefMultiplicity(refMult);
 413   //
 414   header->SetRefMultiplicityComb05(AliESDtrackCuts::GetReferenceMultiplicity(&esd,AliESDtrackCuts::kTrackletsITSTPC,0.5));
 415   header->SetRefMultiplicityComb08(AliESDtrackCuts::GetReferenceMultiplicity(&esd,AliESDtrackCuts::kTrackletsITSTPC,0.8));
 416   //
 417   Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};
 418   Float_t diamcov[3];
 419   esd.GetDiamondCovXY(diamcov);
 420   header->SetDiamond(diamxy,diamcov);
 421   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
 422
 423   // VZERO channel equalization factors for event-plane reconstruction
 424   header->SetVZEROEqFactors(esd.GetVZEROEqFactors());
 425
 426   // T0 Resolution information
 427   const AliESDRun* esdRun = esd.GetESDRun();
 428   for (Int_t i=0;i<AliESDRun::kT0spreadSize;i++) header->SetT0spread(i,esdRun->GetT0spread(i));
 429
 430   return header;
 431 }
 432
 433 //______________________________________________________________________________
 434 void AliAnalysisTaskESDfilter::ConvertCascades(const AliESDEvent& esd)
 435 {
 436
 437   // Convert the cascades part of the ESD.
 438   // Return the number of cascades
 439
 440   AliCodeTimerAuto("",0);
 441
 442   // Create vertices starting from the most complex objects
 443   Double_t chi2 = 0.;
 444
 445   const AliESDVertex* vtx = esd.GetPrimaryVertex();
 446   Double_t pos[3] = { 0. };
 447   Double_t covVtx[6] = { 0. };
 448   Double_t momBach[3]={0.};
 449   Double_t covTr[21]={0.};
 450   Double_t pid[10]={0.};
 451   AliAODPid* detpid(0x0);
 452   AliAODVertex* vV0FromCascade(0x0);
 453   AliAODv0* aodV0(0x0);
 454   AliAODcascade* aodCascade(0x0);
 455   AliAODTrack* aodTrack(0x0);
 456   Double_t momPos[3]={0.};
 457   Double_t momNeg[3] = { 0. };
 458   Double_t momPosAtV0vtx[3]={0.};
 459   Double_t momNegAtV0vtx[3]={0.};
 460
 461   TClonesArray& verticesArray = Vertices();
 462   TClonesArray& tracksArray = Tracks();
 463   TClonesArray& cascadesArray = Cascades();
 464
 465   // Cascades (Modified by A.Maire - February 2009)
 466   for (Int_t nCascade = 0; nCascade < esd.GetNumberOfCascades(); ++nCascade) {
 467
 468     // 0- Preparation
 469     //
 470     AliESDcascade *esdCascade = esd.GetCascade(nCascade);
 471                 Int_t  idxPosFromV0Dghter  = esdCascade->GetPindex();
 472                 Int_t  idxNegFromV0Dghter  = esdCascade->GetNindex();
 473                 Int_t  idxBachFromCascade  = esdCascade->GetBindex();
 474
 475     AliESDtrack  *esdCascadePos  = esd.GetTrack( idxPosFromV0Dghter);
 476     AliESDtrack  *esdCascadeNeg  = esd.GetTrack( idxNegFromV0Dghter);
 477     AliESDtrack  *esdCascadeBach = esd.GetTrack( idxBachFromCascade);
 478
 479     // Identification of the V0 within the esdCascade (via both daughter track indices)
 480     AliESDv0 * currentV0   = 0x0;
 481     Int_t      idxV0FromCascade = -1;
 482
 483     for (Int_t iV0=0; iV0<esd.GetNumberOfV0s(); ++iV0) {
 484
 485       currentV0 = esd.GetV0(iV0);
 486       Int_t posCurrentV0 = currentV0->GetPindex();
 487       Int_t negCurrentV0 = currentV0->GetNindex();
 488
 489       if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {
 490         idxV0FromCascade = iV0;
 491         break;
 492       }
 493     }
 494
 495     if(idxV0FromCascade < 0){
 496       printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");
 497       continue;
 498     }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"
 499
 500     AliESDv0 *esdV0FromCascade   = esd.GetV0(idxV0FromCascade);
 501
 502     // 1 - Cascade selection
 503
 504     //  AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd.GetPrimaryVertex()));
 505     //  TList cascadeObjects;
 506     //  cascadeObjects.AddAt(esdV0FromCascade, 0);
 507     //  cascadeObjects.AddAt(esdCascadePos,    1);
 508     //  cascadeObjects.AddAt(esdCascadeNeg,    2);
 509     //  cascadeObjects.AddAt(esdCascade,       3);
 510     //  cascadeObjects.AddAt(esdCascadeBach,   4);
 511     //  cascadeObjects.AddAt(esdPrimVtx,       5);
 512     //
 513     //  UInt_t selectCascade = 0;
 514     //  if (fCascadeFilter) {
 515     //    // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects);
 516     //          // FIXME AliESDCascadeCuts to be implemented ...
 517     //
 518     //          // Here we may encounter a moot point at the V0 level
 519     //          // between the cascade selections and the V0 ones :
 520     //          // the V0 selected along with the cascade (secondary V0) may
 521     //          // usually be removed from the dedicated V0 selections (prim V0) ...
 522     //          // -> To be discussed !
 523     //
 524     //    // this is a little awkward but otherwise the
 525     //    // list wants to access the pointer (delete it)
 526     //    // again when going out of scope
 527     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 528     //    esdPrimVtx = 0;
 529     //    if (!selectCascade)
 530     //      continue;
 531     //  }
 532     //  else{
 533     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
 534     //    esdPrimVtx = 0;
 535     //  }
 536
 537     // 2 - Add the cascade vertex
 538
 539     esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);
 540     esdCascade->GetPosCovXi(covVtx);
 541     chi2 = esdCascade->GetChi2Xi();
 542
 543     AliAODVertex *vCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex( pos,
 544                                                                      covVtx,
 545                                                                      chi2, // FIXME = Chi2/NDF will be needed
 546                                                                      fPrimaryVertex,
 547                                                                      nCascade, // id
 548                                                                      AliAODVertex::kCascade);
 549     fPrimaryVertex->AddDaughter(vCascade);
 550
 551 //    if (fDebug > 2) {
 552 //      printf("---- Cascade / Cascade Vertex (AOD) : \n");
 553 //      vCascade->Print();
 554 //    }
 555
 556 //    if(esd.GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(const_cast<AliESDEvent*>(&esd), (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production starting form LHC10e without Tender.
 557
 558
 559     // 3 - Add the bachelor track from the cascade
 560
 561     if (!fUsedTrack[idxBachFromCascade]) {
 562
 563       esdCascadeBach->GetPxPyPz(momBach);
 564       esdCascadeBach->GetXYZ(pos);
 565       esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);
 566       esdCascadeBach->GetESDpid(pid);
 567
 568             fUsedTrack[idxBachFromCascade] = kTRUE;
 569             UInt_t selectInfo = 0;
 570             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);
 571             if (fMChandler) fMChandler->SelectParticle(esdCascadeBach->GetLabel());
 572             aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(esdCascadeBach->GetID(),
 573                                                                            esdCascadeBach->GetLabel(),
 574                                                                            momBach,
 575                                                                            kTRUE,
 576                                                                            pos,
 577                                                                            kFALSE, // Why kFALSE for "isDCA" ? FIXME
 578                                                                            covTr,
 579                                                                            (Short_t)esdCascadeBach->GetSign(),
 580                                                                            esdCascadeBach->GetITSClusterMap(),
 581                                                                            pid,
 582                                                                            vCascade,
 583                                                                            kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 584                                                                            vtx->UsesTrack(esdCascadeBach->GetID()),
 585                                                                            AliAODTrack::kSecondary,
 586                                                                            selectInfo);
 587             aodTrack->SetTPCFitMap(esdCascadeBach->GetTPCFitMap());
 588             aodTrack->SetTPCClusterMap(esdCascadeBach->GetTPCClusterMap());
 589             aodTrack->SetTPCSharedMap (esdCascadeBach->GetTPCSharedMap());
 590             aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeBach));
 591             aodTrack->SetTPCPointsF(esdCascadeBach->GetTPCNclsF());
 592             aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeBach->GetTPCCrossedRows()));
 593             fAODTrackRefs->AddAt(aodTrack,idxBachFromCascade);
 594
 595             if (esdCascadeBach->GetSign() > 0) ++fNumberOfPositiveTracks;
 596             aodTrack->ConvertAliPIDtoAODPID();
 597             aodTrack->SetFlags(esdCascadeBach->GetStatus());
 598       SetAODPID(esdCascadeBach,aodTrack,detpid);
 599     }
 600     else {
 601             aodTrack = static_cast<AliAODTrack*>( fAODTrackRefs->At(idxBachFromCascade) );
 602     }
 603
 604     vCascade->AddDaughter(aodTrack);
 605
 606 //    if (fDebug > 4) {
 607 //      printf("---- Cascade / bach dghter : \n");
 608 //      aodTrack->Print();
 609 //    }
 610
 611
 612     // 4 - Add the V0 from the cascade.
 613     // = V0vtx + both pos and neg daughter tracks + the aodV0 itself
 614     //
 615
 616     if ( !fUsedV0[idxV0FromCascade] ) {
 617       // 4.A - if VO structure hasn't been created yet
 618
 619       // 4.A.1 - Create the V0 vertex of the cascade
 620
 621       esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);
 622       esdV0FromCascade->GetPosCov(covVtx);
 623       chi2 = esdV0FromCascade->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3 ?
 624
 625       vV0FromCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex(pos,
 626                                                                          covVtx,
 627                                                                          chi2,
 628                                                                          vCascade,
 629                                                                          idxV0FromCascade, //id of ESDv0
 630                                                                          AliAODVertex::kV0);
 631       // Note:
 632       //    one V0 can be used by several cascades.
 633       // So, one AOD V0 vtx can have several parent vtx.
 634       // This is not directly allowed by AliAODvertex.
 635       // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash
 636       // but to a problem of consistency within AODEvent.
 637       // -> See below paragraph 4.B, for the proposed treatment of such a case.
 638
 639       // Add the vV0FromCascade to the aodVOVtxRefs
 640       fAODV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);
 641
 642
 643       // 4.A.2 - Add the positive tracks from the V0
 644
 645       esdCascadePos->GetPxPyPz(momPos);
 646       esdCascadePos->GetXYZ(pos);
 647       esdCascadePos->GetCovarianceXYZPxPyPz(covTr);
 648       esdCascadePos->GetESDpid(pid);
 649
 650
 651       if (!fUsedTrack[idxPosFromV0Dghter]) {
 652         fUsedTrack[idxPosFromV0Dghter] = kTRUE;
 653
 654         UInt_t selectInfo = 0;
 655         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);
 656         if(fMChandler) fMChandler->SelectParticle(esdCascadePos->GetLabel());
 657         aodTrack = new(tracksArray[fNumberOfTracks++])
 658         AliAODTrack(  esdCascadePos->GetID(),
 659                     esdCascadePos->GetLabel(),
 660                     momPos,
 661                     kTRUE,
 662                     pos,
 663                     kFALSE, // Why kFALSE for "isDCA" ? FIXME
 664                     covTr,
 665                     (Short_t)esdCascadePos->GetSign(),
 666                     esdCascadePos->GetITSClusterMap(),
 667                     pid,
 668                     vV0FromCascade,
 669                     kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 670                     vtx->UsesTrack(esdCascadePos->GetID()),
 671                     AliAODTrack::kSecondary,
 672                     selectInfo);
 673         aodTrack->SetTPCFitMap(esdCascadePos->GetTPCFitMap());
 674         aodTrack->SetTPCClusterMap(esdCascadePos->GetTPCClusterMap());
 675         aodTrack->SetTPCSharedMap (esdCascadePos->GetTPCSharedMap());
 676         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadePos));
 677         aodTrack->SetTPCPointsF(esdCascadePos->GetTPCNclsF());
 678         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadePos->GetTPCCrossedRows()));
 679         fAODTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);
 680
 681         if (esdCascadePos->GetSign() > 0) ++fNumberOfPositiveTracks;
 682         aodTrack->ConvertAliPIDtoAODPID();
 683         aodTrack->SetFlags(esdCascadePos->GetStatus());
 684         SetAODPID(esdCascadePos,aodTrack,detpid);
 685       }
 686       else {
 687         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxPosFromV0Dghter));
 688       }
 689       vV0FromCascade->AddDaughter(aodTrack);
 690
 691
 692       // 4.A.3 - Add the negative tracks from the V0
 693
 694       esdCascadeNeg->GetPxPyPz(momNeg);
 695       esdCascadeNeg->GetXYZ(pos);
 696       esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);
 697       esdCascadeNeg->GetESDpid(pid);
 698
 699
 700       if (!fUsedTrack[idxNegFromV0Dghter]) {
 701         fUsedTrack[idxNegFromV0Dghter] = kTRUE;
 702
 703         UInt_t selectInfo = 0;
 704         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);
 705         if(fMChandler)fMChandler->SelectParticle(esdCascadeNeg->GetLabel());
 706         aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(  esdCascadeNeg->GetID(),
 707                                                       esdCascadeNeg->GetLabel(),
 708                                                       momNeg,
 709                                                       kTRUE,
 710                                                       pos,
 711                                                       kFALSE, // Why kFALSE for "isDCA" ? FIXME
 712                                                       covTr,
 713                                                       (Short_t)esdCascadeNeg->GetSign(),
 714                                                       esdCascadeNeg->GetITSClusterMap(),
 715                                                       pid,
 716                                                       vV0FromCascade,
 717                                                       kTRUE,  // usedForVtxFit = kFALSE ? FIXME
 718                                                       vtx->UsesTrack(esdCascadeNeg->GetID()),
 719                                                       AliAODTrack::kSecondary,
 720                                                       selectInfo);
 721         aodTrack->SetTPCFitMap(esdCascadeNeg->GetTPCFitMap());
 722         aodTrack->SetTPCClusterMap(esdCascadeNeg->GetTPCClusterMap());
 723         aodTrack->SetTPCSharedMap (esdCascadeNeg->GetTPCSharedMap());
 724         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeNeg));
 725         aodTrack->SetTPCPointsF(esdCascadeNeg->GetTPCNclsF());
 726         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeNeg->GetTPCCrossedRows()));
 727         fAODTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);
 728
 729         if (esdCascadeNeg->GetSign() > 0) ++fNumberOfPositiveTracks;
 730         aodTrack->ConvertAliPIDtoAODPID();
 731         aodTrack->SetFlags(esdCascadeNeg->GetStatus());
 732         SetAODPID(esdCascadeNeg,aodTrack,detpid);
 733       }
 734       else {
 735         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxNegFromV0Dghter));
 736       }
 737
 738       vV0FromCascade->AddDaughter(aodTrack);
 739
 740
 741       // 4.A.4 - Add the V0 from cascade to the V0 array
 742
 743       Double_t  dcaV0Daughters      = esdV0FromCascade->GetDcaV0Daughters();
 744       Double_t  dcaV0ToPrimVertex   = esdV0FromCascade->GetD( esd.GetPrimaryVertex()->GetX(),
 745                                                              esd.GetPrimaryVertex()->GetY(),
 746                                                              esd.GetPrimaryVertex()->GetZ() );
 747       esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] );
 748       esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] );
 749
 750       Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
 751       dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD(      esd.GetPrimaryVertex()->GetX(),
 752                                                                   esd.GetPrimaryVertex()->GetY(),
 753                                                                   esd.GetMagneticField())        );
 754       dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD(      esd.GetPrimaryVertex()->GetX(),
 755                                                                   esd.GetPrimaryVertex()->GetY(),
 756                                                                   esd.GetMagneticField())        );
 757
 758       aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0( vV0FromCascade,
 759                                                   dcaV0Daughters,
 760                                                   dcaV0ToPrimVertex,
 761                                                   momPosAtV0vtx,
 762                                                   momNegAtV0vtx,
 763                                                   dcaDaughterToPrimVertex);
 764       // set the aod v0 on-the-fly status
 765       aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());
 766
 767       // Add the aodV0 to the aodVORefs
 768       fAODV0Refs->AddAt(aodV0,idxV0FromCascade);
 769
 770       fUsedV0[idxV0FromCascade] = kTRUE;
 771
 772     } else {
 773       // 4.B - if V0 structure already used
 774
 775       // Note :
 776       //    one V0 can be used by several cascades (frequent in PbPb evts) :
 777       // same V0 which used but attached to different bachelor tracks
 778       // -> aodVORefs and fAODV0VtxRefs are needed.
 779       // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.
 780
 781       vV0FromCascade = static_cast<AliAODVertex*>( fAODV0VtxRefs->At(idxV0FromCascade) );
 782       aodV0          = static_cast<AliAODv0*>    ( fAODV0Refs   ->At(idxV0FromCascade) );
 783
 784       // - Treatment of the parent for such a "re-used" V0 :
 785       // Insert the cascade that reuses the V0 vertex in the lineage chain
 786       // Before : vV0 -> vCascade1 -> vPrimary
 787       //  - Hyp : cascade2 uses the same V0 as cascade1
 788       //  After :  vV0 -> vCascade2 -> vCascade1 -> vPrimary
 789
 790       AliAODVertex *vCascadePreviousParent = static_cast<AliAODVertex*> (vV0FromCascade->GetParent());
 791                         vV0FromCascade->SetParent(vCascade);
 792                         vCascade      ->SetParent(vCascadePreviousParent);
 793
 794 //      if(fDebug > 2)
 795 //        printf("---- Cascade / Lineage insertion\n"
 796 //               "Parent of V0 vtx                 = Cascade vtx %p\n"
 797 //               "Parent of the cascade vtx        = Cascade vtx %p\n"
 798 //               "Parent of the parent cascade vtx = Cascade vtx %p\n",
 799 //               static_cast<void*> (vV0FromCascade->GetParent()),
 800 //               static_cast<void*> (vCascade->GetParent()),
 801 //               static_cast<void*> (vCascadePreviousParent->GetParent()) );
 802
 803     }// end if V0 structure already used
 804
 805 //    if (fDebug > 2) {
 806 //      printf("---- Cascade / V0 vertex: \n");
 807 //      vV0FromCascade->Print();
 808 //    }
 809 //
 810 //    if (fDebug > 4) {
 811 //      printf("---- Cascade / pos dghter : \n");
 812 //                      aodTrack->Print();
 813 //      printf("---- Cascade / neg dghter : \n");
 814 //                      aodTrack->Print();
 815 //      printf("---- Cascade / aodV0 : \n");
 816 //                      aodV0->Print();
 817 //    }
 818
 819     // In any case (used V0 or not), add the V0 vertex to the cascade one.
 820     vCascade->AddDaughter(vV0FromCascade);
 821
 822
 823     // 5 - Add the primary track of the cascade (if any)
 824
 825
 826     // 6 - Add the cascade to the AOD array of cascades
 827
 828     Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd.GetPrimaryVertex()->GetX(),
 829                                                                      esd.GetPrimaryVertex()->GetY(),
 830                                                                      esd.GetMagneticField())        );
 831
 832     Double_t momBachAtCascadeVtx[3]={0.};
 833
 834     esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);
 835
 836     aodCascade = new(cascadesArray[fNumberOfCascades++]) AliAODcascade(  vCascade,
 837                                                           esdCascade->Charge(),
 838                                                           esdCascade->GetDcaXiDaughters(),
 839                                                           -999.,
 840                                                           // DCAXiToPrimVtx -> needs to be calculated   ----|
 841                                                           // doesn't exist at ESD level;
 842                                                           // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)
 843                                                           dcaBachToPrimVertexXY,
 844                                                           momBachAtCascadeVtx,
 845                                                           *aodV0);
 846
 847     if (fDebug > 10) {
 848       printf("---- Cascade / AOD cascade : \n\n");
 849       aodCascade->PrintXi(fPrimaryVertex->GetX(), fPrimaryVertex->GetY(), fPrimaryVertex->GetZ());
 850     }
 851
 852   } // end of the loop on cascades
 853
 854   Cascades().Expand(fNumberOfCascades);
 855 }
 856
 857 //______________________________________________________________________________
 858 void AliAnalysisTaskESDfilter::ConvertV0s(const AliESDEvent& esd)
 859 {
 860   // Access to the AOD container of V0s
 861
 862   AliCodeTimerAuto("",0);
 863
 864   //
 865   // V0s
 866   //
 867
 868   Double_t pos[3] = { 0. };
 869   Double_t chi2(0.0);
 870   Double_t covVtx[6] = { 0. };
 871   Double_t momPos[3]={0.};
 872   Double_t covTr[21]={0.};
 873   Double_t pid[10]={0.};
 874   AliAODTrack* aodTrack(0x0);
 875   AliAODPid* detpid(0x0);
 876   Double_t momNeg[3]={0.};
 877   Double_t momPosAtV0vtx[3]={0.};
 878   Double_t momNegAtV0vtx[3]={0.};
 879
 880   for (Int_t nV0 = 0; nV0 < esd.GetNumberOfV0s(); ++nV0)
 881   {
 882     if (fUsedV0[nV0]) continue; // skip if already added to the AOD
 883
 884     AliESDv0 *v0 = esd.GetV0(nV0);
 885     Int_t posFromV0 = v0->GetPindex();
 886     Int_t negFromV0 = v0->GetNindex();
 887
 888     // V0 selection
 889     //
 890     AliESDVertex *esdVtx   = new AliESDVertex(*(esd.GetPrimaryVertex()));
 891     AliESDtrack  *esdV0Pos = esd.GetTrack(posFromV0);
 892     AliESDtrack  *esdV0Neg = esd.GetTrack(negFromV0);
 893     TList v0objects;
 894     v0objects.AddAt(v0,                      0);
 895     v0objects.AddAt(esdV0Pos,                1);
 896     v0objects.AddAt(esdV0Neg,                2);
 897     v0objects.AddAt(esdVtx,                  3);
 898     UInt_t selectV0 = 0;
 899     if (fV0Filter) {
 900       selectV0 = fV0Filter->IsSelected(&v0objects);
 901       // this is a little awkward but otherwise the
 902       // list wants to access the pointer (delete it)
 903       // again when going out of scope
 904       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 905       esdVtx = 0;
 906       if (!selectV0)
 907         continue;
 908     }
 909     else{
 910       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
 911       esdVtx = 0;
 912     }
 913
 914     v0->GetXYZ(pos[0], pos[1], pos[2]);
 915
 916     if (!fOldESDformat) {
 917             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
 918             v0->GetPosCov(covVtx);
 919     } else {
 920             chi2 = -999.;
 921             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 922     }
 923
 924
 925     AliAODVertex * vV0 =
 926           new(Vertices()[fNumberOfVertices++]) AliAODVertex(pos,
 927                                                       covVtx,
 928                                                       chi2,
 929                                                       fPrimaryVertex,
 930                                                       nV0,
 931                                                       AliAODVertex::kV0);
 932     fPrimaryVertex->AddDaughter(vV0);
 933
 934
 935     // Add the positive tracks from the V0
 936
 937
 938     esdV0Pos->GetPxPyPz(momPos);
 939     esdV0Pos->GetXYZ(pos);
 940     esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
 941     esdV0Pos->GetESDpid(pid);
 942
 943     const AliESDVertex *vtx = esd.GetPrimaryVertex();
 944
 945     if (!fUsedTrack[posFromV0]) {
 946             fUsedTrack[posFromV0] = kTRUE;
 947             UInt_t selectInfo = 0;
 948             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);
 949             if(fMChandler)fMChandler->SelectParticle(esdV0Pos->GetLabel());
 950             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Pos->GetID(),
 951                                                     esdV0Pos->GetLabel(),
 952                                                     momPos,
 953                                                     kTRUE,
 954                                                     pos,
 955                                                     kFALSE,
 956                                                     covTr,
 957                                                     (Short_t)esdV0Pos->GetSign(),
 958                                                     esdV0Pos->GetITSClusterMap(),
 959                                                     pid,
 960                                                     vV0,
 961                                                     kTRUE,  // check if this is right
 962                                                     vtx->UsesTrack(esdV0Pos->GetID()),
 963                                                     AliAODTrack::kSecondary,
 964                                                     selectInfo);
 965             aodTrack->SetTPCFitMap(esdV0Pos->GetTPCFitMap());
 966             aodTrack->SetTPCClusterMap(esdV0Pos->GetTPCClusterMap());
 967             aodTrack->SetTPCSharedMap (esdV0Pos->GetTPCSharedMap());
 968             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Pos));
 969             aodTrack->SetTPCPointsF(esdV0Pos->GetTPCNclsF());
 970             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Pos->GetTPCCrossedRows()));
 971             fAODTrackRefs->AddAt(aodTrack,posFromV0);
 972             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
 973             if (esdV0Pos->GetSign() > 0) ++fNumberOfPositiveTracks;
 974             aodTrack->ConvertAliPIDtoAODPID();
 975             aodTrack->SetFlags(esdV0Pos->GetStatus());
 976       SetAODPID(esdV0Pos,aodTrack,detpid);
 977     }
 978     else {
 979             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));
 980             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
 981     }
 982     vV0->AddDaughter(aodTrack);
 983
 984     // Add the negative tracks from the V0
 985
 986     esdV0Neg->GetPxPyPz(momNeg);
 987     esdV0Neg->GetXYZ(pos);
 988     esdV0Neg->GetCovarianceXYZPxPyPz(covTr);
 989     esdV0Neg->GetESDpid(pid);
 990
 991     if (!fUsedTrack[negFromV0]) {
 992             fUsedTrack[negFromV0] = kTRUE;
 993             UInt_t selectInfo = 0;
 994             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);
 995             if(fMChandler)fMChandler->SelectParticle(esdV0Neg->GetLabel());
 996             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Neg->GetID(),
 997                                                     esdV0Neg->GetLabel(),
 998                                                     momNeg,
 999                                                     kTRUE,
1000                                                     pos,
1001                                                     kFALSE,
1002                                                     covTr,
1003                                                     (Short_t)esdV0Neg->GetSign(),
1004                                                     esdV0Neg->GetITSClusterMap(),
1005                                                     pid,
1006                                                     vV0,
1007                                                     kTRUE,  // check if this is right
1008                                                     vtx->UsesTrack(esdV0Neg->GetID()),
1009                                                     AliAODTrack::kSecondary,
1010                                                     selectInfo);
1011             aodTrack->SetTPCFitMap(esdV0Neg->GetTPCFitMap());
1012             aodTrack->SetTPCClusterMap(esdV0Neg->GetTPCClusterMap());
1013             aodTrack->SetTPCSharedMap (esdV0Neg->GetTPCSharedMap());
1014             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Neg));
1015             aodTrack->SetTPCPointsF(esdV0Neg->GetTPCNclsF());
1016             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Neg->GetTPCCrossedRows()));
1017
1018             fAODTrackRefs->AddAt(aodTrack,negFromV0);
1019             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
1020             if (esdV0Neg->GetSign() > 0) ++fNumberOfPositiveTracks;
1021             aodTrack->ConvertAliPIDtoAODPID();
1022             aodTrack->SetFlags(esdV0Neg->GetStatus());
1023       SetAODPID(esdV0Neg,aodTrack,detpid);
1024     }
1025     else {
1026             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(negFromV0));
1027             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
1028     }
1029     vV0->AddDaughter(aodTrack);
1030
1031
1032     // Add the V0 the V0 array as well
1033
1034     Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();
1035     Double_t  dcaV0ToPrimVertex   = v0->GetD(esd.GetPrimaryVertex()->GetX(),
1036                                              esd.GetPrimaryVertex()->GetY(),
1037                                              esd.GetPrimaryVertex()->GetZ());
1038
1039     v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]);
1040     v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]);
1041
1042     Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
1043     dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD(  esd.GetPrimaryVertex()->GetX(),
1044                                                            esd.GetPrimaryVertex()->GetY(),
1045                                                            esd.GetMagneticField()) );
1046     dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD(  esd.GetPrimaryVertex()->GetX(),
1047                                                            esd.GetPrimaryVertex()->GetY(),
1048                                                            esd.GetMagneticField()) );
1049
1050     AliAODv0* aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0(vV0,
1051                                                 dcaV0Daughters,
1052                                                 dcaV0ToPrimVertex,
1053                                                 momPosAtV0vtx,
1054                                                 momNegAtV0vtx,
1055                                                 dcaDaughterToPrimVertex);
1056
1057     // set the aod v0 on-the-fly status
1058     aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
1059   }//End of loop on V0s
1060
1061   V0s().Expand(fNumberOfV0s);
1062 }
1063
1064 //______________________________________________________________________________
1065 void AliAnalysisTaskESDfilter::ConvertTPCOnlyTracks(const AliESDEvent& esd)
1066 {
1067   // Convert TPC only tracks
1068   // Here we have wo hybrid appraoch to remove fakes
1069   // ******* ITSTPC ********
1070   // Uses a cut on the ITS properties to select global tracks
1071   // which are than marked as HybdridITSTPC for the remainder
1072   // the TPC only tracks are flagged as HybridITSTPConly.
1073   // Note, in order not to get fakes back in the TPC cuts, one needs
1074   // two "ITS" cuts one tight (1) (to throw out fakes) and one lose (2) (to NOT flag the trakcs in the TPC only)
1075   // using cut number (3)
1076   // so fHybridFilterMask == (1)|(2) fTPCFilterMask = (3), Usercode needs to slect with mask = (1)|(3) and track->IsHybridITSTPC()
1077   // ******* TPC ********
1078   // Here, only TPC tracks are flagged that pass the tight ITS cuts and tracks that pass the TPC cuts and NOT the loose ITS cuts
1079   // the ITS cuts neeed to be added to the filter as extra cuts, since here the selections info is reset in the global and put to the TPC only track
1080
1081   AliCodeTimerAuto("",0);
1082
1083   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1084   for(int it = 0;it < fNumberOfTracks;++it)
1085   {
1086     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1087     if(!tr)continue;
1088     UInt_t map = tr->GetFilterMap();
1089     if(map&fTPCConstrainedFilterMask){
1090       // we only reset the track select ionfo, no deletion...
1091       tr->SetFilterMap(map&~fTPCConstrainedFilterMask);
1092     }
1093     if(map&fHybridFilterMaskTPCCG){
1094       // this is one part of the hybrid tracks
1095       // the others not passing the selection will be TPC only selected below
1096       tr->SetIsHybridTPCConstrainedGlobal(kTRUE);
1097     }
1098   }
1099   // Loop over the ESD trcks and pick out the tracks passing TPC only cuts
1100
1101
1102   const AliESDVertex *vtxSPD = esd.GetPrimaryVertexSPD();
1103   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1104
1105   Double_t pos[3] = { 0. };
1106   Double_t covTr[21]={0.};
1107   Double_t pid[10]={0.};
1108
1109
1110   Double_t p[3] = { 0. };
1111
1112   Double_t pDCA[3] = { 0. }; // momentum at DCA
1113   Double_t rDCA[3] = { 0. }; // position at DCA
1114   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1115   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1116
1117
1118   AliAODTrack* aodTrack(0x0);
1119   //  AliAODPid* detpid(0x0);
1120
1121   // account for change in pT after the constraint
1122   Float_t ptMax = 1E10;
1123   Float_t ptMin = 0;
1124   for(int i = 0;i<32;i++){
1125     if(fTPCConstrainedFilterMask&(1<<i)){
1126       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1127       Float_t tmp1= 0,tmp2 = 0;
1128       cuts->GetPtRange(tmp1,tmp2);
1129       if(tmp1>ptMin)ptMin=tmp1;
1130       if(tmp2<ptMax)ptMax=tmp2;
1131     }
1132   }
1133
1134   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1135   {
1136     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1137
1138     UInt_t selectInfo = 0;
1139     Bool_t isHybridITSTPC = false;
1140     //
1141     // Track selection
1142     if (fTrackFilter) {
1143       selectInfo = fTrackFilter->IsSelected(esdTrack);
1144     }
1145
1146     if(!(selectInfo&fHybridFilterMaskTPCCG)){
1147       // not already selected tracks, use second part of hybrid tracks
1148       isHybridITSTPC = true;
1149       // too save space one could only store these...
1150     }
1151
1152     selectInfo &= fTPCConstrainedFilterMask;
1153     if (!selectInfo)continue;
1154     if (fWriteHybridTPCCOnly&&!isHybridITSTPC)continue; // write only complementary tracks
1155     // create a tpc only tracl
1156     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(&esd),esdTrack->GetID());
1157     if(!track) continue;
1158
1159     if(track->Pt()>0.)
1160     {
1161       // only constrain tracks above threshold
1162       AliExternalTrackParam exParam;
1163       // take the B-field from the ESD, no 3D fieldMap available at this point
1164       Bool_t relate = false;
1165       relate = track->RelateToVertexTPC(vtxSPD,esd.GetMagneticField(),kVeryBig,&exParam);
1166       if(!relate){
1167         delete track;
1168         continue;
1169       }
1170       // fetch the track parameters at the DCA (unconstraint)
1171       if(track->GetTPCInnerParam()){
1172         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1173         track->GetTPCInnerParam()->GetXYZ(rDCA);
1174       }
1175       // get the DCA to the vertex:
1176       track->GetImpactParametersTPC(dDCA,cDCA);
1177       // set the constrained parameters to the track
1178       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1179     }
1180
1181     track->GetPxPyPz(p);
1182
1183     Float_t pT = track->Pt();
1184     if(pT<ptMin||pT>ptMax){
1185       delete track;
1186       continue;
1187     }
1188
1189     //
1190
1191
1192     track->GetXYZ(pos);
1193     track->GetCovarianceXYZPxPyPz(covTr);
1194     esdTrack->GetESDpid(pid);// original PID
1195
1196     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1197     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((track->GetID()+1)*-1,
1198                                                             track->GetLabel(),
1199                                                             p,
1200                                                             kTRUE,
1201                                                             pos,
1202                                                             kFALSE,
1203                                                             covTr,
1204                                                             (Short_t)track->GetSign(),
1205                                                             track->GetITSClusterMap(),
1206                                                             pid,
1207                                                             fPrimaryVertex,
1208                                                             kTRUE, // check if this is right
1209                                                             vtx->UsesTrack(track->GetID()),
1210                                                             AliAODTrack::kPrimary,
1211                                                             selectInfo);
1212     aodTrack->SetIsHybridTPCConstrainedGlobal(isHybridITSTPC);
1213     aodTrack->SetTPCFitMap(track->GetTPCFitMap());
1214     aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
1215     aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
1216     aodTrack->SetIsTPCConstrained(kTRUE);
1217     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack)); // original track
1218     // set the DCA values to the AOD track
1219     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1220     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1221     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1222
1223     aodTrack->SetFlags(track->GetStatus());
1224     aodTrack->SetTPCPointsF(track->GetTPCNclsF());
1225     aodTrack->SetTPCNCrossedRows(UShort_t(track->GetTPCCrossedRows()));
1226
1227     //Perform progagation of tracks if needed
1228     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1229     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1230
1231     // do not duplicate PID information
1232     //    aodTrack->ConvertAliPIDtoAODPID();
1233     //    SetAODPID(esdTrack,aodTrack,detpid);
1234
1235     delete track;
1236   } // end of loop on tracks
1237
1238 }
1239
1240 //______________________________________________________________________________
1241 void AliAnalysisTaskESDfilter::ConvertGlobalConstrainedTracks(const AliESDEvent& esd)
1242 {
1243
1244   // Here we have the option to store the complement from global constraint information
1245   // to tracks passing tight cuts (1) in order not to get fakes back in, one needs
1246   // two sets of cuts one tight (1) (to throw out fakes) and one lose (2) (fakes/bad tracks would pass (2) but not (1))
1247   // using cut number (3) selects the tracks that complement (1) e.g. tracks witout ITS refit or cluster requirement
1248
1249
1250   AliCodeTimerAuto("",0);
1251
1252   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1253   for(int it = 0;it < fNumberOfTracks;++it)
1254   {
1255     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1256     if(!tr)continue;
1257     UInt_t map = tr->GetFilterMap();
1258     if(map&fGlobalConstrainedFilterMask){
1259       // we only reset the track select info, no deletion...
1260       // mask reset mask in case track is already taken
1261       tr->SetFilterMap(map&~fGlobalConstrainedFilterMask);
1262     }
1263     if(map&fHybridFilterMaskGCG){
1264       // this is one part of the hybrid tracks
1265       // the others not passing the selection will be the ones selected below
1266       tr->SetIsHybridGlobalConstrainedGlobal(kTRUE);
1267     }
1268   }
1269   // Loop over the ESD trcks and pick out the tracks passing the GlobalConstraint cuts
1270
1271
1272   Double_t pos[3] = { 0. };
1273   Double_t covTr[21]={0.};
1274   Double_t pid[10]={0.};
1275   Double_t p[3] = { 0. };
1276
1277   Double_t pDCA[3] = { 0. }; // momentum at DCA
1278   Double_t rDCA[3] = { 0. }; // position at DCA
1279   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1280   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1281
1282
1283   AliAODTrack* aodTrack(0x0);
1284   AliAODPid* detpid(0x0);
1285   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1286
1287   // account for change in pT after the constraint
1288   Float_t ptMax = 1E10;
1289   Float_t ptMin = 0;
1290   for(int i = 0;i<32;i++){
1291     if(fGlobalConstrainedFilterMask&(1<<i)){
1292       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1293       Float_t tmp1= 0,tmp2 = 0;
1294       cuts->GetPtRange(tmp1,tmp2);
1295       if(tmp1>ptMin)ptMin=tmp1;
1296       if(tmp2<ptMax)ptMax=tmp2;
1297     }
1298   }
1299
1300
1301
1302  for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1303   {
1304     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy
1305     const AliExternalTrackParam * exParamGC = esdTrack->GetConstrainedParam();
1306     if(!exParamGC)continue;
1307
1308     UInt_t selectInfo = 0;
1309     Bool_t isHybridGC = false;
1310
1311     //
1312     // Track selection
1313     if (fTrackFilter) {
1314       selectInfo = fTrackFilter->IsSelected(esdTrack);
1315     }
1316
1317
1318     if(!(selectInfo&fHybridFilterMaskGCG))isHybridGC = true;
1319     if (fWriteHybridGCOnly&&!isHybridGC)continue; // write only complementary tracks
1320
1321     selectInfo &= fGlobalConstrainedFilterMask;
1322     if (!selectInfo)continue;
1323     // fetch the track parameters at the DCA (unconstrained)
1324     esdTrack->GetPxPyPz(pDCA);
1325     esdTrack->GetXYZ(rDCA);
1326     // get the DCA to the vertex:
1327     esdTrack->GetImpactParameters(dDCA,cDCA);
1328
1329     if (!esdTrack->GetConstrainedPxPyPz(p)) continue;
1330
1331
1332     Float_t pT = exParamGC->Pt();
1333     if(pT<ptMin||pT>ptMax){
1334       continue;
1335     }
1336
1337
1338     esdTrack->GetConstrainedXYZ(pos);
1339     exParamGC->GetCovarianceXYZPxPyPz(covTr);
1340     esdTrack->GetESDpid(pid);
1341     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1342     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((esdTrack->GetID()+1)*-1,
1343                                                             esdTrack->GetLabel(),
1344                                                             p,
1345                                                             kTRUE,
1346                                                             pos,
1347                                                             kFALSE,
1348                                                             covTr,
1349                                                             (Short_t)esdTrack->GetSign(),
1350                                                             esdTrack->GetITSClusterMap(),
1351                                                             pid,
1352                                                             fPrimaryVertex,
1353                                                             kTRUE, // check if this is right
1354                                                             vtx->UsesTrack(esdTrack->GetID()),
1355                                                             AliAODTrack::kPrimary,
1356                                                             selectInfo);
1357     aodTrack->SetIsHybridGlobalConstrainedGlobal(isHybridGC);
1358     aodTrack->SetIsGlobalConstrained(kTRUE);
1359     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1360     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1361     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1362     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1363
1364
1365     // set the DCA values to the AOD track
1366     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1367     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1368     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1369
1370     aodTrack->SetFlags(esdTrack->GetStatus());
1371     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1372     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1373
1374     if(isHybridGC){
1375       // only copy AOD information for hybrid, no duplicate information
1376       aodTrack->ConvertAliPIDtoAODPID();
1377       SetAODPID(esdTrack,aodTrack,detpid);
1378     }
1379
1380     //Perform progagation of tracks if needed
1381     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1382     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1383   } // end of loop on tracks
1384
1385 }
1386
1387
1388 //______________________________________________________________________________
1389 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)
1390 {
1391   // Tracks (primary and orphan)
1392
1393   AliCodeTimerAuto("",0);
1394
1395   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));
1396
1397   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1398   Double_t p[3] = { 0. };
1399   Double_t pos[3] = { 0. };
1400   Double_t covTr[21] = { 0. };
1401   Double_t pid[10] = { 0. };
1402   AliAODTrack* aodTrack(0x0);
1403   AliAODPid* detpid(0x0);
1404
1405   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack)
1406   {
1407     if (fUsedTrack[nTrack]) continue;
1408
1409     AliESDtrack *esdTrack = esd.GetTrack(nTrack);
1410     UInt_t selectInfo = 0;
1411     //
1412     // Track selection
1413     if (fTrackFilter) {
1414             selectInfo = fTrackFilter->IsSelected(esdTrack);
1415             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;
1416     }
1417
1418
1419     esdTrack->GetPxPyPz(p);
1420     esdTrack->GetXYZ(pos);
1421     esdTrack->GetCovarianceXYZPxPyPz(covTr);
1422     esdTrack->GetESDpid(pid);
1423     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1424     fPrimaryVertex->AddDaughter(aodTrack =
1425                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),
1426                                                             esdTrack->GetLabel(),
1427                                                             p,
1428                                                             kTRUE,
1429                                                             pos,
1430                                                             kFALSE,
1431                                                             covTr,
1432                                                             (Short_t)esdTrack->GetSign(),
1433                                                             esdTrack->GetITSClusterMap(),
1434                                                             pid,
1435                                                             fPrimaryVertex,
1436                                                             kTRUE, // check if this is right
1437                                                             vtx->UsesTrack(esdTrack->GetID()),
1438                                                             AliAODTrack::kPrimary,
1439                                                             selectInfo)
1440                          );
1441     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1442     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1443     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1444     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1445     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1446     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1447     if(esdTrack->IsEMCAL()) aodTrack->SetEMCALcluster(esdTrack->GetEMCALcluster());
1448     if(esdTrack->IsPHOS())  aodTrack->SetPHOScluster(esdTrack->GetPHOScluster());
1449
1450     //Perform progagation of tracks if needed
1451     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1452     aodTrack->SetTrackPhiEtaOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal());
1453
1454     fAODTrackRefs->AddAt(aodTrack, nTrack);
1455
1456
1457     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;
1458     aodTrack->SetFlags(esdTrack->GetStatus());
1459     aodTrack->ConvertAliPIDtoAODPID();
1460     SetAODPID(esdTrack,aodTrack,detpid);
1461   } // end of loop on tracks
1462 }
1463
1464 //______________________________________________________________________________
1465 void AliAnalysisTaskESDfilter::PropagateTrackToEMCal(AliESDtrack *esdTrack)
1466 {
1467   Double_t trkPos[3] = {0.,0.,0.};
1468   Double_t EMCalEta=-999, EMCalPhi=-999;
1469   Double_t trkphi = esdTrack->Phi()*TMath::RadToDeg();
1470   if(TMath::Abs(esdTrack->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
1471     {
1472       AliExternalTrackParam *trkParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1473       if(trkParam)
1474         {
1475           AliExternalTrackParam trkParamTmp(*trkParam);
1476           if(AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, 430, esdTrack->GetMass(), 20, kTRUE, 0.8, -1))
1477             {
1478               trkParamTmp.GetXYZ(trkPos);
1479               TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1480               EMCalEta = trkPosVec.Eta();
1481               EMCalPhi = trkPosVec.Phi();
1482               if(EMCalPhi<0)  EMCalPhi += 2*TMath::Pi();
1483               esdTrack->SetTrackPhiEtaOnEMCal(EMCalPhi,EMCalEta);
1484             }
1485         }
1486     }
1487 }
1488
1489 //______________________________________________________________________________
1490 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)
1491 {
1492 // Convert PMD Clusters
1493   AliCodeTimerAuto("",0);
1494   Int_t jPmdClusters=0;
1495   // Access to the AOD container of PMD clusters
1496   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
1497   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {
1498     // file pmd clusters, to be revised!
1499     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);
1500     Int_t nLabel = 0;
1501     Int_t *label = 0x0;
1502     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
1503     Double_t pidPmd[13] = { 0.}; // to be revised!
1504     // type not set!
1505     // assoc cluster not set
1506     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
1507   }
1508 }
1509
1510
1511 //______________________________________________________________________________
1512 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)
1513 {
1514 // Convert Calorimeter Clusters
1515   AliCodeTimerAuto("",0);
1516
1517   // Access to the AOD container of clusters
1518   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
1519   Int_t jClusters(0);
1520
1521   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
1522
1523     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
1524
1525     Int_t  id        = cluster->GetID();
1526     Int_t  nLabel    = cluster->GetNLabels();
1527     Int_t *labels    = cluster->GetLabels();
1528     if(labels){
1529                   for(int i = 0;i < nLabel;++i){
1530                           if(fMChandler)fMChandler->SelectParticle(labels[i]);
1531                   }
1532           }
1533
1534     Float_t energy = cluster->E();
1535     Float_t posF[3] = { 0.};
1536     cluster->GetPosition(posF);
1537
1538     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
1539                                                                                       nLabel,
1540                                                                                       labels,
1541                                                                                       energy,
1542                                                                                       posF,
1543                                                                                       NULL,
1544                                                                                       cluster->GetType(),0);
1545
1546     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
1547                                 cluster->GetDispersion(),
1548                                 cluster->GetM20(), cluster->GetM02(),
1549                                 cluster->GetEmcCpvDistance(),
1550                                 cluster->GetNExMax(),cluster->GetTOF()) ;
1551
1552     caloCluster->SetPIDFromESD(cluster->GetPID());
1553     caloCluster->SetNCells(cluster->GetNCells());
1554     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
1555     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
1556
1557     caloCluster->SetTrackDistance(cluster->GetTrackDx(), cluster->GetTrackDz());
1558
1559     Int_t nMatchCount = 0;
1560     TArrayI* matchedT =         cluster->GetTracksMatched();
1561     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
1562       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
1563         Int_t iESDtrack = matchedT->At(im);;
1564         if (fAODTrackRefs->At(iESDtrack) != 0) {
1565           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
1566           nMatchCount++;
1567         }
1568       }
1569     }
1570     if(nMatchCount==0)
1571       caloCluster->SetTrackDistance(-999,-999);
1572
1573   }
1574   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
1575 }
1576
1577 //______________________________________________________________________________
1578 void AliAnalysisTaskESDfilter::ConvertCaloTrigger(TString calo, const AliESDEvent& esd)
1579 {
1580         AliCodeTimerAuto("",0);
1581
1582         if (calo == "PHOS")
1583         {
1584           AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1585           AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1586
1587           aodTrigger.Allocate(esdTrigger.GetEntries());
1588           esdTrigger.Reset();
1589
1590           Float_t a;
1591           Int_t tmod,tabsId;
1592
1593           while (esdTrigger.Next()) {
1594             esdTrigger.GetPosition(tmod,tabsId);
1595             esdTrigger.GetAmplitude(a);
1596             aodTrigger.Add(tmod,tabsId,a,0.,(Int_t*)NULL,0,0,0);
1597           }
1598
1599           return;
1600         }
1601
1602         AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1603
1604         if (aodHandler)
1605         {
1606                 TTree *aodTree = aodHandler->GetTree();
1607
1608                 if (aodTree)
1609                 {
1610                         Int_t *type = esd.GetCaloTriggerType();
1611
1612                         for (Int_t i = 0; i < 15; i++)
1613                         {
1614                                 aodTree->GetUserInfo()->Add(new TParameter<int>(Form("EMCALCaloTrigger%d",i), type[i]));
1615                         }
1616                 }
1617         }
1618
1619         AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1620
1621         AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1622
1623         aodTrigger.Allocate(esdTrigger.GetEntries());
1624
1625         esdTrigger.Reset();
1626         while (esdTrigger.Next())
1627         {
1628                 Int_t px, py, ts, nTimes, times[10], b;
1629                 Float_t a, t;
1630
1631                 esdTrigger.GetPosition(px, py);
1632
1633                 esdTrigger.GetAmplitude(a);
1634                 esdTrigger.GetTime(t);
1635
1636                 esdTrigger.GetL0Times(times);
1637                 esdTrigger.GetNL0Times(nTimes);
1638
1639                 esdTrigger.GetL1TimeSum(ts);
1640
1641                 esdTrigger.GetTriggerBits(b);
1642
1643                 aodTrigger.Add(px, py, a, t, times, nTimes, ts, b);
1644         }
1645
1646         for (int i = 0; i < 4; i++) aodTrigger.SetL1Threshold(i, esdTrigger.GetL1Threshold(i));
1647
1648         Int_t v0[2] =
1649         {
1650                 esdTrigger.GetL1V0(0),
1651                 esdTrigger.GetL1V0(1)
1652         };
1653
1654         aodTrigger.SetL1V0(v0);
1655         aodTrigger.SetL1FrameMask(esdTrigger.GetL1FrameMask());
1656 }
1657
1658 //______________________________________________________________________________
1659 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)
1660 {
1661 // Convert EMCAL Cells
1662   AliCodeTimerAuto("",0);
1663   // fill EMCAL cell info
1664   if (esd.GetEMCALCells()) { // protection against missing ESD information
1665     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
1666     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
1667
1668     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1669     aodEMcells.CreateContainer(nEMcell);
1670     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
1671     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
1672       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),
1673                          esdEMcells.GetTime(iCell), esdEMcells.GetMCLabel(iCell), esdEMcells.GetEFraction(iCell));
1674     }
1675     aodEMcells.Sort();
1676   }
1677 }
1678
1679 //______________________________________________________________________________
1680 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)
1681 {
1682 // Convert PHOS Cells
1683   AliCodeTimerAuto("",0);
1684   // fill PHOS cell info
1685   if (esd.GetPHOSCells()) { // protection against missing ESD information
1686     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
1687     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
1688
1689     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
1690     aodPHcells.CreateContainer(nPHcell);
1691     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
1692     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
1693       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),
1694                          esdPHcells.GetTime(iCell), esdPHcells.GetMCLabel(iCell), esdPHcells.GetEFraction(iCell));
1695     }
1696     aodPHcells.Sort();
1697   }
1698 }
1699
1700 //______________________________________________________________________________
1701 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)
1702 {
1703   // tracklets
1704   AliCodeTimerAuto("",0);
1705
1706   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());
1707   const AliMultiplicity *mult = esd.GetMultiplicity();
1708   if (mult) {
1709     if (mult->GetNumberOfTracklets()>0) {
1710       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
1711
1712       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
1713         if(fMChandler){
1714           fMChandler->SelectParticle(mult->GetLabel(n, 0));
1715           fMChandler->SelectParticle(mult->GetLabel(n, 1));
1716         }
1717         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1718       }
1719     }
1720   } else {
1721     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
1722   }
1723 }
1724
1725 //______________________________________________________________________________
1726 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)
1727 {
1728   AliCodeTimerAuto("",0);
1729
1730   // Kinks: it is a big mess the access to the information in the kinks
1731   // The loop is on the tracks in order to find the mother and daugther of each kink
1732
1733   Double_t covTr[21]={0.};
1734   Double_t pid[10]={0.};
1735   AliAODPid* detpid(0x0);
1736
1737   fNumberOfKinks = esd.GetNumberOfKinks();
1738
1739   const AliESDVertex* vtx = esd.GetPrimaryVertex();
1740
1741   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
1742   {
1743     AliESDtrack * esdTrack = esd.GetTrack(iTrack);
1744
1745     Int_t ikink = esdTrack->GetKinkIndex(0);
1746
1747     if (ikink && fNumberOfKinks) {
1748             // Negative kink index: mother, positive: daughter
1749
1750             // Search for the second track of the kink
1751
1752             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {
1753
1754         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);
1755
1756         Int_t jkink = esdTrack1->GetKinkIndex(0);
1757
1758         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
1759
1760           // The two tracks are from the same kink
1761
1762           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
1763
1764           Int_t imother = -1;
1765           Int_t idaughter = -1;
1766
1767           if (ikink<0 && jkink>0) {
1768
1769             imother = iTrack;
1770             idaughter = jTrack;
1771           }
1772           else if (ikink>0 && jkink<0) {
1773
1774             imother = jTrack;
1775             idaughter = iTrack;
1776           }
1777           else {
1778             //                  cerr << "Error: Wrong combination of kink indexes: "
1779             //                       << ikink << " " << jkink << endl;
1780             continue;
1781           }
1782
1783           // Add the mother track if it passed primary track selection cuts
1784
1785           AliAODTrack * mother = NULL;
1786
1787           UInt_t selectInfo = 0;
1788           if (fTrackFilter) {
1789             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));
1790             if (!selectInfo) continue;
1791           }
1792
1793           if (!fUsedTrack[imother]) {
1794
1795             fUsedTrack[imother] = kTRUE;
1796
1797             AliESDtrack *esdTrackM = esd.GetTrack(imother);
1798             Double_t p[3] = { 0. };
1799             Double_t pos[3] = { 0. };
1800             esdTrackM->GetPxPyPz(p);
1801             esdTrackM->GetXYZ(pos);
1802             esdTrackM->GetCovarianceXYZPxPyPz(covTr);
1803             esdTrackM->GetESDpid(pid);
1804             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());
1805             mother =
1806             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),
1807                                                esdTrackM->GetLabel(),
1808                                                p,
1809                                                kTRUE,
1810                                                pos,
1811                                                kFALSE,
1812                                                covTr,
1813                                                (Short_t)esdTrackM->GetSign(),
1814                                                esdTrackM->GetITSClusterMap(),
1815                                                pid,
1816                                                fPrimaryVertex,
1817                                                kTRUE, // check if this is right
1818                                                vtx->UsesTrack(esdTrack->GetID()),
1819                                                AliAODTrack::kPrimary,
1820                                                selectInfo);
1821             mother->SetTPCFitMap(esdTrackM->GetTPCFitMap());
1822             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());
1823             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());
1824             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));
1825             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());
1826             mother->SetTPCNCrossedRows(UShort_t(esdTrackM->GetTPCCrossedRows()));
1827
1828             fAODTrackRefs->AddAt(mother, imother);
1829
1830             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;
1831             mother->SetFlags(esdTrackM->GetStatus());
1832             mother->ConvertAliPIDtoAODPID();
1833             fPrimaryVertex->AddDaughter(mother);
1834             mother->ConvertAliPIDtoAODPID();
1835             SetAODPID(esdTrackM,mother,detpid);
1836           }
1837           else {
1838             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1839             //                       << " track " << imother << " has already been used!" << endl;
1840           }
1841
1842           // Add the kink vertex
1843           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);
1844
1845           AliAODVertex * vkink =
1846           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),
1847                                                   NULL,
1848                                                   0.,
1849                                                   mother,
1850                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!
1851                                                   AliAODVertex::kKink);
1852           // Add the daughter track
1853
1854           AliAODTrack * daughter = NULL;
1855
1856           if (!fUsedTrack[idaughter]) {
1857
1858             fUsedTrack[idaughter] = kTRUE;
1859
1860             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);
1861             Double_t p[3] = { 0. };
1862             Double_t pos[3] = { 0. };
1863
1864             esdTrackD->GetPxPyPz(p);
1865             esdTrackD->GetXYZ(pos);
1866             esdTrackD->GetCovarianceXYZPxPyPz(covTr);
1867             esdTrackD->GetESDpid(pid);
1868             selectInfo = 0;
1869             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);
1870             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());
1871             daughter =
1872             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),
1873                                                esdTrackD->GetLabel(),
1874                                                p,
1875                                                kTRUE,
1876                                                pos,
1877                                                kFALSE,
1878                                                covTr,
1879                                                (Short_t)esdTrackD->GetSign(),
1880                                                esdTrackD->GetITSClusterMap(),
1881                                                pid,
1882                                                vkink,
1883                                                kTRUE, // check if this is right
1884                                                vtx->UsesTrack(esdTrack->GetID()),
1885                                                AliAODTrack::kSecondary,
1886                                                selectInfo);
1887             daughter->SetTPCFitMap(esdTrackD->GetTPCFitMap());
1888             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());
1889             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());
1890             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());
1891             daughter->SetTPCNCrossedRows(UShort_t(esdTrackD->GetTPCCrossedRows()));
1892             fAODTrackRefs->AddAt(daughter, idaughter);
1893
1894             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;
1895             daughter->SetFlags(esdTrackD->GetStatus());
1896             daughter->ConvertAliPIDtoAODPID();
1897             vkink->AddDaughter(daughter);
1898             daughter->ConvertAliPIDtoAODPID();
1899             SetAODPID(esdTrackD,daughter,detpid);
1900           }
1901           else {
1902             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1903             //                       << " track " << idaughter << " has already been used!" << endl;
1904           }
1905         }
1906             }
1907     }
1908   }
1909 }
1910
1911 //______________________________________________________________________________
1912 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)
1913 {
1914   AliCodeTimerAuto("",0);
1915
1916   // Access to the AOD container of vertices
1917   fNumberOfVertices = 0;
1918
1919   Double_t pos[3] = { 0. };
1920   Double_t covVtx[6] = { 0. };
1921
1922   // Add primary vertex. The primary tracks will be defined
1923   // after the loops on the composite objects (V0, cascades, kinks)
1924   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1925
1926   vtx->GetXYZ(pos); // position
1927   vtx->GetCovMatrix(covVtx); //covariance matrix
1928
1929   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])
1930   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1931   fPrimaryVertex->SetName(vtx->GetName());
1932   fPrimaryVertex->SetTitle(vtx->GetTitle());
1933   fPrimaryVertex->SetBC(vtx->GetBC());
1934
1935   TString vtitle = vtx->GetTitle();
1936   if (!vtitle.Contains("VertexerTracks"))
1937     fPrimaryVertex->SetNContributors(vtx->GetNContributors());
1938
1939   if (fDebug > 0) fPrimaryVertex->Print();
1940
1941   // Add SPD "main" vertex
1942   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
1943   vtxS->GetXYZ(pos); // position
1944   vtxS->GetCovMatrix(covVtx); //covariance matrix
1945   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])
1946   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
1947   mVSPD->SetName(vtxS->GetName());
1948   mVSPD->SetTitle(vtxS->GetTitle());
1949   mVSPD->SetNContributors(vtxS->GetNContributors());
1950
1951   // Add SPD pileup vertices
1952   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
1953   {
1954     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
1955     vtxP->GetXYZ(pos); // position
1956     vtxP->GetCovMatrix(covVtx); //covariance matrix
1957     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])
1958     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
1959     pVSPD->SetName(vtxP->GetName());
1960     pVSPD->SetTitle(vtxP->GetTitle());
1961     pVSPD->SetNContributors(vtxP->GetNContributors());
1962     pVSPD->SetBC(vtxP->GetBC());
1963   }
1964
1965   // Add TRK pileup vertices
1966   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
1967   {
1968     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
1969     vtxP->GetXYZ(pos); // position
1970     vtxP->GetCovMatrix(covVtx); //covariance matrix
1971     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])
1972     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
1973     pVTRK->SetName(vtxP->GetName());
1974     pVTRK->SetTitle(vtxP->GetTitle());
1975     pVTRK->SetNContributors(vtxP->GetNContributors());
1976     pVTRK->SetBC(vtxP->GetBC());
1977   }
1978 }
1979
1980 //______________________________________________________________________________
1981 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)
1982 {
1983   // Convert VZERO data
1984   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
1985   *vzeroData = *(esd.GetVZEROData());
1986 }
1987
1988 //______________________________________________________________________________
1989 void AliAnalysisTaskESDfilter::ConvertTZERO(const AliESDEvent& esd)
1990 {
1991   // Convert TZERO data
1992   const AliESDTZERO* esdTzero = esd.GetESDTZERO();
1993   AliAODTZERO* aodTzero = AODEvent()->GetTZEROData();
1994
1995   for (Int_t icase=0; icase<3; icase++){
1996     aodTzero->SetT0TOF(    icase, esdTzero->GetT0TOF(icase));
1997     aodTzero->SetT0TOFbest(icase, esdTzero->GetT0TOFbest(icase));
1998   }
1999   aodTzero->SetBackgroundFlag(esdTzero->GetBackgroundFlag());
2000   aodTzero->SetPileupFlag(esdTzero->GetPileupFlag());
2001   aodTzero->SetSatelliteFlag(esdTzero->GetSatellite());
2002
2003   Float_t rawTime[24];
2004   for(Int_t ipmt=0; ipmt<24; ipmt++)
2005     rawTime[ipmt] = esdTzero->GetTimeFull(ipmt,0);
2006
2007   Int_t idxOfFirstPmtA = -1,       idxOfFirstPmtC = -1;
2008   Float_t timeOfFirstPmtA = 9999, timeOfFirstPmtC = 9999;
2009   for(int ipmt=0;  ipmt<12; ipmt++){
2010     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtC && rawTime[ipmt]!=0){
2011       timeOfFirstPmtC = rawTime[ipmt];
2012       idxOfFirstPmtC  = ipmt;
2013     }
2014   }
2015   for(int ipmt=12; ipmt<24; ipmt++){
2016     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtA && rawTime[ipmt]!=0 ){
2017       timeOfFirstPmtA = rawTime[ipmt];
2018       idxOfFirstPmtA  = ipmt;
2019     }
2020   }
2021
2022   if(idxOfFirstPmtA != -1 && idxOfFirstPmtC != -1){
2023     //speed of light in cm/ns   TMath::C()*1e-7
2024     Float_t vertexraw = TMath::C()*1e-7 * (rawTime[idxOfFirstPmtA] - rawTime[idxOfFirstPmtC])/2;
2025     aodTzero->SetT0VertexRaw( vertexraw );
2026   }else{
2027     aodTzero->SetT0VertexRaw(99999);
2028   }
2029
2030   aodTzero->SetT0zVertex(esdTzero->GetT0zVertex());
2031 }
2032
2033
2034 //______________________________________________________________________________
2035 void AliAnalysisTaskESDfilter::ConvertZDC(const AliESDEvent& esd)
2036 {
2037   // Convert ZDC data
2038   AliESDZDC* esdZDC = esd.GetZDCData();
2039
2040   const Double_t zem1Energy = esdZDC->GetZEM1Energy();
2041   const Double_t zem2Energy = esdZDC->GetZEM2Energy();
2042
2043   const Double_t *towZNC = esdZDC->GetZNCTowerEnergy();
2044   const Double_t *towZPC = esdZDC->GetZPCTowerEnergy();
2045   const Double_t *towZNA = esdZDC->GetZNATowerEnergy();
2046   const Double_t *towZPA = esdZDC->GetZPATowerEnergy();
2047   const Double_t *towZNCLG = esdZDC->GetZNCTowerEnergyLR();
2048   const Double_t *towZNALG = esdZDC->GetZNATowerEnergyLR();
2049
2050   AliAODZDC* zdcAOD = AODEvent()->GetZDCData();
2051
2052   zdcAOD->SetZEM1Energy(zem1Energy);
2053   zdcAOD->SetZEM2Energy(zem2Energy);
2054   zdcAOD->SetZNCTowers(towZNC, towZNCLG);
2055   zdcAOD->SetZNATowers(towZNA, towZNALG);
2056   zdcAOD->SetZPCTowers(towZPC);
2057   zdcAOD->SetZPATowers(towZPA);
2058
2059   zdcAOD->SetZDCParticipants(esdZDC->GetZDCParticipants(), esdZDC->GetZDCPartSideA(), esdZDC->GetZDCPartSideC());
2060   zdcAOD->SetZDCImpactParameter(esdZDC->GetImpactParameter(), esdZDC->GetImpactParamSideA(),
2061         esdZDC->GetImpactParamSideC());
2062   zdcAOD->SetZDCTDCSum(esdZDC->GetZNTDCSum(0));
2063   zdcAOD->SetZDCTDCDiff(esdZDC->GetZNTDCDiff(0));
2064
2065 }
2066
2067 //_______________________________________________________________________________________________________________________________________
2068 Int_t AliAnalysisTaskESDfilter::ConvertHMPID(const AliESDEvent& esd) // clm
2069 {
2070   //
2071   // Convtert ESD HMPID info to AOD and return the number of good tracks with HMPID signal.
2072   // We need to return an int since there is no signal counter in the ESD.
2073   //
2074
2075   AliCodeTimerAuto("",0);
2076
2077   Int_t cntHmpidGoodTracks = 0;
2078
2079   Float_t  xMip = 0;
2080   Float_t  yMip = 0;
2081   Int_t    qMip = 0;
2082   Int_t    nphMip = 0;
2083
2084   Float_t xTrk = 0;
2085   Float_t yTrk = 0;
2086   Float_t thetaTrk = 0;
2087   Float_t phiTrk = 0;
2088
2089   Double_t hmpPid[5]={0};
2090   Double_t hmpMom[3]={0};
2091
2092   TClonesArray &hmpidRings = *(AODEvent()->GetHMPIDrings());
2093
2094  for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack)
2095   {
2096     if(! esd.GetTrack(iTrack) ) continue;
2097
2098     if(esd.GetTrack(iTrack)->GetHMPIDsignal() > -20 ) {                  //
2099
2100       (esd.GetTrack(iTrack))->GetHMPIDmip(xMip, yMip, qMip, nphMip);     // Get MIP properties
2101       (esd.GetTrack(iTrack))->GetHMPIDtrk(xTrk,yTrk,thetaTrk,phiTrk);
2102       (esd.GetTrack(iTrack))->GetHMPIDpid(hmpPid);
2103       if((esd.GetTrack(iTrack))->GetOuterHmpParam()) (esd.GetTrack(iTrack))->GetOuterHmpPxPyPz(hmpMom);
2104
2105      if(esd.GetTrack(iTrack)->GetHMPIDsignal() == 0 && thetaTrk == 0 && qMip == 0 && nphMip ==0 ) continue;  //
2106
2107      new(hmpidRings[cntHmpidGoodTracks++]) AliAODHMPIDrings(
2108                                                                  (esd.GetTrack(iTrack))->GetID(),             // Unique track id to attach the ring to
2109                                                                  1000000*nphMip+qMip,                         // MIP charge and number of photons
2110                                                                  (esd.GetTrack(iTrack))->GetHMPIDcluIdx(),    // 1000000*chamber id + cluster idx of the assigned MIP cluster
2111                                                                  thetaTrk,                                    // track inclination angle theta
2112                                                                  phiTrk,                                      // track inclination angle phi
2113                                                                  (esd.GetTrack(iTrack))->GetHMPIDsignal(),    // Cherenkov angle
2114                                                                  (esd.GetTrack(iTrack))->GetHMPIDoccupancy(), // Occupancy claculated for the given chamber
2115                                                                  (esd.GetTrack(iTrack))->GetHMPIDchi2(),      // Ring resolution squared
2116                                                                   xTrk,                                       // Track x coordinate (LORS)
2117                                                                   yTrk,                                       // Track y coordinate (LORS)
2118                                                                   xMip,                                       // MIP x coordinate (LORS)
2119                                                                   yMip,                                       // MIP y coordinate (LORS)
2120                                                                   hmpPid,                                     // PID probablities from ESD, remove later once it is in CombinedPid
2121                                                                   hmpMom                                      // Track momentum in HMPID at ring reconstruction
2122                                                                );
2123
2124       //  Printf(Form("+++++++++ yes/no: %d  %lf %lf %lf %lf %lf %lf ",(esd.GetTrack(iTrack))->IsHMPID(),thetaTrk, (esd.GetTrack(iTrack))->GetHMPIDchi2(),xTrk, yTrk , xMip, yMip));
2125
2126
2127    }// HMPID signal > -20
2128   }//___esd track loop
2129
2130   return cntHmpidGoodTracks;
2131 }
2132
2133 //______________________________________________________________________________
2134 void AliAnalysisTaskESDfilter::ConvertESDtoAOD()
2135 {
2136   // ESD Filter analysis task executed for each event
2137
2138   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
2139
2140   if(!esd)return;
2141
2142   AliCodeTimerAuto("",0);
2143
2144   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );
2145
2146       // Reconstruct cascades and V0 here
2147   if (fIsV0CascadeRecoEnabled) {
2148     esd->ResetCascades();
2149     esd->ResetV0s();
2150
2151     AliV0vertexer lV0vtxer;
2152     AliCascadeVertexer lCascVtxer;
2153
2154     lV0vtxer.SetCuts(fV0Cuts);
2155     lCascVtxer.SetCuts(fCascadeCuts);
2156
2157
2158     lV0vtxer.Tracks2V0vertices(esd);
2159     lCascVtxer.V0sTracks2CascadeVertices(esd);
2160   }
2161
2162
2163   fNumberOfTracks = 0;
2164   fNumberOfPositiveTracks = 0;
2165   fNumberOfV0s = 0;
2166   fNumberOfVertices = 0;
2167   fNumberOfCascades = 0;
2168   fNumberOfKinks = 0;
2169
2170   AliAODHeader* header = ConvertHeader(*esd);
2171
2172   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);
2173   if ( fIsTZEROEnabled ) ConvertTZERO(*esd);
2174
2175   // Fetch Stack for debuggging if available
2176   fMChandler=0x0;
2177   if(MCEvent())
2178   {
2179     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
2180   }
2181
2182   // loop over events and fill them
2183   // Multiplicity information needed by the header (to be revised!)
2184   Int_t nTracks    = esd->GetNumberOfTracks();
2185   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);
2186
2187   // Update the header
2188
2189   Int_t nV0s      = esd->GetNumberOfV0s();
2190   Int_t nCascades = esd->GetNumberOfCascades();
2191   Int_t nKinks    = esd->GetNumberOfKinks();
2192   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
2193   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
2194   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
2195   nVertices+=nPileSPDVertices;
2196   nVertices+=nPileTrkVertices;
2197   Int_t nJets     = 0;
2198   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
2199   Int_t nFmdClus  = 0;
2200   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
2201   Int_t nHmpidRings = 0;
2202
2203   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));
2204
2205   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus,nHmpidRings);
2206
2207   if (nV0s > 0)
2208   {
2209     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0
2210     fAODV0VtxRefs = new TRefArray(nV0s);
2211     // RefArray to store the mapping between esd V0 number and newly created AOD-V0
2212     fAODV0Refs = new TRefArray(nV0s);
2213     // Array to take into account the V0s already added to the AOD (V0 within cascades)
2214     fUsedV0 = new Bool_t[nV0s];
2215     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;
2216   }
2217
2218   if (nTracks>0)
2219   {
2220     // RefArray to store the mapping between esd track number and newly created AOD-Track
2221
2222     fAODTrackRefs = new TRefArray(nTracks);
2223
2224     // Array to take into account the tracks already added to the AOD
2225     fUsedTrack = new Bool_t[nTracks];
2226     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;
2227   }
2228
2229   // Array to take into account the kinks already added to the AOD
2230   if (nKinks>0)
2231   {
2232     fUsedKink = new Bool_t[nKinks];
2233     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;
2234   }
2235
2236   ConvertPrimaryVertices(*esd);
2237
2238   //setting best TOF PID
2239   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
2240   if (esdH)
2241       fESDpid = esdH->GetESDpid();
2242
2243   if (fIsPidOwner && fESDpid){
2244     delete fESDpid;
2245     fESDpid = 0;
2246   }
2247   if(!fESDpid)
2248   { //in case of no Tender attached
2249     fESDpid = new AliESDpid;
2250     fIsPidOwner = kTRUE;
2251   }
2252
2253   if(!esd->GetTOFHeader())
2254   { //protection in case the pass2 LHC10b,c,d have been processed without tender.
2255     Float_t t0spread[10];
2256     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!!
2257     for (Int_t i=0; i<10; i++) t0spread[i] = (TMath::Sqrt(esd->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
2258     fESDpid->GetTOFResponse().SetT0resolution(t0spread);
2259     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
2260       //    fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);
2261     AliTOFHeader tmpTOFHeader(0,t0spread[0],0,NULL,NULL,NULL,intrinsicTOFres,t0spread[0]);
2262     AODEvent()->SetTOFHeader(&tmpTOFHeader);         // write dummy TOF header in AOD
2263   } else {
2264     AODEvent()->SetTOFHeader(esd->GetTOFHeader());    // write TOF header in AOD
2265   }
2266
2267   //  if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender.
2268
2269   if ( fAreCascadesEnabled ) ConvertCascades(*esd);
2270
2271   if ( fAreV0sEnabled ) ConvertV0s(*esd);
2272
2273   if ( fAreKinksEnabled ) ConvertKinks(*esd);
2274
2275   if ( fAreTracksEnabled ) ConvertTracks(*esd);
2276
2277   // Update number of AOD tracks in header at the end of track loop (M.G.)
2278   header->SetRefMultiplicity(fNumberOfTracks);
2279   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);
2280   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);
2281
2282   if ( fTPCConstrainedFilterMask ) ConvertTPCOnlyTracks(*esd);
2283   if( fGlobalConstrainedFilterMask) ConvertGlobalConstrainedTracks(*esd);
2284
2285   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);
2286
2287   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);
2288
2289   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);
2290
2291   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);
2292
2293         if ( fAreEMCALTriggerEnabled )ConvertCaloTrigger(TString("EMCAL"), *esd);
2294
2295         if ( fArePHOSTriggerEnabled )ConvertCaloTrigger(TString("PHOS"), *esd);
2296
2297   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);
2298   if ( fIsZDCEnabled ) ConvertZDC(*esd);
2299
2300   if(fIsHMPIDEnabled) nHmpidRings = ConvertHMPID(*esd);
2301
2302   delete fAODTrackRefs; fAODTrackRefs=0x0;
2303   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;
2304   delete fAODV0Refs; fAODV0Refs=0x0;
2305
2306   delete[] fUsedTrack; fUsedTrack=0x0;
2307   delete[] fUsedV0; fUsedV0=0x0;
2308   delete[] fUsedKink; fUsedKink=0x0;
2309
2310   if ( fIsPidOwner){
2311     delete fESDpid;
2312     fESDpid = 0x0;
2313   }
2314
2315
2316 }
2317
2318
2319 //______________________________________________________________________________
2320 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)
2321 {
2322   //
2323   // Setter for the raw PID detector signals
2324   //
2325
2326   // Save PID object for candidate electrons
2327     Bool_t pidSave = kFALSE;
2328     if (fTrackFilter) {
2329         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");
2330         if (selectInfo)  pidSave = kTRUE;
2331     }
2332
2333
2334     // Tracks passing pt cut
2335     if(esdtrack->Pt()>fHighPthreshold) {
2336         pidSave = kTRUE;
2337     } else {
2338         if(fPtshape){
2339             if(esdtrack->Pt()> fPtshape->GetXmin()){
2340                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
2341                 if(gRandom->Rndm(0)<1./y){
2342                     pidSave = kTRUE;
2343                 }//end rndm
2344             }//end if p < pmin
2345         }//end if p function
2346     }// end else
2347
2348     if (pidSave) {
2349       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track
2350         detpid = new AliAODPid();
2351         SetDetectorRawSignals(detpid,esdtrack);
2352         aodtrack->SetDetPID(detpid);
2353       }
2354     }
2355 }
2356
2357 //______________________________________________________________________________
2358 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)
2359 {
2360 //
2361 //assignment of the detector signals (AliXXXesdPID inspired)
2362 //
2363  if(!track){
2364  AliInfo("no ESD track found. .....exiting");
2365  return;
2366  }
2367  // TPC momentum
2368  const AliExternalTrackParam *in=track->GetInnerParam();
2369  if (in) {
2370    aodpid->SetTPCmomentum(in->GetP());
2371  }else{
2372    aodpid->SetTPCmomentum(-1.);
2373  }
2374
2375
2376  aodpid->SetITSsignal(track->GetITSsignal());
2377  Double_t itsdedx[4]; // dE/dx samples for individual ITS layers
2378  track->GetITSdEdxSamples(itsdedx);
2379  aodpid->SetITSdEdxSamples(itsdedx);
2380
2381  aodpid->SetTPCsignal(track->GetTPCsignal());
2382  aodpid->SetTPCsignalN(track->GetTPCsignalN());
2383  if(track->GetTPCdEdxInfo()) aodpid->SetTPCdEdxInfo(track->GetTPCdEdxInfo());
2384
2385  //n TRD planes = 6
2386  Int_t nslices = track->GetNumberOfTRDslices()*6;
2387  TArrayD trdslices(nslices);
2388  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
2389    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
2390  }
2391
2392 //TRD momentum
2393  for(Int_t iPl=0;iPl<6;iPl++){
2394    Double_t trdmom=track->GetTRDmomentum(iPl);
2395    aodpid->SetTRDmomentum(iPl,trdmom);
2396  }
2397
2398  aodpid->SetTRDslices(track->GetNumberOfTRDslices()*6,trdslices.GetArray());
2399  aodpid->SetTRDsignal(track->GetTRDsignal());
2400
2401  //TRD clusters and tracklets
2402  aodpid->SetTRDncls(track->GetTRDncls());
2403  aodpid->SetTRDntrackletsPID(track->GetTRDntrackletsPID());
2404
2405  aodpid->SetTRDChi2(track->GetTRDchi2());
2406
2407  //TOF PID
2408  Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
2409  aodpid->SetIntegratedTimes(times);
2410
2411    //  Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
2412    //  aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);
2413    aodpid->SetTOFsignal(track->GetTOFsignal());
2414
2415   Double_t tofRes[5];
2416   for (Int_t iMass=0; iMass<5; iMass++){
2417     //    tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));
2418     tofRes[iMass]=0; //backward compatibility
2419   }
2420   aodpid->SetTOFpidResolution(tofRes);
2421
2422 //  aodpid->SetHMPIDsignal(0); // set to zero for compression but it will be removed later
2423
2424 }
2425
2426 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)
2427 {
2428     // Calculate chi2 per ndf for track
2429     Int_t  nClustersTPC = track->GetTPCNcls();
2430
2431     if ( nClustersTPC > 5) {
2432        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));
2433     } else {
2434        return (-1.);
2435     }
2436  }
2437
2438
2439 //______________________________________________________________________________
2440 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
2441 {
2442 // Terminate analysis
2443 //
2444     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
2445 }
2446
2447 //______________________________________________________________________________
2448 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){
2449 // Print MC info
2450   if(!pStack)return;
2451   label = TMath::Abs(label);
2452   TParticle *part = pStack->Particle(label);
2453   Printf("########################");
2454   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
2455   part->Print();
2456   TParticle* mother = part;
2457   Int_t imo = part->GetFirstMother();
2458   Int_t nprim = pStack->GetNprimary();
2459   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
2460   while((imo >= nprim)) {
2461     mother =  pStack->Particle(imo);
2462     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
2463     mother->Print();
2464     imo =  mother->GetFirstMother();
2465   }
2466   Printf("########################");
2467 }
2468
2469 //______________________________________________________
2470