PWG4/JCORRAN/AliJCORRANTask.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 //______________________________________________________________________________
  17 // Analysis task for high pt particle correlations
  18 // author: R.Diaz, J. Rak,  D.J. Kim
  19 // ALICE Group University of Jyvaskyla
  20 // Finland
  21 // Fill the analysis containers for ESD or AOD
  22 // Adapted for AliAnalysisTaskSE and AOD objects
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25
  26 #include <Riostream.h>
  27 #include <TChain.h>
  28 #include <TVector3.h>
  29 #include <TFile.h>
  30 #include <TH1.h>
  31 #include <TClonesArray.h>
  32 #include <TObjArray.h>
  33 #include <TObjString.h>
  34 #include <TFormula.h>
  35 #include <TString.h>
  36 #include <TRefArray.h>
  37 #include <TNtuple.h>
  38
  39 #include "AliAnalysisTaskSE.h"
  40 #include "AliAODHandler.h"
  41
  42 #include "AliJCORRANTask.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliESDEvent.h"
  45 #include "AliMCEvent.h"
  46 #include "AliStack.h"
  47 #include "AliGenEventHeader.h"
  48 #include "AliGenCocktailEventHeader.h"
  49 #include "AliGenPythiaEventHeader.h"
  50 #include "AliInputEventHandler.h"
  51 #include "AliESDCaloCluster.h"
  52 #include "AliAODEvent.h"
  53 #include "AliAODHeader.h"
  54 #include "AliLog.h"
  55 #include "AliESDVertex.h"
  56 #include "AliESDtrack.h"
  57 #include "AliAODTrack.h"
  58 #include "AliAnalysisFilter.h"
  59 #include "AliPHOSGeoUtils.h"
  60 #include "AliEMCALGeometry.h"
  61 #include "AliESDtrackCuts.h"
  62 #include "AliAODVertex.h"
  63 #include "AliAODTracklets.h"
  64 #include "AliAODPid.h"
  65
  66 #include "AliPhJTrackList.h"
  67 #include "AliPhJMCTrackList.h"
  68 #include "AliPhJPhotonList.h"
  69 #include "AliPhJHeaderList.h"
  70 #include "AliJTrack.h"
  71 #include "AliJPhoton.h"
  72 #include "AliJHeader.h"
  73 #include "AliAODTracklets.h"
  74 #include "AliMultiplicity.h"
  75 #include "JConst.h"
  76 #include "AliESDRun.h"
  77 #include "AliJRunHeader.h"
  78
  79
  80
  81
  82
  83 //______________________________________________________________________________
  84 AliJCORRANTask::AliJCORRANTask() :
  85   AliAnalysisTaskSE("PWG4JCORRAN"),
  86   fInputFormat("ESD"),
  87   fEsdTrackCuts(0x0),
  88   fDownscaling(1),
  89   fLowerCutOnLPMom(0),
  90   fLowerCutOnLeadingCaloClusterE(0),
  91   fLowerCutOnCaloClusterE(0.2),
  92   fIsRealOrMC(0),
  93   fAODName("jcorran.root"),
  94   f1CutMaxDCAToVertexXYPtDep(0x0),
  95   fTrackList(0x0),
  96   fMCTrackList(0x0),
  97   fPhotonList(0x0),
  98   fHeaderList(0x0),
  99   fAliRunHeader(0x0),
 100   fPHOSGeom(0x0)
 101 {
 102   //Default constructor
 103   for(Int_t i=0;i<kRangeTriggerTableAlice;i++)   fActiveTriggers[i]=" ";
 104   for(Int_t i=0;i<kRangeTriggerTableJCorran;i++) fTriggerTableJCorran[i]=" ";
 105
 106   DefineInput (0, TChain::Class());
 107 }
 108
 109 //______________________________________________________________________________
 110 AliJCORRANTask::AliJCORRANTask(const char *name, TString inputformat):
 111   AliAnalysisTaskSE(name),
 112   fInputFormat(inputformat),
 113   fEsdTrackCuts(0x0),
 114   fDownscaling(1),
 115   fLowerCutOnLPMom(0),
 116   fLowerCutOnLeadingCaloClusterE(0),
 117   fLowerCutOnCaloClusterE(0.2),
 118   fIsRealOrMC(0),
 119   fAODName("jcorran.root"),
 120   f1CutMaxDCAToVertexXYPtDep(0x0),
 121   fTrackList(0x0),
 122   fMCTrackList(0x0),
 123   fPhotonList(0x0),
 124   fHeaderList(0x0),
 125   fAliRunHeader(0x0),
 126   fPHOSGeom(0x0)
 127 {
 128   // Constructor
 129   if(fDebug > 5) cout << "---- AliJCORRANTask Constructor ----"<<endl;
 130
 131   for(Int_t i=0;i<kRangeTriggerTableAlice;i++)   fActiveTriggers[i]=" ";
 132   for(Int_t i=0;i<kRangeTriggerTableJCorran;i++) fTriggerTableJCorran[i]=" ";
 133
 134   DefineInput (0, TChain::Class());
 135 }
 136
 137 //____________________________________________________________________________
 138 AliJCORRANTask::AliJCORRANTask(const AliJCORRANTask& ap) :
 139   AliAnalysisTaskSE(ap.GetName()),
 140   fInputFormat(ap.fInputFormat),
 141   fEsdTrackCuts(ap.fEsdTrackCuts),
 142   fDownscaling(ap.fDownscaling),
 143   fLowerCutOnLPMom(ap.fLowerCutOnLPMom),
 144   fLowerCutOnLeadingCaloClusterE(ap.fLowerCutOnLeadingCaloClusterE),
 145   fLowerCutOnCaloClusterE(ap.fLowerCutOnCaloClusterE),
 146   fIsRealOrMC(ap.fIsRealOrMC),
 147   fAODName(ap.fAODName),
 148   f1CutMaxDCAToVertexXYPtDep(ap.f1CutMaxDCAToVertexXYPtDep),
 149   fTrackList(ap.fTrackList),
 150   fMCTrackList(ap.fMCTrackList),
 151   fPhotonList(ap.fPhotonList),
 152   fHeaderList(ap.fHeaderList),
 153   fAliRunHeader(ap.fAliRunHeader),
 154   fPHOSGeom(ap.fPHOSGeom)
 155 {
 156   // cpy ctor
 157   for(int k=0; k < kRangeTriggerTableAlice; k++)
 158     fActiveTriggers[k] = ap.fActiveTriggers[k];
 159
 160   for(int j=0; j < kRangeTriggerTableJCorran; j++)
 161     fTriggerTableJCorran[j] = ap.fTriggerTableJCorran[j];
 162 }
 163
 164 //_____________________________________________________________________________
 165 AliJCORRANTask& AliJCORRANTask::operator = (const AliJCORRANTask& ap)
 166 {
 167 // assignment operator
 168
 169   this->~AliJCORRANTask();
 170   new(this) AliJCORRANTask(ap);
 171   return *this;
 172 }
 173
 174 //______________________________________________________________________________
 175 AliJCORRANTask::~AliJCORRANTask()
 176 {
 177   // destructor
 178   if(f1CutMaxDCAToVertexXYPtDep) delete f1CutMaxDCAToVertexXYPtDep;
 179   if(fTrackList)    delete fTrackList;
 180   if(fMCTrackList)  delete fMCTrackList;
 181   if(fPhotonList)   delete fPhotonList;
 182   if(fHeaderList)   delete fHeaderList;
 183   if(fAliRunHeader) delete fAliRunHeader;
 184   if(fPHOSGeom)     delete fPHOSGeom  ;
 185   GetEMCALGeoUtils(kTRUE);
 186 }
 187
 188 //________________________________________________________________________
 189 void AliJCORRANTask::UserCreateOutputObjects()
 190 {
 191   // create the jcorran outputs objects
 192   fTrackList    = new AliPhJTrackList(kALICE);
 193   if(fIsRealOrMC) fMCTrackList  = new AliPhJMCTrackList(kALICE);
 194   fPhotonList   = new AliPhJPhotonList(kALICE);
 195   fHeaderList   = new AliPhJHeaderList(kALICE);
 196
 197   fAliRunHeader = new AliJRunHeader();
 198
 199   fPHOSGeom  = new AliPHOSGeoUtils("PHOSgeo") ;
 200
 201   // create the jcorran output deltaAOD
 202   //if(fDebug > 5) cout << "AliJCORRANTask UserCreateOutputObjects----------------------"<<endl;
 203
 204   if(fDebug > 1) printf("AliJCORRANTask::UserCreateOutPutData() \n");
 205   if(!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) {
 206     Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
 207     return;
 208   }
 209
 210   AddAODBranch("AliPhJTrackList", &fTrackList,fAODName.Data());
 211   AddAODBranch("AliPhJPhotonList", &fPhotonList,fAODName.Data());
 212   AddAODBranch("AliPhJHeaderList", &fHeaderList,fAODName.Data());
 213
 214   if(fIsRealOrMC){
 215     AddAODBranch("AliPhJMCTrackList", &fMCTrackList, fAODName.Data());
 216   }
 217
 218 }
 219
 220 //______________________________________________________________________________
 221 void AliJCORRANTask::UserExec(Option_t* /*option*/)
 222 {
 223   // Processing of one event
 224   //if(fDebug > 5) cout << "------- AliJCORRANTask Exec-------"<<endl;
 225   if(!((Entry()-1)%100))
 226       AliInfo(Form(" Processing event # %lld",  Entry()));
 227
 228
 229   Bool_t storeEvent = kFALSE;//based on offline trigger decide whether to store the event or not
 230   if(fIsRealOrMC){
 231     storeEvent = kTRUE; //store all MC events
 232   }else{ //when we are processing real events store only selected events
 233     if(StoreDownscaledMinBiasEvent() ||
 234       ContainsESDHighPtTrack()   ||
 235       ContainsESDHighECaloClusters()){
 236         storeEvent = kTRUE;
 237     }
 238   }
 239
 240   fTrackList->Reset();
 241   if(fIsRealOrMC) fMCTrackList->Reset();
 242   fPhotonList->Reset();
 243   fHeaderList->Reset();
 244
 245   static int runId=-1;
 246
 247   if(fInputFormat=="ESD"){   //   Reading ESD
 248     AliESDEvent* esd = (AliESDEvent*)InputEvent();
 249     if(!esd) return;
 250                        //cout << storeEvent<<"    "<<esd->GetFiredTriggerClasses() << endl;
 251     AliMCEvent* mcEvent = NULL;
 252     if(fIsRealOrMC)  mcEvent = MCEvent();
 253
 254      //=========== FILL AND STORE RUN HEADER   (ONLY ONCE PER RUN) =============
 255      if(esd->GetRunNumber() != runId){ //new run has started
 256
 257        const AliESDRun* esdRun = esd->GetESDRun();
 258
 259        for(Int_t triggerBit=0; triggerBit<kRangeTriggerTableAlice; triggerBit++){
 260          fActiveTriggers[triggerBit] = esdRun->GetTriggerClass(triggerBit);
 261        }
 262        runId = esd->GetRunNumber();
 263
 264        //Polarity of magnetic field in L3 solenoid
 265        Short_t l3MgFieldPolarity=0; // (LHC convention: +current -> +Bz)
 266        if(esd->GetCurrentL3() >0) l3MgFieldPolarity =  1;
 267        if(esd->GetCurrentL3() <0) l3MgFieldPolarity = -1;
 268
 269        //Create internal JCorran trigger mask.  Mapping between trigger and trigger bit
 270        fTriggerTableJCorran[kMinBiasTriggerBitJCorran]="Minimum Bias";//minimum bias occupies first trigger bit
 271        fTriggerTableJCorran[kHighMultTriggerBitJCorran]="High Multiplicity";//high multiplicity trigger => second trigger bit
 272
 273        //=========== Fill Run header object ===============
 274        fAliRunHeader->SetRunNumber(runId);
 275        fAliRunHeader->SetActiveTriggersAlice(fActiveTriggers);
 276        fAliRunHeader->SetL3Field(l3MgFieldPolarity, esd->GetMagneticField());
 277        fAliRunHeader->SetActiveTriggersJCorran(fTriggerTableJCorran,kRangeTriggerTableJCorran);
 278
 279       //Store Run header
 280       ( OutputTree()->GetUserInfo())->Add(fAliRunHeader);
 281     }//end RunHeader
 282
 283     if(storeEvent){ //FK//
 284       //-------------- reset all the arrays -------------
 285       //store event only when it is downscaled min bias
 286       // or contais high pt hadron
 287       // or contains high energy cluster in EMCAL or PHOS
 288       ReadESDTracks(esd);
 289       ReadESDCaloClusters(esd);
 290       ReadESDHeader(esd);
 291       if(fIsRealOrMC) ReadMCTracks(mcEvent);
 292     }
 293   }else if( fInputFormat == "AOD") {
 294
 295     AliAODEvent* aod = AODEvent();
 296     if(!aod) return;
 297
 298     //=========== FILL AND STORE RUN HEADER   (ONLY ONCE PER RUN) =============
 299     if(aod->GetRunNumber() != runId){ //new run has started
 300
 301       runId = aod->GetRunNumber();
 302       // trigger names???
 303
 304       //Polarity of magnetic field in L3 solenoid
 305       Short_t l3MgFieldPolarity=0; // (LHC convention: +current -> +Bz)
 306       if( aod->GetMagneticField() > 0 ) l3MgFieldPolarity =  1;
 307       if( aod->GetMagneticField() < 0 ) l3MgFieldPolarity = -1;
 308
 309       fAliRunHeader->SetRunNumber(runId);
 310       fAliRunHeader->SetL3Field(l3MgFieldPolarity, aod->GetMagneticField());
 311       //Store Run header
 312       (OutputTree()->GetUserInfo())->Add(fAliRunHeader);
 313
 314     }
 315
 316     ReadAODTracks(aod);
 317     ReadAODCaloClusters(aod);
 318     ReadAODHeader(aod);
 319
 320   }else{
 321     cout << "Error: Not correct InputDataFormat especified " << endl;
 322     return;
 323   }
 324
 325
 326   if(fTrackList->GetNTracks() > 0 ||
 327     fPhotonList->GetNPhotons() >0 ||
 328     ( fIsRealOrMC && fMCTrackList->GetNTracks()>0)){
 329
 330     AliAODHandler* outputHandler =
 331     (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 332     outputHandler->SetFillAOD(kTRUE);
 333   }
 334 }
 335
 336 //______________________________________________________________________________
 337 void AliJCORRANTask::Init()
 338 {
 339   // Intialisation of parameters
 340   AliInfo("Doing initialization") ;
 341
 342   TString formula(fEsdTrackCuts->GetMaxDCAToVertexXYPtDep());
 343   if(formula.Length()>0){ // momentum dep DCA cut for AOD
 344     formula.ReplaceAll("pt","x");
 345     if(f1CutMaxDCAToVertexXYPtDep) delete f1CutMaxDCAToVertexXYPtDep;
 346     f1CutMaxDCAToVertexXYPtDep = new TFormula("f1CutMaxDCAToVertexXYPtDep",formula.Data());
 347   }
 348 }
 349
 350 //______________________________________________________________________________
 351 void AliJCORRANTask::Terminate(Option_t *)
 352 {
 353   // Processing when the event loop is ended
 354  // OutputTree()->Print();
 355
 356   // ((AliJRunHeader *) (( OutputTree()->GetUserInfo())->First()))->PrintOut();
 357
 358   cout<<"PWG4JCORRAN Analysis DONE !!"<<endl;
 359 }
 360
 361 //______________________________________________________________________________
 362 void AliJCORRANTask::ReadESDTracks(const AliESDEvent * esd)
 363 {
 364     // Read the AliESDtrack and fill the list of AliJTrack containers
 365     Int_t nt = esd->GetNumberOfTracks();
 366    // if(fDebug < 5) cout << "ESD::NumberOfTracks = " << nt << endl;
 367     Short_t ntrk = 0;
 368     Double_t pid[10];
 369
 370     //loop over tracks
 371     for(Int_t it = 0; it < nt; it++) {
 372
 373         AliESDtrack *track = esd->GetTrack(it);
 374         if(! fEsdTrackCuts->IsSelected(track)) continue; // apply track selection criteria
 375
 376         //------------ T P C ------------
 377         Float_t tpcDca[2], tpcCov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 378         track->GetImpactParametersTPC(tpcDca,tpcCov);
 379
 380         Int_t nClust         = track->GetTPCclusters(0);
 381         Int_t nFindableClust = track->GetTPCNclsF();
 382         Float_t tpcChi2PerCluster = 0.;
 383         if(nClust>0.) tpcChi2PerCluster = track->GetTPCchi2()/Float_t(nClust);
 384
 385         Float_t tpcClustPerFindClust = 0.;
 386         if(nFindableClust>0.) tpcClustPerFindClust = Float_t(nClust)/nFindableClust;
 387         //--------------------------------
 388
 389              //create a new AliJTrack and fill the track info
 390         fTrackList->AddAliJTrack(ntrk);
 391         AliJTrack *ctrack = fTrackList->GetAliJTrack(ntrk);
 392
 393         ctrack->SetPtot(track->P());//here
 394         ctrack->SetPt(track->Pt());
 395         ctrack->SetTheta(track->Theta());
 396         ctrack->SetPhi(atan2(sin(track->Phi()), cos(track->Phi())));
 397         track->GetESDpid(pid);
 398         ctrack->SetPID(pid);
 399
 400         ctrack->SetFlavor(kNone);
 401         ctrack->SetCharge(track->Charge());
 402         ctrack->ConvertAliPID();
 403         ctrack->SetEta(track->Eta());
 404
 405         Int_t itof = track->GetTOFcluster(); // index of the assigned TOF cluster
 406         if(itof>=0) ctrack->SetTof(track->GetTOFsignal());
 407
 408         ctrack->SetTPCdEdx(track->GetTPCsignal());
 409         ctrack->SetTPCnClust(track->GetTPCNcls());
 410         ctrack->SetTPCDCAXY(tpcDca[0]);
 411         ctrack->SetTPCDCAZ(tpcDca[1]);
 412         ctrack->SetTPCClustPerFindClust(tpcClustPerFindClust);
 413         ctrack->SetTPCChi2PerClust(tpcChi2PerCluster);
 414
 415         Float_t impactXY, impactZ; //get track impact parameters
 416         track->GetImpactParameters(impactXY,impactZ);
 417         ctrack->SetImapactXY(impactXY);
 418         ctrack->SetImapactZ(impactZ);
 419
 420         ctrack->SetKinkIndex(track->GetKinkIndex(0));
 421         ctrack->SetStatus(track->GetStatus()); // ULong_t
 422
 423         if(fIsRealOrMC){
 424           ctrack->SetITSLabel(track->GetITSLabel());
 425           ctrack->SetTPCLabel(track->GetTPCLabel());
 426         }
 427
 428         fTrackList->SetNTracks(++ntrk);
 429
 430      } // end tracks loop
 431 }
 432
 433 //______________________________________________________________________________
 434 void AliJCORRANTask::ReadAODTracks(const AliAODEvent * aod)
 435 {
 436     // Read the AliAODtrack and fill the list of AliJTrack containers
 437     Int_t nt = aod->GetNumberOfTracks();
 438     if(fDebug < 5) cout << "AOD::NumberOfTracks = " << nt << endl;
 439     Short_t ntrk = 0;
 440     //loop over tracks
 441     for(Int_t it = 0; it < nt; it++) {
 442
 443         AliAODTrack *track = aod->GetTrack(it);
 444         if(!AcceptAODTrack(track)) continue;
 445         //if(! fEsdTrackCuts->IsSelected(track)) continue; //apply loose selection criteria
 446        //FK//if(track->GetType() != AliAODTrack::kPrimary) continue; // only primaries
 447
 448         //create a new AliJTrack and fill the track info
 449         fTrackList->AddAliJTrack(ntrk);
 450         AliJTrack *ctrack = fTrackList->GetAliJTrack(ntrk);
 451
 452         ctrack->SetPtot(track->P());
 453         ctrack->SetPt(track->Pt());
 454         ctrack->SetTheta(track->Theta());
 455         //ctrack->SetPhi(track->Phi());
 456         ctrack->SetPhi(atan2(sin(track->Phi()), cos(track->Phi())));
 457         ctrack->SetPID((Double_t*)track->PID());
 458         ctrack->SetFlavor(kNone);
 459         ctrack->SetCharge(track->Charge());
 460         ctrack->SetEta(track->Eta());
 461
 462         AliAODPid* pidAOD = track->GetDetPid();
 463         if(pidAOD){
 464           ctrack->SetTof(pidAOD->GetTOFsignal());
 465           ctrack->SetTPCdEdx(pidAOD->GetTPCsignal());
 466         }
 467
 468         Double_t impactDCA[3];
 469         if( track->GetPosition(impactDCA)){
 470           ctrack->SetImapactXY(sqrt(impactDCA[0]*impactDCA[0] + impactDCA[1]*impactDCA[1]));//impactXY);
 471           ctrack->SetImapactZ(impactDCA[2]);//impactZ);
 472         }
 473
 474         ctrack->SetChi2perNDF(track->Chi2perNDF());
 475         ctrack->SetTPCnClust(track->GetTPCNcls());
 476         ctrack->SetChi2Trig(track->GetChi2MatchTrigger());
 477
 478         ctrack->SetStatus(track->GetStatus());//
 479         //ctrack->SetRecFlags(track->GetFlags());//?status
 480
 481         if(fIsRealOrMC){
 482           Int_t  label = track->GetLabel();
 483           ctrack->SetITSLabel(label);
 484           ctrack->SetTPCLabel(label);
 485         }
 486
 487         fTrackList->SetNTracks(++ntrk);
 488
 489      } // end tracks loop
 490 }
 491
 492 //______________________________________________________________________________
 493 void AliJCORRANTask::ReadESDCaloClusters(const AliESDEvent* esd)
 494 {
 495   // Read the AliESDCaloClusters and fill the list of AliJPhoton containers
 496   Short_t nPhotons = 0;
 497   Int_t numberOfCaloClusters = esd->GetNumberOfCaloClusters() ;
 498   if(fDebug < 5) cout << "ESD::number of ESD caloclusters " << numberOfCaloClusters << endl;
 499
 500   // loop over all the Calo Clusters
 501   for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 502
 503     AliESDCaloCluster *caloCluster = esd->GetCaloCluster(icluster) ;
 504     if(!caloCluster) continue;
 505     if(caloCluster->GetNTracksMatched()<=0){
 506       if(caloCluster->E()<fLowerCutOnCaloClusterE) continue;
 507
 508       // we will not implement any PID cut here
 509       fPhotonList->AddAliJPhoton(nPhotons);
 510       AliJPhoton *pht = fPhotonList->GetAliJPhoton(nPhotons);
 511       pht->SetFlavor(kNone);//kPhoton);
 512       pht->SetE(caloCluster->E());
 513       pht->SetChi2(caloCluster->Chi2());
 514       pht->SetPID(caloCluster->GetPID());
 515       Float_t pos[3]; caloCluster->GetPosition(pos) ;
 516       pht->SetX(pos[0]);
 517       pht->SetY(pos[1]);
 518       pht->SetZ(pos[2]);
 519       pht->SetPhi(atan2(pos[1],pos[0]));
 520       pht->SetTheta(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2]));
 521       pht->SetPt(caloCluster->E()*sin(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2])));
 522       pht->SetCharge(0);
 523       if(caloCluster->IsEMCAL()) pht->SetCaloType(AliJPhoton::kEMCALCalo);
 524       if(caloCluster->IsPHOS())  pht->SetCaloType(AliJPhoton::kPHOSCalo);
 525       pht->SetDistToBadChannel(caloCluster->GetDistanceToBadChannel());
 526       pht->SetDispersion(caloCluster->GetDispersion());
 527       pht->SetM20(caloCluster->GetM20());
 528       pht->SetM02(caloCluster->GetM02());
 529       pht->SetEmcCpvDist(caloCluster->GetEmcCpvDistance());
 530       pht->SetNCells(caloCluster->GetNCells());
 531       pht->SetCellsAbsId(caloCluster->GetCellsAbsId());
 532       Int_t imoduleID = GetSuperModuleNumber(caloCluster->IsEMCAL(), caloCluster->GetCellAbsId(0));
 533       pht->SetSuperModuleID(imoduleID);
 534
 535       fPhotonList->SetNPhotons(++nPhotons);
 536     } // end if
 537   } //PHOS and EMCAL clusters
 538
 539 }
 540
 541 //______________________________________________________________________________
 542 void AliJCORRANTask::ReadAODCaloClusters(const AliAODEvent* aod)
 543 {
 544   // Read the AliAODCaloClusters and fill the list of AliJPhoton containers
 545   Int_t numberOfCaloClusters = aod->GetNumberOfCaloClusters();
 546   if(fDebug < 5) cout << "AOD::number of ESD caloclusters " << numberOfCaloClusters << endl;
 547   Short_t nPhotons = 0;
 548
 549   for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 550
 551       AliAODCaloCluster * caloCluster = aod->GetCaloCluster(icluster) ;
 552       if(!caloCluster) continue;
 553       if(caloCluster->GetNTracksMatched() > 0) continue;
 554       if(caloCluster->E() < fLowerCutOnCaloClusterE) continue;
 555       // we will not implement any PID cut here
 556       fPhotonList->AddAliJPhoton(nPhotons);
 557
 558       AliJPhoton *pht = fPhotonList->GetAliJPhoton(nPhotons);
 559       pht->SetFlavor(kNone);
 560       pht->SetE(caloCluster->E());
 561       pht->SetChi2(caloCluster->Chi2());
 562       pht->SetPID((Double_t*)caloCluster->GetPID());
 563       Float_t pos[3]; caloCluster->GetPosition(pos);
 564       pht->SetX(pos[0]);
 565       pht->SetY(pos[1]);
 566       pht->SetZ(pos[2]);
 567       pht->SetPhi(atan2(pos[1],pos[0]));
 568       pht->SetTheta(atan2(sqrt(pos[0]*pos[1]+pos[1]*pos[1]),pos[2]));
 569       pht->SetPt(caloCluster->E()*sin(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2])));
 570       pht->SetCharge(0);
 571       if(caloCluster->IsEMCAL()) pht->SetCaloType(AliJPhoton::kEMCALCalo);
 572       if(caloCluster->IsPHOS())  pht->SetCaloType(AliJPhoton::kPHOSCalo);
 573       pht->SetDistToBadChannel(caloCluster->GetDistanceToBadChannel());
 574       pht->SetDispersion(caloCluster->GetDispersion());
 575       pht->SetM20(caloCluster->GetM20());
 576       pht->SetM02(caloCluster->GetM02());
 577       pht->SetEmcCpvDist(caloCluster->GetEmcCpvDistance());
 578       pht->SetNCells(int(caloCluster->GetNCells()));
 579       pht->SetCellsAbsId(caloCluster->GetCellsAbsId());
 580       Int_t imoduleID = GetSuperModuleNumber(caloCluster->IsEMCAL(), caloCluster->GetCellAbsId(0));
 581       pht->SetSuperModuleID(imoduleID);
 582
 583       fPhotonList->SetNPhotons(++nPhotons);
 584
 585   } // clusters
 586 }
 587
 588 //______________________________________________________________________________
 589 void AliJCORRANTask::ReadESDHeader(const AliESDEvent *esd)
 590 {
 591     // Read the AliESDEvent and fill the list of AliJHeader containers
 592     Short_t nHeaders = 0;
 593     //create a header and fill it
 594     fHeaderList->AddAliJHeader(nHeaders);
 595     AliJHeader *hdr = fHeaderList->GetAliJHeader(nHeaders);
 596
 597     AliMultiplicity *fSPDMult =(AliMultiplicity *) esd->GetMultiplicity();
 598     if(fSPDMult) hdr->SetSPDTrackletMult(fSPDMult->GetNumberOfTracklets());
 599
 600     hdr->SetEventID( esd->GetEventNumberInFile());
 601
 602     hdr->SetTriggerMaskAlice(esd->GetTriggerMask()); //ULong64_t
 603     hdr->SetTriggerMaskJCorran(ConvertTriggerMask()); //UInt_t
 604
 605     const AliESDVertex * vtxESD = esd->GetVertex();
 606     if(vtxESD){
 607       hdr->SetZVertex(vtxESD->GetZv());
 608       hdr->SetZVertexErr(vtxESD->GetZRes());
 609     }
 610     hdr->SetEventType(esd->GetEventType());
 611
 612     fHeaderList->SetNHeaders(++nHeaders);
 613 }
 614
 615 //______________________________________________________________________________
 616 void AliJCORRANTask::ReadAODHeader(const AliAODEvent *aod)
 617 {
 618    //Read the AliAODEvent and fill the list of AliJHeader containers
 619    static int eventID = 0; //FK//?? dummy indexing of events (I cannot see how to get evet ID from AOD)
 620
 621     //read AOD event header
 622     Short_t nHeaders = 0;
 623
 624     //create a header and fill it
 625     fHeaderList->AddAliJHeader(nHeaders);
 626     AliJHeader *hdr = fHeaderList->GetAliJHeader(nHeaders);
 627
 628     const AliAODTracklets *trackletsSPD = aod->GetTracklets();
 629     if(trackletsSPD){
 630       hdr->SetSPDTrackletMult(trackletsSPD->GetNumberOfTracklets());
 631     }
 632     hdr->SetEventID( eventID++ );  //FK//?? I cannot see how to get evet ID from AOD
 633     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
 634     if(vtxAOD){
 635       hdr->SetZVertex((float) vtxAOD->GetZ());
 636       float sigmaVtx[3];
 637       vtxAOD->GetSigmaXYZ(sigmaVtx);
 638       hdr->SetZVertexErr((float) sqrt(sigmaVtx[2]));//TMath::Sqrt(fCovZZ)
 639     }
 640
 641     //load aod event header
 642     AliAODHeader * aodh = aod->GetHeader();
 643     if(aodh){
 644       hdr->SetCentrality(int(aodh->GetCentrality()));
 645       hdr->SetEventType(aodh->GetEventType());
 646       hdr->SetTriggerMaskAlice(aodh->GetTriggerMask()); //ULong64_t
 647     }
 648
 649     hdr->SetTriggerMaskJCorran(ConvertTriggerMask()); //UInt_t
 650
 651     fHeaderList->SetNHeaders(++nHeaders);
 652 }
 653
 654 //______________________________________________________________________________
 655 Int_t AliJCORRANTask::GetSuperModuleNumber(bool isemcal, Int_t absId)
 656 {
 657   //get super module number
 658   if(isemcal){
 659     return GetEMCALGeoUtils()->GetSuperModuleNumber(absId) ;
 660   } else {
 661     Int_t    relId[4];
 662     if ( absId >= 0) {
 663       fPHOSGeom->AbsToRelNumbering(absId,relId);
 664       return relId[0]-1;
 665     } else return -1;
 666   }//PHOS
 667
 668   return -1;
 669 }
 670
 671 //_____________________________________________________________________________
 672
 673 UInt_t AliJCORRANTask::ConvertTriggerMask(){
 674
 675   //convert alice trigger mask to jcorran trigger mask
 676   UInt_t triggerMaskJC=0;
 677   if(!fIsRealOrMC){  //REAL data
 678     if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 679        ->IsEventSelected() & AliVEvent::kMB){
 680       // minimum bias TBit 0
 681       triggerMaskJC += (1<<kMinBiasTriggerBitJCorran);
 682     }
 683
 684     if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 685        ->IsEventSelected() & AliVEvent::kHighMult){
 686       //high multiplicity trigger TBit 1
 687       triggerMaskJC += (1<<kHighMultTriggerBitJCorran);
 688     }
 689   }else{
 690     triggerMaskJC=1; //MC data, at the moment all events filled as MB
 691   }
 692
 693   return triggerMaskJC;
 694 }
 695
 696
 697 //______________________________________________________________________________
 698 bool AliJCORRANTask::StoreDownscaledMinBiasEvent(){
 699   //Decide whether to downscale this MinBiasEvent and store it
 700   bool isThisEventToBeStored = kFALSE;
 701   static Long_t evtNumber=0;
 702
 703   if( ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() ){
 704     //collision candidate
 705
 706     if(evtNumber % fDownscaling ==0){ isThisEventToBeStored = kTRUE; } //store every Xth collision candidate event
 707     evtNumber++;
 708   }
 709   return isThisEventToBeStored;
 710 }
 711 //______________________________________________________________________________
 712 bool AliJCORRANTask::ContainsESDHighPtTrack(){
 713   //If there was an identified high pT particle above threshold reutrn flag to store this event
 714   bool isThisEventToBeStored = kFALSE;
 715
 716   if(fInputFormat=="ESD"){  // E S D
 717
 718     AliESDEvent* esd = (AliESDEvent*)InputEvent();
 719     if(!esd) return kFALSE;
 720     Int_t nt = esd->GetNumberOfTracks();
 721
 722     for(Int_t it = 0; it < nt; it++){
 723       AliESDtrack *track = esd->GetTrack(it);
 724       //Does event contain high pt particle above thereshold in GeV
 725       if(track->Pt() > fLowerCutOnLPMom && fEsdTrackCuts->IsSelected(track)){
 726         isThisEventToBeStored = kTRUE;
 727         break;
 728       }
 729     }
 730   }else{  // A O D
 731     AliAODEvent* aod=(AliAODEvent*) InputEvent();
 732     if(!aod) return kFALSE;
 733     Int_t nt = aod->GetNumberOfTracks();
 734
 735     for(Int_t it = 0; it < nt; it++) {
 736       AliAODTrack *track = aod->GetTrack(it);
 737       //Does event contain high pt particle above threshold in GeV
 738       if(track->Pt() > fLowerCutOnLPMom && AcceptAODTrack(track)){
 739         isThisEventToBeStored = kTRUE;
 740         break;
 741       }
 742     }
 743   }
 744
 745   return isThisEventToBeStored;
 746 }
 747
 748 //______________________________________________________________________________
 749 bool AliJCORRANTask::ContainsESDHighECaloClusters(){
 750   //Check whether there was in the event high E calo clustre and renturn a flag whether to store this event.
 751   bool isThisEventToBeStored = kFALSE;
 752
 753   if(fInputFormat=="ESD"){
 754
 755     AliESDEvent* esd = (AliESDEvent*)InputEvent();
 756     if(!esd) return kFALSE;
 757     Int_t numberOfCaloClusters = esd->GetNumberOfCaloClusters() ;
 758     // loop over all the Calo Clusters
 759     for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 760       AliESDCaloCluster *caloCluster = esd->GetCaloCluster(icluster) ;
 761       if(!caloCluster) continue;
 762       if(caloCluster->GetNTracksMatched() ==-1){
 763         //sotre calo clusters above 1 GeV
 764         if( caloCluster->E() > fLowerCutOnLeadingCaloClusterE){
 765           isThisEventToBeStored = kTRUE;
 766           break;
 767         }
 768       }
 769     }
 770   }else{
 771
 772     AliAODEvent* aod = (AliAODEvent*)InputEvent();
 773     if(!aod) return kFALSE;
 774     Int_t numberOfCaloClusters = aod->GetNumberOfCaloClusters();
 775     // loop over all the Calo Clusters
 776     for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 777       AliAODCaloCluster *caloCluster = aod->GetCaloCluster(icluster) ;
 778       if(!caloCluster) continue;
 779       if(caloCluster->GetNTracksMatched() == -1){
 780         //sotre calo clusters above 1 GeV
 781         if( caloCluster->E() > fLowerCutOnLeadingCaloClusterE){
 782           isThisEventToBeStored = kTRUE;
 783           break;
 784         }
 785       }
 786     }
 787   }
 788
 789   return isThisEventToBeStored;
 790 }
 791
 792
 793 //--------------------------------------------------------------------
 794
 795 void AliJCORRANTask::ReadMCTracks(AliMCEvent *fMC){
 796   //store MC information from AliStack
 797   AliStack *stack = fMC->Stack();
 798   Int_t np    = fMC->GetNumberOfTracks();
 799
 800   //  AliGenEventHeader* genHeader = fMC->GenEventHeader();
 801   //  AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 802   //  Double_t ptHard = 0;
 803   //  Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 804   //  nTrials = pythiaGenHeader->Trials();
 805   //  ptHard  = pythiaGenHeader->GetPtHard();
 806   //  Int_t nprim = stack->GetNtrack();
 807
 808   Short_t ntrack = 0;
 809
 810   for(Int_t iTrack = 0; iTrack < np; iTrack++){
 811     AliMCParticle *track = (AliMCParticle*) fMC->GetTrack(iTrack);
 812     if(!track){
 813       Printf("ERROR: Could not receive track %d", iTrack);
 814       continue;
 815     }
 816     Bool_t isPrimary = fMC->Stack()->IsPhysicalPrimary(iTrack);
 817     if(isPrimary){
 818       //create a new JMCTrack and fill the track info
 819       fMCTrackList->AddJMCTrack(ntrack);
 820       AliJMCTrack *ctrack = fMCTrackList->GetTrack(ntrack);
 821
 822       TParticle *partStack = stack->Particle(iTrack);
 823       Int_t   pdg  = partStack->GetPdgCode();
 824       Float_t engy = partStack->Energy();
 825       Float_t pt   = partStack->Pt();
 826       Float_t ptot = partStack->P();
 827       Float_t eta  = partStack->Eta();
 828       Float_t theta  = partStack->Theta();
 829       Float_t phi    = atan2(sin( partStack->Phi()), cos(partStack->Phi()));
 830       Short_t ch     = (Short_t) partStack->GetPDG()->Charge();
 831       Int_t label    = track->GetLabel();
 832
 833       ctrack->SetLabel(label);
 834       ctrack->SetPdgCode(pdg);
 835       ctrack->SetPt(pt);
 836       ctrack->SetTheta(theta);
 837       ctrack->SetEta(eta);
 838       ctrack->SetPhi(phi);
 839       ctrack->SetE(engy);
 840       ctrack->SetCharge(ch);
 841       ctrack->SetPtot(ptot);
 842       ctrack->SetIsPrimary(isPrimary);
 843
 844       ctrack->SetProductionVertex(partStack->Vx(),partStack->Vy(),partStack->Vz());
 845    /*
 846       Int_t   status  = partStack->GetStatusCode();
 847       ctrack->SetStatusCode(status);
 848
 849       //ctrack->SetPtHard(ptHard);
 850
 851       //bool isInc = (status ==  1 && icode ==  22); //Inclusive
 852       bool ispi0 = (status == 11 && pdg == 111); //kPizero
 853       bool isDgamma = (status == 6 || status == 7) && pdg == 22; // Direct photon
 854       bool inPHOS  = (ispi0||isDgamma)&&fabs(eta)<0.12;
 855       bool inEMCAL = (ispi0||isDgamma)&&fabs(eta)<0.7;
 856       bool inTPC   = fabs(eta)<0.9;
 857       ctrack->SetMother(0,partStack->GetFirstMother());
 858       ctrack->SetMother(1,partStack->GetSecondMother());
 859       ctrack->SetDaughter(0,partStack->GetFirstDaughter());
 860       ctrack->SetDaughter(1,partStack->GetLastDaughter());
 861       ctrack->SetIsInPHOS(inPHOS);
 862       ctrack->SetIsInEMCAL(inEMCAL);
 863       ctrack->SetIsInTPC(inTPC);
 864 */
 865       fMCTrackList->SetNTracks(++ntrack);
 866     }// loop for al primary tracks
 867   }
 868 }
 869
 870 //--------------------------------------------------------------------
 871 bool AliJCORRANTask::AcceptAODTrack(AliAODTrack* aodTrack){
 872   //This function mimics for the AliAODTracks object the AliESDtrackCut function IsSelected
 873   //Cuts are taken from fEsdTrackCuts object
 874   if(fEsdTrackCuts->GetMinNClusterTPC() > aodTrack->GetTPCNcls()) return kFALSE;
 875
 876   //if(fEsdTrackCuts->GetMaxChi2PerClusterTPC() <    );//<-------- how to check?
 877   //      ctrack->SetChi2perNDF(track->Chi2perNDF());
 878
 879   //         C h e c k    r e f i t
 880
 881   if(fEsdTrackCuts->GetRequireTPCRefit()  &&
 882     ((aodTrack->GetStatus() & AliJTrack::kTPCrefit) == 0)) return kFALSE;
 883   if(fEsdTrackCuts->GetRequireITSRefit()  &&
 884     ((aodTrack->GetStatus() & AliJTrack::kITSrefit) == 0)) return kFALSE;
 885
 886
 887    //            C u t s   o n    D C A
 888    Float_t impactDCA[3];
 889    if( aodTrack->GetPosition(impactDCA)){
 890      if((fEsdTrackCuts->GetMaxDCAToVertexXY()>0) &&
 891        (fEsdTrackCuts->GetMaxDCAToVertexXY() < sqrt(impactDCA[0]*impactDCA[0] + impactDCA[1]*impactDCA[1]))) return kFALSE;
 892      if((fEsdTrackCuts->GetMaxDCAToVertexZ()>0) &&
 893        (fEsdTrackCuts->GetMaxDCAToVertexZ() < TMath::Abs(impactDCA[2]))) return kFALSE;
 894    } else return kFALSE;
 895
 896    if(f1CutMaxDCAToVertexXYPtDep){
 897      if(f1CutMaxDCAToVertexXYPtDep->Eval(aodTrack->Pt()) < sqrt(impactDCA[0]*impactDCA[0] + impactDCA[1]*impactDCA[1])) return kFALSE;
 898    }
 899
 900
 901   // if(fEsdTrackCuts->GetAcceptKinkDaughters()) //<--------how to check ?
 902   //  esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
 903
 904
 905   return kTRUE;
 906 }
 907
 908 //--------------------------------------------------------------------
 909 AliEMCALGeometry * AliJCORRANTask::GetEMCALGeoUtils (bool doDelete){
 910   //include EMCAL singleton modul
 911   static AliEMCALGeometry* emcalgeo = 0x0;
 912   if( ! emcalgeo ){
 913     emcalgeo = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
 914   }
 915   if( emcalgeo && doDelete ){ //FK// !emcalgeo
 916     delete emcalgeo;
 917     emcalgeo=0x0;
 918   }
 919   return emcalgeo;
 920 }
 921
 922
 923