PWGCF/Correlations/JCORRAN/AliJFilter.cxx

   1 // #include <Riostream.h>
   2 // #include <TChain.h>
   3 // #include <TVectorT.h>
   4 // #include <TVector3.h>
   5 // #include <TFile.h>
   6 // #include <TH1.h>
   7 // #include <TClonesArray.h>
   8 // #include <TObjArray.h>
   9 // #include <TObjString.h>
  10 // #include <TFormula.h>
  11 // #include <TString.h>
  12 // #include <TRefArray.h>
  13 // #include <TNtuple.h>
  14 // #include <TArrayF.h>
  15
  16
  17 #include "AliJFilter.h"
  18 #include "AliAnalysisManager.h"
  19 #include "AliAnalysisTaskSE.h"
  20 #include "AliESDEvent.h"
  21 #include "AliMCEvent.h"
  22 #include "AliStack.h"
  23 #include "AliGenEventHeader.h"
  24 #include "AliGenCocktailEventHeader.h"
  25 #include "AliGenPythiaEventHeader.h"
  26 #include "AliInputEventHandler.h"
  27 #include "AliESDCaloCluster.h"
  28 #include "AliAODEvent.h"
  29 #include "AliAODHeader.h"
  30 #include "AliAODHandler.h"
  31 #include "AliLog.h"
  32 #include "AliESDVertex.h"
  33 #include "AliESDtrack.h"
  34 #include "AliAODTrack.h"
  35 #include "AliAnalysisFilter.h"
  36 #include "AliESDtrackCuts.h"
  37 #include "AliAODVertex.h"
  38 #include "AliAODTracklets.h"
  39 #include "AliAODPid.h"
  40 #include "AliAODMCHeader.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliESDUtils.h"
  43 //#include "AliESDVZERO.h"
  44 #include "AliCentrality.h"
  45 #include "AliAODTracklets.h"
  46 #include "AliMultiplicity.h"
  47 #include "AliJConst.h"
  48 #include "AliESDRun.h"
  49 #include "AliDAQ.h"
  50 #include "AliESDVZERO.h"
  51 #include "AliExternalTrackParam.h"
  52 #include "AliHeader.h"
  53 //== EMCAL
  54 #include "AliESDCaloCluster.h"
  55 #include "AliEMCALGeometry.h"
  56 #include "AliVCluster.h"
  57 #include "AliVCaloCells.h"
  58 #include "AliEMCALRecoUtils.h"
  59 #include "AliEMCALPIDUtils.h"
  60
  61 #include "AliJTrack.h"
  62 #include "AliJMCTrack.h"
  63 #include "AliJPhoton.h"
  64 //#include "AliJCaloCell.h"
  65 #include "AliJEventHeader.h"
  66 #include "AliJRunHeader.h"
  67
  68 #include "AliPIDResponse.h"
  69 #include "AliPIDCombined.h"
  70 #include "AliPHOSGeoUtils.h"
  71 #include "AliAnalysisUtils.h"
  72
  73
  74 ClassImp(AliJFilter);
  75
  76 //______________________________________________________________________________
  77 AliJFilter::AliJFilter() :
  78         TNamed(),
  79         fEsdTrackCuts(0x0),
  80         fESDFilter(0x0),
  81         fIsRealOrMC(0),
  82         fStoreEventPlaneSource(0),
  83         fOADBPath(),
  84         fCaloClustersArr(0),
  85         fClusterThreshold(0),
  86         fTrackThreshold(0),
  87         fEventSuccess(0),
  88         fMcMap(0),
  89         fTrackList(0),
  90         fMCTrackList(0x0),
  91         fPhotonList(0x0),
  92         fCaloCellList(0x0),
  93         fHeaderList(0x0),
  94         fRunInfoList(0x0),
  95         fPIDResponse(0x0),
  96         fPIDCombined(0x0),
  97         fVZEROData(0x0),
  98         fTZEROData(0x0),
  99         //fFMDData(0x0),
 100         fZDCData(0x0),
 101         fEMCLabels(0),
 102         fEMCTreeLabels(0),
 103         fAliJRunHeader(0x0),
 104         fEMCALGeometry(0x0),
 105         fEMCALRecoUtils(0x0),
 106         fPHOSGeom(0x0),
 107         fAnaUtils(0x0),
 108         fMyTask(0x0)
 109 {
 110         //Default constructor
 111 }
 112
 113 //______________________________________________________________________________
 114 AliJFilter::AliJFilter(const char *name,AliAnalysisTaskSE *task):
 115         TNamed(name,name),
 116         fEsdTrackCuts(0x0),
 117         fESDFilter(0x0),
 118         fIsRealOrMC(0),
 119         fStoreEventPlaneSource(0),
 120         fOADBPath(),
 121         fCaloClustersArr(0),
 122         fClusterThreshold(0),
 123         fTrackThreshold(0),
 124         fEventSuccess(0),
 125         fMcMap(0),
 126         fTrackList(0),
 127         fMCTrackList(0x0),
 128         fPhotonList(0x0),
 129         fCaloCellList(0x0),
 130         fHeaderList(0x0),
 131         fRunInfoList(0x0),
 132         fPIDResponse(0x0),
 133         fPIDCombined(0x0),
 134         fVZEROData(0x0),
 135         fTZEROData(0x0),
 136         //fFMDData(0x0),
 137         fZDCData(0x0),
 138         fEMCLabels(0),
 139         fEMCTreeLabels(0),
 140         fAliJRunHeader(0x0),
 141         fEMCALGeometry(0x0),
 142         fEMCALRecoUtils(0x0),
 143         fPHOSGeom(0x0),
 144         fAnaUtils(0x0),
 145         fMyTask(0x0)
 146 {
 147         // Constructor
 148         if(task->DebugLevel() > 5) cout << "---- AliJFilter Constructor ----"<<endl;
 149
 150 }
 151
 152 //____________________________________________________________________________
 153 AliJFilter::AliJFilter(const AliJFilter& ap) :
 154         TNamed(ap.GetName(), ap.GetTitle()),
 155         fEsdTrackCuts(ap.fEsdTrackCuts),
 156         fESDFilter(ap.fESDFilter),
 157         fIsRealOrMC(ap.fIsRealOrMC),
 158         fStoreEventPlaneSource(ap.fStoreEventPlaneSource),
 159         fOADBPath(ap.fOADBPath),
 160         fCaloClustersArr(ap.fCaloClustersArr),
 161         fClusterThreshold(ap.fClusterThreshold),
 162         fTrackThreshold(ap.fTrackThreshold),
 163         fEventSuccess(ap.fEventSuccess),
 164         fMcMap(ap.fMcMap),
 165         fTrackList(ap.fTrackList),
 166         fMCTrackList(ap.fMCTrackList),
 167         fPhotonList(ap.fPhotonList),
 168         fCaloCellList(ap.fCaloCellList),
 169         fHeaderList(ap.fHeaderList),
 170         fRunInfoList(ap.fRunInfoList),
 171         fPIDResponse(ap.fPIDResponse),
 172         fPIDCombined(ap.fPIDCombined),
 173         fVZEROData(ap.fVZEROData),
 174         fTZEROData(ap.fTZEROData),
 175         //fFMDData(ap.fFMDData),
 176         fZDCData(ap.fZDCData),
 177         fEMCLabels(ap.fEMCLabels),
 178         fEMCTreeLabels(ap.fEMCTreeLabels),
 179         fAliJRunHeader(ap.fAliJRunHeader),
 180         fEMCALGeometry(ap.fEMCALGeometry),
 181         fEMCALRecoUtils(ap.fEMCALRecoUtils),
 182         fPHOSGeom(ap.fPHOSGeom),
 183         fAnaUtils(ap.fAnaUtils),
 184         fMyTask(ap.fMyTask)
 185 {
 186         // cpy ctor
 187 }
 188
 189 //_____________________________________________________________________________
 190 AliJFilter& AliJFilter::operator = (const AliJFilter& ap)
 191 {
 192         // assignment operator
 193
 194         this->~AliJFilter();
 195         new(this) AliJFilter(ap);
 196         return *this;
 197 }
 198
 199 //______________________________________________________________________________
 200 AliJFilter::~AliJFilter()
 201 {
 202         // destructor
 203         delete fMcMap;
 204         delete fTrackList;
 205         delete fMCTrackList;
 206         delete fPhotonList;
 207         delete fCaloCellList;
 208         delete fHeaderList;
 209         delete fAliJRunHeader;
 210         delete fRunInfoList;
 211         delete fPIDResponse;
 212         delete fPIDCombined;
 213         delete fEMCALRecoUtils;
 214         delete fEMCALGeometry;
 215         delete fPHOSGeom;
 216         delete fAnaUtils;
 217         delete fVZEROData;
 218         delete fTZEROData;
 219         delete fZDCData;
 220         //  delete fFMDData;
 221
 222
 223 }
 224
 225 //________________________________________________________________________
 226
 227 void AliJFilter::UserCreateOutputObjects()
 228 {
 229         //=== create the jcorran outputs objects
 230         if(fMyTask->DebugLevel() > 1) printf("AliJFilter::UserCreateOutPutData() \n");
 231
 232         //== RUN HEADER
 233         cout<<"TEST2 "<<fAliJRunHeader<<endl;
 234         if(!fAliJRunHeader) fAliJRunHeader = new AliJRunHeader();
 235         fRunInfoList  = new TList();
 236         fRunInfoList->SetName("RunInfoList");
 237         fRunInfoList->SetOwner();
 238         fRunInfoList->Clear();
 239         fRunInfoList->Add(fAliJRunHeader);
 240
 241         //=== Other Objects
 242         fCaloClustersArr = new TRefArray();
 243         fEMCALGeometry = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 244         fEMCALRecoUtils = new AliEMCALRecoUtils();
 245         fPHOSGeom = new AliPHOSGeoUtils();
 246         fAnaUtils = new AliAnalysisUtils();
 247         fAnaUtils->SetUseOutOfBunchPileUp( kTRUE );
 248         fMcMap = new TArrayI();
 249
 250         //=== Set Tree and TClonesArray
 251         //== TRACKS
 252         AddList("AliJTrackList", "AliJTrack", &fTrackList, 1000);
 253         if( fAliJRunHeader->GetStoreEMCalInfo() ){
 254                 AddList("AliJPhotonList", "AliJPhoton", &fPhotonList, 1000);
 255                 //BS AddList("AliJCaloCell", "AliJCaloCell", &fCaloCellList, 1000);
 256         }
 257         if( IsMC() )
 258                 AddList("AliJMCTrackList", "AliJMCTrack", &fMCTrackList, 1000);
 259         //== Event Header
 260         AddList("AliJEventHeaderList", "AliJEventHeader", &fHeaderList, 1000);
 261
 262         //== EventPlane SRC
 263         if( fAliJRunHeader->GetStoreEventPlaneSource() ){
 264                 fVZEROData = new AliESDVZERO;
 265                 fTZEROData = new AliESDTZERO;
 266                 fZDCData   = new AliESDZDC;
 267         }
 268         //== PID
 269         //   fPIDCombined = new AliPIDCombined;
 270         //   fPIDCombined->SetDefaultTPCPriors();
 271         //   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
 272         //   fPIDResponse = ((AliInputEventHandler*) (man->GetInputEventHandler()))->GetPIDResponse();
 273         //   fPIDResponse->SetOADBPath(AliAnalysisManager::GetOADBPath());
 274         //   if (!fOADBPath.IsNull()) fPIDResponse->SetOADBPath(fOADBPath.Data());
 275
 276         cout << "Add(fAliJRunHeader) in UserCreateObject() ======= " << endl;
 277
 278 }
 279
 280 //______________________________________________________________________________
 281 void AliJFilter::UserExec(Option_t* /*option*/)
 282 {
 283         // user loop
 284         AliJRunHeader *runh = fAliJRunHeader;
 285         Bool_t hasGoodTrack, hasGoodCluster;
 286
 287         fEventSuccess = kFALSE;
 288
 289         // Processing of one event
 290         DEBUG( 5, 1, "------- AliJFilter Exec-------" );
 291         if(!((fMyTask->Entry()-1)%100))  AliInfo(Form(" Processing event # %lld",  fMyTask->Entry()));
 292
 293         //=== Init Variables
 294         fTrackList->Clear();
 295         if( IsMC() ){
 296                 fMCTrackList->Clear();
 297                 fEMCLabels.clear();
 298                 fEMCTreeLabels.clear();
 299         }
 300
 301         if( fAliJRunHeader->GetStoreEMCalInfo() ){
 302                 fPhotonList->Clear("C");
 303                 fCaloCellList->Clear();
 304         }
 305         fHeaderList->Clear();
 306
 307         hasGoodCluster = kTRUE;
 308         hasGoodTrack = kTRUE;
 309
 310         //=== CHECK ESD, AOD, MC event
 311         if( !Event() ) return;
 312
 313         if( FromESD() ) {   //Reading ESD
 314                 DEBUG( 5, 1, "\t------- Start ESD " );
 315                 if( !ESDEvent() ) return;
 316                 if( runh->GetWithoutSDD() && !(ESDEvent()->GetTriggerMask()  & (1<<13)) ) return;
 317
 318                 if( IsMC() ){
 319                         if( ! MCEvent() ) return;
 320                 }
 321         }
 322
 323         if( FromAOD() ) {
 324                 DEBUG( 5, 1, "\t------- Start AOD " );
 325                 if( !AODEvent() ) return;
 326         }
 327
 328
 329         // pileup rejection
 330         if( fAnaUtils->IsPileUpEvent( Event() ))
 331                 return;
 332
 333         //---------------------------------------------------------------
 334         // RUN Header
 335         //---------------------------------------------------------------
 336         if(!runh->GetRunNumber()){ //new run has started : I suppose no change of run in process
 337                 runh->SetRunNumber( Event()->GetRunNumber() );
 338                 if( FromESD() ){
 339                         //==== General ====//
 340                         runh->SetBeamEnergy( ESDEvent()->GetBeamEnergy() );
 341                         runh->SetBeamType( ESDEvent()->GetBeamType() );
 342                         //==== Detector status ==//
 343                         if( ESDEvent()->GetCurrentL3() > 0 ) runh->SetL3MagnetFieldPolarity(1);
 344                         if( ESDEvent()->GetCurrentL3() < 0 ) runh->SetL3MagnetFieldPolarity(-1);
 345                         runh->SetL3MagnetFieldIntensity( ESDEvent()->GetMagneticField() );
 346                         runh->SetCurrentL3( ESDEvent()->GetCurrentL3() );
 347                         runh->SetCurrentDip( ESDEvent()->GetCurrentDip() );
 348                         runh->SetUniformBMap( ESDEvent()->IsUniformBMap() );
 349                         //==== Triggers ====//
 350                         const AliESDRun* esdRun = ESDEvent()->GetESDRun();
 351                         for(Int_t triggerBit=0; triggerBit<kRangeTriggerTableAlice; triggerBit++){
 352                                 runh->SetActiveTriggersAlice( triggerBit, esdRun->GetTriggerClass(triggerBit) );
 353                         }
 354                 }
 355                 else if( FromAOD() ){
 356                         //==== General ====//
 357                         cout << "Run # = "<< AODEvent()->GetRunNumber() << endl;
 358                         runh->SetRunNumber( AODEvent()->GetRunNumber() );
 359                         //TODO runh->SetBeamEnergy( ESDEvent()->GetBeamEnergy() );
 360                         //TODO runh->SetBeamType( ESDEvent()->GetBeamType() );
 361                         //==== Detector status ==//
 362                         //TODO runh->Setl3MgFieldPolarity(1);
 363                         runh->SetL3MagnetFieldIntensity( AODEvent()->GetMagneticField() );
 364                         runh->SetCurrentL3( AODEvent()->GetMagneticField()*30000.0/5.00668 );
 365                         runh->SetCurrentDip( AODEvent()->GetMuonMagFieldScale()*6000.0 );
 366                         runh->SetUniformBMap( kFALSE ); // TODO is this?
 367                 }
 368                 cout << "Add(fAliJRunHeader) is done =============" << endl;
 369         }
 370
 371         //---------------------------------------------------------------
 372         // EventHeader and read Others
 373         //---------------------------------------------------------------
 374         if( FromESD() ){   //Reading ESD
 375                 DEBUG( 5, 1, "\t------- Start READ ESD " );
 376
 377                 ReadESDHeader( ESDEvent() );
 378                 ReadESDTracks( ESDEvent() );
 379
 380                 if( fAliJRunHeader->GetStoreEMCalInfo() ){
 381                         ReadESDCaloClusters( ESDEvent() );
 382                         ReadESDCaloCells( ESDEvent() );
 383                 }
 384                 if( IsMC() ){
 385                         ReadMCTracksFromESD();
 386                         //RemapMCLabels();
 387                 }
 388         }else if( FromAOD() ){
 389                 DEBUG( 5, 1, "\t------- Start READ AOD " );
 390                 ReadAODHeader( AODEvent() );
 391                 ReadAODTracks( AODEvent() );
 392                 if( fAliJRunHeader->GetStoreEMCalInfo() ){
 393                         ReadAODCaloClusters( AODEvent() );
 394                         ReadAODCaloCells( AODEvent() );
 395                 }
 396                 if( IsMC() ){
 397                         ReadMCTracksFromAOD();
 398                         //RemapMCLabels();
 399                 }
 400
 401         }else{
 402                 cout << "Error: Not correct InputDataFormat especified " << endl;
 403                 return;
 404         }
 405
 406         if( hasGoodCluster || hasGoodTrack ){
 407                 //=== TODO : need this?
 408                 AliAODHandler* outputHandler =
 409                         (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 410                 outputHandler->SetFillAOD(kTRUE);
 411                 outputHandler->SetFillExtension(kTRUE);
 412                 fEventSuccess = kTRUE;
 413         }
 414         else{
 415                 fTrackList->Clear();
 416                 if( IsMC() ){
 417                         fMCTrackList->Clear();
 418                         fEMCLabels.clear();
 419                         fEMCTreeLabels.clear();
 420                 }
 421
 422                 if( fAliJRunHeader->GetStoreEMCalInfo() ){
 423                         fPhotonList->Clear("C");
 424                         fCaloCellList->Clear();
 425                 }
 426                 fHeaderList->Clear();
 427         }
 428
 429         DEBUG( 5, 1, "\t------- End UserExec " );
 430 }
 431
 432 //______________________________________________________________________________
 433 void AliJFilter::Init()
 434 {
 435         // Intialisation of parameters
 436         AliInfo("Doing initialization") ;
 437
 438         //   TString formula(fEsdTrackCuts->GetMaxDCAToVertexXYPtDep());
 439         //   if(formula.Length()>0){ // momentum dep DCA cut for AOD
 440         //     formula.ReplaceAll("pt","x");
 441         //   }
 442 }
 443
 444 //______________________________________________________________________________
 445 void AliJFilter::Terminate(Option_t *)
 446 {
 447         // termination
 448         fTrackList->Clear();
 449         if( IsMC() ) fMCTrackList->Clear();
 450         if( fAliJRunHeader->GetStoreEMCalInfo() ){
 451                 fPhotonList->Clear();
 452                 fCaloCellList->Clear();
 453         }
 454         fHeaderList->Clear();
 455
 456         // Processing when the event loop is ended
 457         cout<<"PWG4JCORRAN Analysis DONE !!"<<endl;
 458
 459 }
 460
 461 //______________________________________________________________________________
 462 void AliJFilter::ReadESDTracks(AliESDEvent * esd)
 463         //void AliJFilter::ReadESDTracks(const AliESDEvent * esd)
 464 {
 465         // Read the AliESDtrack and fill the list of AliJTrack containers
 466         Int_t nt = esd->GetNumberOfTracks();
 467         DEBUG( 5, 1 , Form("ESD::NumberOfTracks = %d",nt), "AliJFilter::ReadESDTracks" );
 468
 469         //==== Prepare TPC, GCG track ====//
 470         Float_t ptMaxTPC = 0;
 471         Float_t ptMinTPC = 1E10;
 472         Float_t ptMaxGCG = 0;
 473         Float_t ptMinGCG = 1E10;
 474         for(int i = 0;i<32;i++){
 475                 AliESDtrackCuts* cuts = (AliESDtrackCuts*)fESDFilter->GetCuts()->At(i);
 476                 if(!cuts) continue;
 477                 Float_t tmp1= 0,tmp2 = 0;
 478                 cuts->GetPtRange(tmp1,tmp2);
 479                 if( TESTBIT ( fAliJRunHeader->GetStoreTPCTrackBitMask(), i ) ){
 480                         if(tmp1<ptMinTPC)ptMinTPC=tmp1;
 481                         if(tmp2>ptMaxTPC)ptMaxTPC=tmp2;
 482                 }
 483                 if( TESTBIT(fAliJRunHeader->GetStoreGCGTrackBitMask() , i ) ){
 484                         if(tmp1<ptMinGCG)ptMinGCG=tmp1;
 485                         if(tmp2>ptMaxGCG)ptMaxGCG=tmp2;
 486                 }
 487         }
 488
 489         //==== loop over tracks ====//
 490         for(Int_t it = 0; it < nt; it++) {
 491
 492                 AliESDtrack *track = esd->GetTrack(it);
 493                 if( !track ) continue;
 494                 UInt_t filterMap = fESDFilter->IsSelected( track );
 495                 if(! filterMap ) continue; // apply track selection criteria
 496
 497                 //====create a new AliJTrack and fill the track info
 498                 AliJTrack * ctrack = new( (*fTrackList)[fTrackList->GetEntriesFast()] ) AliJTrack;
 499                 ctrack->SetPxPyPzE(track->Px(), track->Py(), track->Pz(), 0 );
 500                 Double32_t pos[3];
 501                 track->GetXYZ(pos);
 502                 ctrack->SetTrackPos( pos );
 503                 ctrack->SetTPCdEdx( track->GetTPCsignal()  );
 504                 ctrack->SetParticleType(kJNone);
 505                 ctrack->SetCharge(track->Charge());
 506                 ctrack->SetFilterMap( filterMap );
 507                 ctrack->SetLabel( track->GetLabel() );
 508
 509                 ReadESDPID( track, ctrack );
 510                 //==== TPC Tracks ====//
 511                 if( filterMap & fAliJRunHeader->GetStoreTPCTrackBitMask() ) {
 512                         ConvertESDTPCOnlyTracks( esd, it, ctrack, ptMinTPC, ptMaxTPC );
 513                 }
 514                 //==== GCG Tracks ====//
 515                 if( filterMap & fAliJRunHeader->GetStoreGCGTrackBitMask() ) {
 516                         ConvertESDGCGTracks( esd, it, ctrack, ptMinGCG, ptMaxGCG );
 517                 }
 518
 519                 Float_t b[2];
 520                 Float_t bCov[3];
 521                 track->GetImpactParameters(b,bCov);
 522                 //         ctrack->SetDCAtoVertexXY( b[0] );
 523                 //         ctrack->SetDCAtoVertexZ( b[1] );
 524
 525                 //if( track->P()>1 ) DEBUG( 5, 1, Form("P = %f", track->P() ) ) ;
 526
 527         } // end tracks loop
 528 }
 529
 530 //______________________________________________________________________________
 531 void AliJFilter::ConvertESDTPCOnlyTracks(AliESDEvent* esd, int iTrack, AliJTrack * ctrack, double ptmin, double ptmax)
 532 {
 533
 534         const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 535
 536         Double_t pos[3] = { 0. };
 537         Double_t covTr[21]={0.};
 538         //Double_t pid[10]={0.};
 539
 540         Double_t p[3] = { 0. };
 541
 542         Double_t pDCA[3] = { 0. }; // momentum at DCA
 543         Double_t rDCA[3] = { 0. }; // position at DCA
 544         Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 545         Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 546
 547
 548         AliESDtrack* esdTrack = esd->GetTrack(iTrack); //carefull do not modify it othwise  need to work with a copy
 549
 550         // Track selection
 551
 552         AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
 553         if(!track) return;
 554
 555         if(track->Pt()>0.)
 556         {
 557                 // only constrain tracks above threshold
 558                 AliExternalTrackParam exParam;
 559                 // take the B-field from the ESD, no 3D fieldMap available at this point
 560                 Bool_t relate = false;
 561                 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 562                 if(!relate){
 563                         delete track;
 564                         return;
 565                 }
 566                 // fetch the track parameters at the DCA (unconstraint)
 567                 if(track->GetTPCInnerParam()){
 568                         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
 569                         track->GetTPCInnerParam()->GetXYZ(rDCA);
 570                 }
 571                 // get the DCA to the vertex:
 572                 track->GetImpactParametersTPC(dDCA,cDCA);
 573                 // set the constrained parameters to the track
 574                 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 575         }
 576
 577         track->GetPxPyPz(p);
 578
 579         double p2[3];
 580         esdTrack->GetInnerPxPyPz(p2);
 581
 582         Float_t pT = track->Pt();
 583         if(pT<ptmin||pT>ptmax){
 584                 delete track;
 585                 return;
 586         }
 587
 588         track->GetXYZ(pos);
 589         track->GetCovarianceXYZPxPyPz(covTr);
 590
 591         ctrack->SetTPCTrack(p[0], p[1], p[2]);
 592
 593         delete track;
 594 }
 595
 596
 597 void AliJFilter::ConvertESDGCGTracks(AliESDEvent *esd, int iTrack, AliJTrack *ctrack, double ptMin, double ptMax)
 598 {
 599
 600         Double_t pos[3] = { 0. };
 601         Double_t covTr[21]={0.};
 602         Double_t p[3] = { 0. };
 603
 604         Double_t pDCA[3] = { 0. }; // momentum at DCA
 605         Double_t rDCA[3] = { 0. }; // position at DCA
 606         Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 607         Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 608
 609
 610         AliESDtrack* esdTrack = esd->GetTrack(iTrack); //carefull do not modify it othwise  need to work with a copy
 611         const AliExternalTrackParam * exParamGC = esdTrack->GetConstrainedParam();
 612         if(!exParamGC) return;
 613
 614         // fetch the track parameters at the DCA (unconstrained)
 615         esdTrack->GetPxPyPz(pDCA);
 616         esdTrack->GetXYZ(rDCA);
 617         // get the DCA to the vertex:
 618         esdTrack->GetImpactParameters(dDCA,cDCA);
 619
 620         if (!esdTrack->GetConstrainedPxPyPz(p)) return;
 621
 622
 623         Float_t pT = exParamGC->Pt();
 624         if(pT<ptMin||pT>ptMax){
 625                 return;
 626         }
 627
 628         esdTrack->GetConstrainedXYZ(pos);
 629         exParamGC->GetCovarianceXYZPxPyPz(covTr);
 630
 631         ctrack->SetGCGTrack(p[0], p[1], p[2]);
 632 }
 633
 634
 635
 636 //_________________________________________________________________________________-
 637 void AliJFilter::ReadESDPID(AliESDtrack *track, AliJTrack *ctrack)
 638 {
 639         // To reduce the size of output, the variables which cannot be calculated later are only kept
 640         // expected TOF signal, TPC momentum for expected TPC signal. Measured values are stored in ReadESDTrack()
 641         // 1. expected TOF signal
 642         Double_t times[AliPID::kSPECIES];
 643         track->GetIntegratedTimes(times);
 644         for(int ip=0; ip < (AliJTrack::kNAliJTrkPID); ip++) {
 645                 ctrack->SetExpectedTOFsignal(AliJTrack::AliJTrkPID(ip), times[ip]);
 646
 647         }
 648         // 2. TPC momentum
 649         Double_t momTPC = track->GetTPCmomentum();
 650         ctrack->SetTPCmomentum(momTPC);
 651 }
 652
 653 //______________________________________________________________________________
 654 Bool_t AliJFilter::ReadAODTracks(const AliAODEvent * aod)
 655 {
 656         // AOD track reader
 657         Bool_t hasGoodTrack;
 658         hasGoodTrack = kFALSE;
 659
 660         // Read the AliAODtrack and fill the list of AliJTrack containers
 661         Int_t nt = aod->GetNumberOfTracks();
 662         Int_t listnt = 0;
 663
 664         DEBUG(5, 1, Form("AOD::NumberOfTracks = %d",nt) );
 665
 666         //==== loop over tracks ====//
 667         for(Int_t it = 0; it < nt; it++) {
 668
 669                 AliAODTrack *track = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
 670                 if(!track) AliFatal("Not a standard AOD track");
 671                 //if(track->GetFilterMap() & (1 << 7) ) continue;
 672                 //if(!AcceptAODTrack(track)) continue;
 673                 //if(! fEsdTrackCuts->IsSelected(track)) continue; //apply loose selection criteria
 674                 //FK//if(track->GetType() != AliAODTrack::kPrimary) continue; // only primaries
 675                 //
 676
 677                 AliJTrack * ctrack = new( (*fTrackList)[listnt++] ) AliJTrack;
 678                 ctrack->SetID( track->GetID() );
 679                 ctrack->SetPxPyPzE(track->Px(), track->Py(), track->Pz(), 0 );
 680                 Double32_t pos[3];
 681                 track->GetXYZ(pos);
 682                 ctrack->SetTrackPos( pos );
 683                 //TODO if( fStoreTPCTrack )
 684                 ctrack->SetParticleType(kJNone);
 685                 ctrack->SetCharge(track->Charge());
 686                 ctrack->SetStatus(track->GetStatus());//
 687                 ctrack->SetFlags( track->GetFlags() );
 688                 ctrack->SetLabel( track->GetLabel() );
 689                 //     //FilterMap
 690                 //     UInt_t filterMap=0;
 691                 //     for( unsigned int i=0;i<sizeof(filterMap)*8;i++ ){
 692                 //       if( track->TestFilterBit( BIT(i) )){
 693                 //         SETBIT( filterMap ,  i);
 694                 //       }
 695                 //     }
 696                 ctrack->SetFilterMap( track->GetFilterMap() );
 697
 698                 //PID TODO
 699                 double const * pid = track->PID();
 700                 ctrack->SetPID(AliJTrack::kElectronAliJ,pid[AliAODTrack::kElectron],AliJTrack::kTOF);
 701                 ctrack->SetPID(AliJTrack::kMuonAliJ,    pid[AliAODTrack::kMuon],    AliJTrack::kTOF);
 702                 ctrack->SetPID(AliJTrack::kPionAliJ,    pid[AliAODTrack::kPion],    AliJTrack::kTOF);
 703                 ctrack->SetPID(AliJTrack::kKaonAliJ,    pid[AliAODTrack::kKaon],    AliJTrack::kTOF);
 704                 ctrack->SetPID(AliJTrack::kProtonAliJ,  pid[AliAODTrack::kProton],  AliJTrack::kTOF);
 705                 //TPC
 706                 ctrack->SetTPCnClust(track->GetTPCNcls());
 707                 ctrack->SetTPCdEdx( track->GetTPCsignal()  );
 708                 ctrack->SetTOFsignal( track->GetTOFsignal() );
 709                 ctrack->SetLabel( track->GetLabel() );
 710                 for( int i=0;i<int(sizeof(UInt_t)*8);i++ ){
 711                         ctrack->SetBit( i, track->TestBit( i ));
 712                 }
 713
 714                 // check track threshold
 715                 if( track->Pt() > fTrackThreshold )
 716                         hasGoodTrack = kTRUE;
 717
 718                 //if(fMyTask->DebugLevel() > 5 && track->P()>1 ) cout << "P = " << track->P() << endl;
 719         } // end tracks loop
 720
 721         return hasGoodTrack;
 722 }
 723
 724
 725 //______________________________________________________________________________
 726 AliJEventHeader* AliJFilter::ReadCommonHeader(AliVEvent *event){
 727         //Read the AliVEvent and fill the list of AliJEventHeader containers
 728         //create a header and fill it
 729         AliJEventHeader *hdr = new( (*fHeaderList)[fHeaderList->GetEntriesFast()] ) AliJEventHeader;
 730
 731
 732         // Get Centrality as a percent from 0% to 100%
 733         AliCentrality *cent = event->GetCentrality();
 734         if( cent ){
 735                 hdr->SetCentrality( cent->GetCentralityPercentile("V0M"));
 736                 hdr->SetCentralityArray(AliJEventHeader::kcV0M, cent->GetCentralityPercentile("V0M"));
 737                 hdr->SetCentralityArray(AliJEventHeader::kcFMD, cent->GetCentralityPercentile("FMD"));
 738                 hdr->SetCentralityArray(AliJEventHeader::kcTRK, cent->GetCentralityPercentile("TRK"));
 739                 hdr->SetCentralityArray(AliJEventHeader::kcTKL, cent->GetCentralityPercentile("TKL"));
 740                 hdr->SetCentralityArray(AliJEventHeader::kcCL0, cent->GetCentralityPercentile("CL0"));
 741                 hdr->SetCentralityArray(AliJEventHeader::kcCL1, cent->GetCentralityPercentile("CL1"));
 742                 hdr->SetCentralityArray(AliJEventHeader::kcV0MvsFMD, cent->GetCentralityPercentile("V0MvsFMD"));
 743                 hdr->SetCentralityArray(AliJEventHeader::kcTKLvsV0, cent->GetCentralityPercentile("TKLvsV0"));
 744                 hdr->SetCentralityArray(AliJEventHeader::kcZEMvsZDC, cent->GetCentralityPercentile("ZEMvsZDC"));
 745                 hdr->SetCentralityArray(AliJEventHeader::kcV0A, cent->GetCentralityPercentile("V0A"));
 746                 hdr->SetCentralityArray(AliJEventHeader::kcV0C, cent->GetCentralityPercentile("V0C"));
 747         }
 748         hdr->SetTriggerMaskAlice(event->GetTriggerMask()); //ULong64_t
 749         hdr->SetTriggerMaskJCorran(ConvertTriggerMask()); //UInt_t
 750         hdr->SetEventType(event->GetEventType());
 751         hdr->SetBunchCrossNumber(event->GetBunchCrossNumber());
 752
 753         int ncontributors = 0;
 754         const AliVVertex * vtxESD = event->GetPrimaryVertex();
 755         if(vtxESD){
 756                 hdr->SetXVertex(vtxESD->GetX()); //FK// EFF
 757                 hdr->SetYVertex(vtxESD->GetY()); //FK// EFF
 758                 hdr->SetZVertex(vtxESD->GetZ());
 759                 //hdr->SetZVertexErr(vtxESD->GetZRes());
 760                 double covMat[6];
 761                 vtxESD->GetCovarianceMatrix(covMat);
 762                 hdr->SetZVertexErr(TMath::Sqrt(covMat[5])); // GetZRes := TMath::Sqrt(fCovZZ)
 763                 ncontributors = vtxESD->GetNContributors(); // get number of contributors to vertex
 764                 hdr->SetVtxMult( vtxESD->GetNContributors() );
 765         }else{
 766                 hdr->SetZVertex(9999);
 767                 hdr->SetZVertexErr(9999);
 768         }
 769         hdr->SetVtxMult(ncontributors); //FK// EFF contrib to vertex
 770         return hdr;
 771 }
 772 //______________________________________________________________________________
 773 void AliJFilter::ReadESDHeader(AliESDEvent *esd)
 774 {
 775         // Read the AliESDEvent and fill the list of AliJEventHeader containers
 776         if(!esd) return;
 777         if( fAliJRunHeader->GetRefitESDVertexTracks() )
 778                 AliESDUtils::RefitESDVertexTracks( esd ); // TODO only for LHC11a right?
 779         AliJEventHeader *hdr = ReadCommonHeader( esd );
 780         //   AliMultiplicity *fSPDMult =(AliMultiplicity *) esd->GetMultiplicity();
 781         //   if(fSPDMult) hdr->SetSPDTrackletMult(fSPDMult->GetNumberOfTracklets());
 782         // This is moved from ReadCommonHeader. AOD should have same.TODO!!
 783         AliESDVZERO *v0 = esd->GetVZEROData();
 784
 785         if( v0 ) hdr->SetV0Mult(v0->GetMTotV0A() + v0->GetMTotV0C());
 786         if( v0 ) hdr->SetV0AMult(v0->GetMTotV0A());
 787         if( v0 ) hdr->SetV0CMult(v0->GetMTotV0C());
 788
 789         const AliESDRun* esdRun = esd->GetESDRun();
 790         //cout <<"========================"<<endl;
 791         //cout << (esdRun->GetDetectorsInReco() & AliDAQ::kSPD) << endl;
 792         //cout << (esdRun->GetDetectorsInReco() & AliDAQ::kSSD) << endl;
 793         //cout << (esdRun->GetDetectorsInReco() & AliDAQ::kSDD) << endl;
 794         //cout << (esdRun->GetDetectorsInReco() & AliDAQ::kTPC) << endl;
 795         if(esdRun->GetDetectorsInReco() & AliDAQ::kSPD) hdr->SetSPDTrackletMult(AliESDtrackCuts::GetReferenceMultiplicity( esd, AliESDtrackCuts::kTracklets, 1.0 ));
 796         if((esdRun->GetDetectorsInReco() & AliDAQ::kSSD) || (esdRun->GetDetectorsInReco() & AliDAQ::kSDD)) hdr->SetITSSATrackletMult(AliESDtrackCuts::GetReferenceMultiplicity( esd, AliESDtrackCuts::kTrackletsITSSA, 1.0 ));
 797         if(esdRun->GetDetectorsInReco() & AliDAQ::kTPC) hdr->SetITSTPCTrackletMult(AliESDtrackCuts::GetReferenceMultiplicity( esd, AliESDtrackCuts::kTrackletsITSTPC, 1.0 ));
 798
 799
 800         //TODO  Store Detector data
 801         if( fAliJRunHeader->GetStoreEventPlaneSource() ){
 802                 *fVZEROData = *esd->GetVZEROData();
 803                 *fTZEROData = AliESDTZERO(*esd->GetESDTZERO());
 804                 *fZDCData  = *esd->GetESDZDC();
 805         }
 806         hdr->SetEventID( esd->GetEventNumberInFile());
 807         //const AliESDVertex * vtxESD = esd->GetPrimaryVertex();
 808         //if( vtxESD->GetStatus() == 0 ) hdr->SetVtxMult( 0 );
 809         // if fNcontributes > 0 then status is always true. do we need this?
 810
 811         //==== MC ====/
 812         if( IsMC() ){
 813                 const AliVVertex * primaryMCVertex = MCEvent()->GetPrimaryVertex();
 814                 //cout<<"AliMCEvent = "<<MCEvent()<<endl;
 815                 //cout<<"AliVVertex = "<<primaryMCVertex<<endl;
 816                 if( primaryMCVertex ){
 817                         hdr->SetXVertexMC( primaryMCVertex->GetX() );
 818                         hdr->SetYVertexMC( primaryMCVertex->GetY() );
 819                         hdr->SetZVertexMC( primaryMCVertex->GetZ() );
 820                 }
 821                 AliESDHeader * esdHeader = esd->GetHeader();
 822                 hdr->SetL0TriggerInputs( esdHeader->GetL0TriggerInputs() );
 823         }
 824 }
 825
 826 //______________________________________________________________________________
 827 void AliJFilter::ReadAODHeader(AliAODEvent *aod)
 828 {
 829         //Read the AliAODEvent and fill the list of AliJEventHeader containers
 830         AliJEventHeader *hdr = ReadCommonHeader( aod );
 831
 832         const AliAODTracklets *trackletsSPD = aod->GetTracklets();
 833         if(trackletsSPD){
 834                 hdr->SetSPDTrackletMult(trackletsSPD->GetNumberOfTracklets());
 835         }
 836
 837         hdr->SetFiredTriggers( aod->GetFiredTriggerClasses() );
 838         //TODO hdr->SetEventID( esd->GetEventNumberInFile());
 839         //==== MC ====//
 840         if( IsMC() ){
 841                 AliAODMCHeader *aodMCheader = (AliAODMCHeader *) aod->FindListObject(AliAODMCHeader::StdBranchName());
 842                 hdr->SetXVertexMC( aodMCheader->GetVtxX() );
 843                 hdr->SetYVertexMC( aodMCheader->GetVtxY() );
 844                 hdr->SetZVertexMC( aodMCheader->GetVtxZ() );
 845         }
 846
 847         AliAODHeader * ah = dynamic_cast<AliAODHeader*>(aod->GetHeader());
 848         if(!ah) AliFatal("Not a standard AOD");
 849         hdr->SetESDFileName( ah->GetESDFileName() );
 850         hdr->SetEventNumberESDFile( ah->GetEventNumberESDFile() );
 851 }
 852
 853 //______________________________________________________________________________
 854 Int_t AliJFilter::GetSuperModuleNumber(bool isemcal, AliVCluster *cluster, AliVCaloCells *cells, Int_t absId)
 855 {
 856         //get super module number
 857         if(isemcal){
 858                 Int_t absIdMax  = -1, iSM =-1, ieta = -1, iphi = -1;
 859                 Bool_t shared = kFALSE;
 860                 fEMCALRecoUtils->GetMaxEnergyCell(fEMCALGeometry, cells, cluster, absIdMax,  iSM, ieta, iphi, shared);
 861
 862                 if(iSM < 0 || iphi < 0 || ieta < 0 )
 863                 {
 864                         AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",
 865                                                 iSM,ieta,iphi));
 866                 }
 867
 868                 return iSM ;
 869
 870         } else {
 871                 Int_t    relId[4];
 872                 if ( absId >= 0) {
 873                         fPHOSGeom->AbsToRelNumbering(absId,relId);
 874                         fPHOSGeom->AbsToRelNumbering(absId,relId);
 875                         return relId[0]-1;
 876                 } else return -1;
 877         }//PHOS
 878
 879         return -1;
 880 }
 881
 882 //______________________________________________________________________________
 883 Double_t * AliJFilter::GetCellsAmplitude( bool isemcal, AliVCluster *cluster, AliVCaloCells *emCells, AliVCaloCells *phoCells )
 884 {
 885         // cell amplitude reader
 886         Int_t iCell, nCell;
 887         UShort_t *cellAddrs;
 888         Double_t *amps;
 889
 890         // get cluster cells
 891         nCell = cluster->GetNCells();
 892
 893         amps = new Double_t[nCell];
 894
 895         // get the cell addresses
 896         cellAddrs = cluster->GetCellsAbsId();
 897
 898         // get the cell amplitudes
 899         for( iCell = 0; iCell < nCell; iCell++ ){
 900                 if( isemcal )
 901                         amps[iCell] = emCells->GetCellAmplitude( cellAddrs[iCell] );
 902                 else
 903                         amps[iCell] = phoCells->GetCellAmplitude( cellAddrs[iCell] );
 904
 905         }
 906
 907         return amps;
 908 }
 909
 910 //_____________________________________________________________________________
 911
 912 UInt_t AliJFilter::ConvertTriggerMask(){
 913
 914         //convert alice trigger mask to jcorran trigger mask
 915         UInt_t triggerMaskJC=0;
 916         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 917                         ->IsEventSelected() & AliVEvent::kMB){
 918                 // minimum bias TBit 0
 919                 triggerMaskJC |= (1<<kMinBiasTriggerBitJCorran);
 920                 fAliJRunHeader->SetActiveTriggersJCorran( kMinBiasTriggerBitJCorran, "MinBiasTriggerBitJCorran");
 921         }
 922
 923         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 924                         ->IsEventSelected() & AliVEvent::kHighMult){
 925                 //high multiplicity trigger TBit 1
 926                 triggerMaskJC |= (1<<kHighMultTriggerBitJCorran);
 927                 fAliJRunHeader->SetActiveTriggersJCorran( kHighMultTriggerBitJCorran,"HighMultTriggerBitJCorran");
 928         }
 929
 930         if((((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 931                                 ->IsEventSelected() & AliVEvent::kEMC1) ||
 932                         (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 933                          ->IsEventSelected() & AliVEvent::kEMC7 )){
 934                 //EMCAL L0   TBit2
 935                 triggerMaskJC |= (1<<kEmc0TriggerBitJCorran);
 936                 fAliJRunHeader->SetActiveTriggersJCorran( kEmc0TriggerBitJCorran,"Emc0TriggerBitJCorran");
 937         }
 938
 939         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 940                         ->IsEventSelected() & AliVEvent::kEMCEGA){
 941                 //EMCAL Gamma TBit3
 942                 triggerMaskJC |= (1<<kEmc1GammaTriggerBitJCorran);
 943                 fAliJRunHeader->SetActiveTriggersJCorran( kEmc1GammaTriggerBitJCorran,"Emc1GammaTriggerBitJCorran");
 944         }
 945
 946         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 947                         ->IsEventSelected() & AliVEvent::kEMCEJE){
 948                 //EMCAL JET TBit4
 949                 triggerMaskJC |= (1<<kEmc1JetTriggerBitJCorran);
 950                 fAliJRunHeader->SetActiveTriggersJCorran( kEmc1JetTriggerBitJCorran,"Emc1JetTriggerBitJCorran");
 951         }
 952
 953         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 954                         ->IsEventSelected() & AliVEvent::kCentral){
 955                 //central trigger TBit 5
 956                 triggerMaskJC |= (1<<kCentralTriggerBitJCorran);
 957                 fAliJRunHeader->SetActiveTriggersJCorran( kCentralTriggerBitJCorran,"CentralTriggerBitJCorran");
 958         }
 959
 960         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 961                         ->IsEventSelected() & AliVEvent::kSemiCentral){
 962                 //semi-central trigger TBit 6
 963                 triggerMaskJC |= (1<<kSemiCentralTriggerBitJCorran);
 964                 fAliJRunHeader->SetActiveTriggersJCorran( kSemiCentralTriggerBitJCorran,"SemiCentralTriggerBitJCorran");
 965         }
 966
 967         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 968                         ->IsEventSelected() & AliVEvent::kFastOnly){
 969                 //semi-central trigger TBit 6
 970                 triggerMaskJC |= (1<<kFastOnlyBitJCorran);
 971                 fAliJRunHeader->SetActiveTriggersJCorran( kFastOnlyBitJCorran ,"FastOnlyBitJCorran");
 972         }
 973
 974         if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 975                         ->IsEventSelected() & AliVEvent::kINT7){
 976                 // minimum bias TBit 0
 977                 triggerMaskJC |= (1<<kINT7TriggerBitJCorran);
 978                 fAliJRunHeader->SetActiveTriggersJCorran( kINT7TriggerBitJCorran, "INT7TriggerBitJCorran");
 979         }
 980
 981         return triggerMaskJC;
 982 }
 983
 984
 985 //______________________________________________________________________________
 986 void AliJFilter::ReadMCTracksFromESD(){
 987         //store MC information from AliStack
 988         if(!MCEvent()) return;
 989         AliStack *stack = MCEvent()->Stack();
 990         if(!stack) return;
 991         Int_t np    = MCEvent()->GetNumberOfTracks();
 992
 993         //  AliGenEventHeader* genHeader = fMC->GenEventHeader();
 994         //  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 995         //  Double_t ptHard = 0;
 996         //  Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 997         //  nTrials = pythiaGenHeader->Trials();
 998         //  ptHard  = pythiaGenHeader->GetPtHard();
 999         //  Int_t nprim = stack->GetNtrack();
1000
1001         Long64_t ntrack = 0;
1002
1003         for(Long64_t iTrack = 0; iTrack < np; iTrack++){
1004                 AliMCParticle *track = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
1005                 if(!track){
1006                         Printf("ERROR: Could not receive track %d",(int) iTrack);
1007                         continue;
1008                 }
1009                 Bool_t isPrimary = stack->IsPhysicalPrimary(iTrack);
1010                 if(isPrimary){
1011                         //create a new JMCTrack and fill the track info
1012                         AliJMCTrack *ctrack = new( (*fMCTrackList)[ntrack++] ) AliJMCTrack;
1013
1014                         TParticle *partStack = stack->Particle(iTrack);
1015                         Int_t   pdg  = partStack->GetPdgCode();
1016
1017                         Char_t ch     = (Char_t) partStack->GetPDG()->Charge();
1018                         Int_t label    = track->GetLabel();
1019
1020                         ctrack->SetLabel(label);
1021                         ctrack->SetPdgCode(pdg);
1022                         ctrack->SetPxPyPzE( partStack->Px(), partStack->Py(), partStack->Pz(), partStack->Energy());
1023                         ctrack->SetCharge(ch);
1024                         ctrack->SetFlag(AliJMCTrack::kPrimary, isPrimary);
1025
1026                         ctrack->SetProductionVertex(partStack->Vx(),partStack->Vy(),partStack->Vz());
1027                 }// loop for al primary tracks
1028         }
1029 }
1030
1031 //--------------------------------------------------------------------
1032 void AliJFilter::ReadMCTracksFromAOD(){
1033         //retreive MC particles from event //FKEFF//
1034         if(!AODEvent()) return;  TClonesArray *mcArray = (TClonesArray*) AODEvent()->
1035                 FindListObject(AliAODMCParticle::StdBranchName());
1036         if(!mcArray){
1037                 Printf("No MC particle branch found");
1038                 return;
1039         }
1040
1041         Long64_t ntrack = 0;
1042         Long64_t np = mcArray->GetEntriesFast();
1043
1044         for(Long64_t it = 0; it < np; it++) {
1045                 AliAODMCParticle *track = (AliAODMCParticle*) mcArray->At(it);
1046                 if(!track){
1047                         Error("ReadEventAODMC", "Could not receive particle %d",(int) it);
1048                         continue;
1049                 }
1050                 bool isPrimary = track->IsPhysicalPrimary();
1051                 if(isPrimary){
1052                         //create a new JMCTrack and fill the track info
1053                         AliJMCTrack *ctrack = new ((*fMCTrackList)[ntrack++]) AliJMCTrack;;
1054
1055                         Int_t   pdg  = track->GetPdgCode();
1056
1057                         Char_t ch     = (Char_t) track->Charge();
1058                         Int_t label    = track->GetLabel();
1059
1060                         ctrack->SetLabel(label);
1061                         ctrack->SetPdgCode(pdg);
1062                         ctrack->SetPxPyPzE( track->Px(), track->Py(), track->Pz(), track->E());
1063                         ctrack->SetCharge(ch);
1064                         ctrack->SetFlag(AliJMCTrack::kPrimary, isPrimary);
1065
1066                         ctrack->SetProductionVertex(track->Xv(),track->Yv(),track->Zv());
1067                 }
1068         }
1069
1070 }
1071
1072
1073 //--------------------------------------------------------------------
1074 void AliJFilter::RemapMCLabels(){
1075         // remaps all MC labels to the new arrays
1076
1077         Int_t i, j, label, mother0, mother1;
1078         AliJTrack *track;
1079         AliJPhoton *cluster;
1080         // BS AliJCaloCell *cell;
1081         AliJMCTrack *mctrack;
1082
1083         // tracks
1084         for( i = 0; i < fTrackList->GetEntries(); i++ ){
1085                 track = (AliJTrack*)fTrackList->At( i );
1086
1087                 track->SetLabel( fMcMap->At( track->GetLabel() ));
1088         }
1089
1090         // clusters
1091         if( fAliJRunHeader->GetStoreEMCalInfo() ){
1092                 for( i = 0; i < fPhotonList->GetEntries(); i++ ){
1093                         cluster = (AliJPhoton*)fPhotonList->At( i );
1094                         for( j = 0; j < cluster->GetNEMCLabel(); j++ ){
1095                                 label = cluster->GetEMCLabel( j );
1096                                 // no label clusters protection
1097                                 if( label >= 0 )
1098                                         cluster->SetEMCLabel( j, fMcMap->At( label ));
1099                         }
1100                 }
1101
1102                 /*  BS
1103                 // cells
1104                 for( i = 0; i < fCaloCellList->GetEntries(); i++ ){
1105                 cell = (AliJCaloCell*)fCaloCellList->At( i );
1106                 label = cell->GetMcLabel();
1107 // no label cells protection
1108 if( label >= 0 )
1109 cell->SetMcLabel( fMcMap->At( cell->GetMcLabel() ));
1110 }
1111 */
1112 }
1113
1114 // MC particles
1115 for( i = 0; i < fMCTrackList->GetEntries(); i++ ){
1116         mctrack = (AliJMCTrack*)fMCTrackList->At( i );
1117
1118         mother0 = mctrack->GetMother( 0 );
1119         mother1 = mctrack->GetMother( 1 );
1120
1121         if( mother0 >= 0 )
1122                 mother0 = fMcMap->At( mother0 );
1123         if( mother1 >= 0 )
1124                 mother1 = fMcMap->At( mother1 );
1125
1126         mctrack->SetMother( mother0, mother1 );
1127 }
1128 }
1129
1130 //--------------------------------------------------------------------
1131
1132
1133 void AliJFilter::PrintOut() const {
1134         //AliJRunHeader * RunInfo = fAliJRunHeader;
1135 }
1136
1137 //********************************************
1138 //    UTILS
1139 //********************************************
1140 void AliJFilter::AddList(const char* aname, const char* cname, TClonesArray **obj, int nlist){
1141         *obj = new TClonesArray(cname, nlist);
1142         (*obj)->SetName(aname);
1143         (*obj)->SetOwner();
1144 }
1145