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
  41 #include "AliJCORRANTask.h"
  42 #include "AliAnalysisManager.h"
  43 #include "AliESDEvent.h"
  44 #include "AliMCEvent.h"
  45 #include "AliStack.h"
  46 #include "AliGenEventHeader.h"
  47 #include "AliGenCocktailEventHeader.h"
  48 #include "AliGenPythiaEventHeader.h"
  49 #include "AliInputEventHandler.h"
  50 #include "AliESDCaloCluster.h"
  51 #include "AliAODEvent.h"
  52 #include "AliAODHeader.h"
  53 #include "AliLog.h"
  54 #include "AliESDVertex.h"
  55 #include "AliESDtrack.h"
  56 #include "AliAODTrack.h"
  57 #include "AliAnalysisFilter.h"
  58 #include "AliESDtrackCuts.h"
  59 #include "AliPHOSGeoUtils.h"
  60 #include "AliEMCALGeoUtils.h"
  61 #include "AliESDtrackCuts.h"
  62
  63 #include "AliPhJTrackList.h"
  64 #include "AliPhJMCTrackList.h"
  65 #include "AliPhJPhotonList.h"
  66 #include "AliPhJHeaderList.h"
  67 #include "AliJTrack.h"
  68 #include "AliJPhoton.h"
  69 #include "AliJHeader.h"
  70 #include "AliAODTracklets.h"
  71 #include "AliMultiplicity.h"
  72 #include "JConst.h"
  73 #include "AliESDRun.h"
  74 #include "AliJRunHeader.h"
  75
  76
  77 //______________________________________________________________________________
  78 AliJCORRANTask::AliJCORRANTask() :
  79   AliAnalysisTaskSE("PWG4JCORRAN"),
  80   fInputFormat(0),
  81   fEsdTrackCuts(0), //FK//
  82   fDownscaling(0),//FK//
  83   fLowerCutOnLPMom(0),//FK//
  84   fLowerCutOnLeadingCaloClusterE(0),//FK//
  85   fTree(0x0),
  86   fTrackList(0x0),
  87   //fMCTrackList(0),
  88   fPhotonList(0x0),
  89   fHeaderList(0x0),
  90   fAliRunHeader(0x0),
  91   fQAList(0x0),
  92   fTrackQACuts(0x0),
  93   fTrackKineCuts(0x0),
  94   fPHOSGeom(0x0),
  95   fEMCALGeom(0x0)
  96 {
  97   //Default constructor
  98   fInputFormat = "ESD";
  99
 100   for(Int_t i=0;i<kRangeTriggerTableAlice;i++)   fActiveTriggers[i]=" ";
 101   for(Int_t i=0;i<kRangeTriggerTableJCorran;i++) fTriggerTableJCorran[i]=" ";
 102
 103   DefineInput (0, TChain::Class());
 104
 105 }
 106
 107 //______________________________________________________________________________
 108 AliJCORRANTask::AliJCORRANTask(const char *name, TString inputformat, AliESDtrackCuts* esdTrackCuts, Int_t downSc, Double_t lowLPmom, Double_t lowCaloE) :
 109   AliAnalysisTaskSE(name),
 110   fInputFormat(inputformat),
 111   fEsdTrackCuts(esdTrackCuts), //FK//
 112   fDownscaling(downSc),//FK//
 113   fLowerCutOnLPMom(lowLPmom),//FK//
 114   fLowerCutOnLeadingCaloClusterE(lowCaloE),//FK//
 115   fTree(0x0),
 116   fTrackList(0x0),
 117   //fMCTrackList(0),
 118   fPhotonList(0x0),
 119   fHeaderList(0x0),
 120   fAliRunHeader(0x0),
 121   fQAList(0x0),
 122   fTrackQACuts(0x0),
 123   fTrackKineCuts(0x0),
 124   fPHOSGeom(0x0),
 125   fEMCALGeom(0x0)
 126 {
 127   // Constructor.
 128   if(fDebug > 5) cout << "---- AliJCORRANTask Constructor ----"<<endl;
 129
 130   for(Int_t i=0;i<kRangeTriggerTableAlice;i++)   fActiveTriggers[i]=" ";
 131   for(Int_t i=0;i<kRangeTriggerTableJCorran;i++) fTriggerTableJCorran[i]=" ";
 132
 133   DefineInput (0, TChain::Class());
 134
 135   fTrackList    = new AliPhJTrackList(kALICE);
 136   //fMCTrackList  = new AliPhJMCTrackList(kALICE);
 137   fPhotonList   = new AliPhJPhotonList(kALICE);
 138   fHeaderList   = new AliPhJHeaderList(kALICE);
 139
 140   fAliRunHeader = new AliJRunHeader();
 141
 142   fPHOSGeom  = new AliPHOSGeoUtils("PHOSgeo") ;
 143   fEMCALGeom = new AliEMCALGeoUtils("EMCAL_COMPLETE");
 144
 145
 146   DefineOutput(1, TTree::Class());
 147   DefineOutput(2, TList::Class());
 148
 149 }
 150
 151 //____________________________________________________________________________
 152 AliJCORRANTask::AliJCORRANTask(const AliJCORRANTask& ap) :
 153   AliAnalysisTaskSE(ap.GetName()),
 154   fInputFormat(ap.fInputFormat),
 155   fEsdTrackCuts(ap.fEsdTrackCuts), //FK//
 156   fDownscaling(ap.fDownscaling),  //FK//
 157   fLowerCutOnLPMom(ap.fLowerCutOnLPMom), //FK//
 158   fLowerCutOnLeadingCaloClusterE(ap.fLowerCutOnLeadingCaloClusterE),  //FK//
 159   fTree(ap.fTree),
 160   fTrackList(ap.fTrackList),
 161   //fMCTrackList(ap.fMCTrackList),
 162   fPhotonList(ap.fPhotonList),
 163   fHeaderList(ap.fHeaderList),
 164   fAliRunHeader(ap.fAliRunHeader),
 165   fQAList(ap.fQAList),
 166   fTrackQACuts(ap.fTrackQACuts),
 167   fTrackKineCuts(ap.fTrackKineCuts),
 168   fPHOSGeom(ap.fPHOSGeom),
 169   fEMCALGeom(ap.fEMCALGeom)
 170 {
 171   // cpy ctor
 172 }
 173
 174 //_____________________________________________________________________________
 175 AliJCORRANTask& AliJCORRANTask::operator = (const AliJCORRANTask& ap)
 176 {
 177 // assignment operator
 178
 179   this->~AliJCORRANTask();
 180   new(this) AliJCORRANTask(ap);
 181   return *this;
 182 }
 183
 184 //______________________________________________________________________________
 185 AliJCORRANTask::~AliJCORRANTask()
 186 {
 187   // destructor
 188
 189   delete fTrackList;
 190   //delete fMCTrackList;
 191   delete fPhotonList;
 192   delete fHeaderList;
 193   delete fAliRunHeader;
 194   delete fTree;
 195
 196   delete fQAList;
 197   delete fTrackQACuts;
 198   delete fTrackKineCuts;
 199   if(fPHOSGeom)  delete fPHOSGeom  ;
 200   if(fEMCALGeom) delete fEMCALGeom ;
 201
 202 }
 203
 204 //________________________________________________________________________
 205 void AliJCORRANTask::UserCreateOutputObjects()
 206 {
 207   // create the jcorran outputs objects
 208   if(fDebug > 5) cout << "AliJCORRANTask UserCreateOutputObjects----------------------"<<endl;
 209   OpenFile(1) ;   // Will open the file for the object to be written at output #1
 210
 211   fTree = new TTree("T","ALICE JYU CORRAN");
 212   fTree->Branch("JTKT/JTrackList","AliPhJTrackList",&fTrackList,500000);
 213   //fTree->Branch("JTKT/JMCTrackList","AliPhJMCTrackList",&fMCTrackList,500000);
 214   fTree->Branch("JTKT/JPhotonList","AliPhJPhotonList",&fPhotonList,500000);
 215   fTree->Branch("JTKT/JHeaderList","AliPhJHeaderList",&fHeaderList,500000);
 216
 217   //QA histgrams
 218   // Offline Trigger QA
 219   OpenFile(2) ; // Will open the file for the object to be written at output #2
 220   TString trackQAlist = "MinNClusterTPC:";
 221           trackQAlist+= "MinNClustersITS:";
 222           trackQAlist+= "MaxChi2PerClusterTPC:";
 223           trackQAlist+= "MaxChi2PerClusterITS:";
 224           trackQAlist+= "RequireTPCRefit:";
 225           trackQAlist+= "RequireITSRefit:";
 226           trackQAlist+= "AcceptKinkDaughters:";
 227           trackQAlist+= "MaxCovMatrixDiag11:";
 228           trackQAlist+= "MaxCovMatrixDiag22:";
 229           trackQAlist+= "MaxCovMatrixDiag33:";
 230           trackQAlist+= "MaxCovMatrixDiag44:";
 231           trackQAlist+= "MaxCovMatrixDiag55:";
 232           trackQAlist+= "MaxNsigmaToVertex:";
 233           trackQAlist+= "RequireSigmaToVertex";
 234
 235   TString trackKinelist = "PMin:PMax:";
 236           trackKinelist+= "PtMin:PtMax:";
 237           trackKinelist+= "PxMin:PxMax:";
 238           trackKinelist+= "PyMin:PyMax:";
 239           trackKinelist+= "PzMin:PzMax:";
 240           trackKinelist+= "EtaMin:EtaMax:";
 241           trackKinelist+= "RapMin:RapMax";
 242
 243   fTrackQACuts    = new TNtuple("TrackQACuts", "track quality cuts applied on ESD",trackQAlist.Data());
 244   fTrackKineCuts  = new TNtuple("TrackKineCuts", "kinematic cuts applied on ESD", trackKinelist.Data());
 245   fQAList = new TList();
 246   fQAList->SetName("ESD QA List") ;
 247   fQAList->AddAt(fTrackQACuts,   0) ;
 248   fQAList->AddAt(fTrackKineCuts, 1) ;
 249
 250 }
 251
 252 //______________________________________________________________________________
 253 void AliJCORRANTask::UserExec(Option_t */*option*/)
 254 {
 255   // Processing of one event
 256   //if(fDebug > 5) cout << "------- AliJCORRANTask Exec-------"<<endl;
 257   if(!((Entry()-1)%100))
 258       AliInfo(Form(" Processing event # %lld",  Entry()));
 259
 260
 261
 262   static int runId=-1;
 263
 264   if(fInputFormat=="ESD"){
 265    //  if(fDebug > 5) cout <<"--------- Reading ESD --------"<< endl;
 266      AliESDEvent* esd = (AliESDEvent*)InputEvent();
 267
 268      //AliMCEvent* fmcEvent = MCEvent();
 269      //ReadMCTracks(fmcEvent);
 270
 271
 272      //=========== FILL AND STORE RUN HEADER   (ONLY ONCE PER RUN) =============
 273      if(esd->GetRunNumber() != runId){ //new run has started
 274
 275        const AliESDRun* esdRun = esd->GetESDRun();
 276
 277        for(Int_t triggerBit=0; triggerBit<kRangeTriggerTableAlice; triggerBit++){
 278          fActiveTriggers[triggerBit] = esdRun->GetTriggerClass(triggerBit);
 279        }
 280        runId = esd->GetRunNumber();
 281
 282        //Polarity of magnetic field in L3 solenoid
 283        Short_t l3MgFieldPolarity=0; // (LHC convention: +current -> +Bz)
 284        if(esd->GetCurrentL3() >0) l3MgFieldPolarity =  1;
 285        if(esd->GetCurrentL3() <0) l3MgFieldPolarity = -1;
 286
 287        //Create JCorran trigger mask mapping
 288        fTriggerTableJCorran[kMinBiasTriggerBitJCorran]="Minimum Bias";//minimum bias occupies first trigger bit
 289
 290        //=========== Fill Run header object ===============
 291        fAliRunHeader->SetRunNumber(runId);
 292        fAliRunHeader->SetActiveTriggersAlice(fActiveTriggers);
 293        fAliRunHeader->SetL3Field(l3MgFieldPolarity, esd->GetMagneticField());
 294        fAliRunHeader->SetActiveTriggersJCorran(fTriggerTableJCorran,kRangeTriggerTableJCorran);
 295
 296       //Store Run header
 297       (fTree->GetUserInfo())->Add(fAliRunHeader);
 298     }
 299
 300
 301     if( StoreDownscaledMinBiasEvent() || ContainsESDHighPtTrack(esd) || ContainsESDHighECaloClusters(esd) ){ //FK//
 302       //-------------- reset all the arrays -------------
 303       fTrackList->Reset();
 304       //fMCTrackList->Reset();
 305       fPhotonList->Reset();
 306       fHeaderList->Reset();
 307
 308       //store event only when it is downscaled min bias
 309       // or contais high pt hadron
 310       // or contains high energy cluster in EMCAL or PHOS
 311         ReadESDTracks(esd);
 312         ReadESDCaloClusters(esd);
 313         ReadESDHeader(esd);
 314
 315         fTree->Fill(); // fill the TTree
 316
 317         PostData(1, fTree);
 318     }
 319   }else{
 320
 321        AliAODEvent* aod;
 322        if(fInputFormat == "AODout") // reading from AOD output handler
 323          aod = AODEvent();
 324        else
 325          if(fInputFormat == "AODin") // reading from AOD input handler
 326             aod = (AliAODEvent*)InputEvent();
 327          else{
 328             cout << "Error: Not correct InputDataFormat especified " << endl;
 329             return;
 330          }
 331
 332      ReadAODTracks(aod);
 333      ReadAODCaloClusters(aod);
 334      ReadAODHeader(aod);
 335
 336      fTree->Fill(); // fill the TTree
 337
 338      PostData(1, fTree);
 339
 340   }
 341
 342
 343 }
 344
 345 //______________________________________________________________________________
 346 void AliJCORRANTask::Init()
 347 {
 348   // Intialisation of parameters
 349   AliInfo("Doing initialization") ;
 350
 351 }
 352
 353 //______________________________________________________________________________
 354 void AliJCORRANTask::Terminate(Option_t *)
 355 {
 356   // Processing when the event loop is ended
 357   fTree->Print();
 358   if(fInputFormat == "AODout") ReadFilter(); // change it to save this info also from AODin !!!!
 359
 360   ((AliJRunHeader *) (fTree->GetUserInfo())->First())->PrintOut();
 361
 362   cout<<"PWG4JCORRAN Analysis DONE !!"<<endl;
 363 }
 364
 365 //______________________________________________________________________________
 366 void AliJCORRANTask::ReadESDTracks(const AliESDEvent * esd)
 367 {
 368
 369     // Read the AliESDtrack and fill the list of AliJTrack containers
 370     Int_t nt = esd->GetNumberOfTracks();
 371    // if(fDebug < 5) cout << "ESD::NumberOfTracks = " << nt << endl;
 372     Float_t ptot, pt, eta;
 373     TVector3 p3(0,0,0);
 374     Short_t ntrk = 0;
 375     Double_t pid[10];
 376
 377     //loop over tracks
 378     for(Int_t it = 0; it < nt; it++) {
 379
 380         AliESDtrack *track = esd->GetTrack(it);
 381         if(! fEsdTrackCuts->IsSelected(track)) continue; //FK// apply loose selection criteria
 382
 383         UInt_t status = track->GetStatus();
 384
 385         Float_t impactXY, impactZ;
 386         track->GetImpactParameters(impactXY,impactZ);
 387         //------------ T P C ------------
 388         Float_t tpcDca[2], tpcCov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 389         track->GetImpactParametersTPC(tpcDca,tpcCov);
 390
 391         Float_t tpcDedx = track->GetTPCsignal();
 392         Int_t   tpcNcls = track->GetTPCNcls();
 393
 394         Int_t nClust         = track->GetTPCclusters(0);
 395         Int_t nFindableClust = track->GetTPCNclsF();
 396         Float_t tpcChi2PerCluster = 0.;
 397         if(nClust>0.) tpcChi2PerCluster = track->GetTPCchi2()/Float_t(nClust);
 398
 399         Float_t tpcClustPerFindClust = 0.;
 400         if(nFindableClust>0.) tpcClustPerFindClust = Float_t(nClust)/nFindableClust;
 401         //--------------------------------
 402
 403         Double_t mom[3];
 404         track->GetPxPyPz(mom);
 405         p3.SetXYZ(mom[0],mom[1],mom[2]);
 406         ptot = track->GetP();
 407         pt = p3.Pt();
 408         eta = p3.Eta();
 409         track->GetESDpid(pid);
 410
 411         Double_t extCov[15];
 412         track->GetExternalCovariance(extCov);
 413
 414         Double_t extDiaCov[5];
 415         extDiaCov[0]=extCov[0];
 416         extDiaCov[1]=extCov[2];
 417         extDiaCov[2]=extCov[5];
 418         extDiaCov[3]=extCov[9];
 419         extDiaCov[4]=extCov[14];
 420
 421         //create a new AliJTrack and fill the track info
 422         fTrackList->AddAliJTrack(ntrk);
 423         AliJTrack *ctrack = fTrackList->GetAliJTrack(ntrk);
 424
 425         ctrack->SetPtot(ptot);
 426         ctrack->SetPt(pt);
 427         ctrack->SetTheta(p3.Theta());
 428         ctrack->SetPhi(p3.Phi());
 429         ctrack->SetPID(pid);
 430         ctrack->SetFlavor(kHadron);
 431         ctrack->SetCharge(track->Charge());
 432         ctrack->ConvertAliPID();
 433         ctrack->SetEta(eta);
 434
 435         Int_t itof = track->GetTOFcluster(); // index of the assigned TOF cluster
 436         if(itof>=0) ctrack->SetTof(track->GetTOFsignal());
 437
 438         ctrack->SetTPCdEdx(tpcDedx);
 439         ctrack->SetTPCnClust(tpcNcls);
 440         ctrack->SetTPCDCAXY(tpcDca[0]);
 441         ctrack->SetTPCDCAZ(tpcDca[1]);
 442         ctrack->SetTPCClustPerFindClust(tpcClustPerFindClust);
 443         ctrack->SetTPCChi2PerClust(tpcChi2PerCluster);
 444         ctrack->SetImapactXY(impactXY);
 445         ctrack->SetImapactZ(impactZ);
 446
 447         ctrack->SetKinkIndex(track->GetKinkIndex(0));
 448         ctrack->SetStatus(status);
 449         ctrack->SetExternalDiaCovariance(extDiaCov);
 450
 451         fTrackList->SetNTracks(++ntrk);
 452
 453      } // end tracks loop
 454 }
 455
 456 //______________________________________________________________________________
 457 void AliJCORRANTask::ReadAODTracks(const AliAODEvent * aod)
 458 {
 459     // Read the AliAODtrack and fill the list of AliJTrack containers
 460     Int_t nt = aod->GetNumberOfTracks();
 461     if(fDebug < 5) cout << "AOD::NumberOfTracks = " << nt << endl;
 462     Short_t ntrk = 0;
 463     //loop over tracks
 464     for(Int_t it = 0; it < nt; it++) {
 465
 466         AliAODTrack *track = aod->GetTrack(it);
 467         if(track->GetType() != AliAODTrack::kPrimary) continue; // only primaries
 468         if(! fEsdTrackCuts->IsSelected(track)) continue; //apply loose selection criteria
 469         //create a new AliJTrack and fill the track info
 470         fTrackList->AddAliJTrack(ntrk);
 471         AliJTrack *ctrack = fTrackList->GetAliJTrack(ntrk);
 472
 473         ctrack->SetPt(track->Pt());
 474         ctrack->SetTheta(track->Theta());
 475         ctrack->SetPhi(track->Phi());
 476         ctrack->SetPID((Double_t*)track->PID());
 477         ctrack->SetFlavor(kHadron);
 478         ctrack->SetCharge(track->Charge());
 479         ctrack->SetChi2perNDF(track->Chi2perNDF());
 480         ctrack->SetChi2Trig(track->GetChi2MatchTrigger());
 481         ctrack->SetRecFlags(track->GetFlags());
 482
 483         fTrackList->SetNTracks(++ntrk);
 484
 485      } // end tracks loop
 486 }
 487
 488 //______________________________________________________________________________
 489 void AliJCORRANTask::ReadESDCaloClusters(const AliESDEvent* esd)
 490 {
 491   // Read the AliESDCaloClusters and fill the list of AliJPhoton containers
 492   Short_t nPhotons = 0;
 493   Int_t numberOfCaloClusters = esd->GetNumberOfCaloClusters() ;
 494   if(fDebug < 5) cout << "ESD::number of ESD caloclusters " << numberOfCaloClusters << endl;
 495   // loop over all the Calo Clusters
 496   for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 497     AliESDCaloCluster *caloCluster = esd->GetCaloCluster(icluster) ;
 498     if(!caloCluster) continue;
 499     if(caloCluster->GetTrackMatched()==-1){
 500       if(caloCluster->E()<0.2) continue;                  //FK//
 501       // we will not implement any PID cut here
 502       fPhotonList->AddAliJPhoton(nPhotons);
 503       AliJPhoton *pht = fPhotonList->GetAliJPhoton(nPhotons);
 504       pht->SetFlavor(kPhoton);
 505       pht->SetE(caloCluster->E());
 506       pht->SetChi2(caloCluster->GetClusterChi2());
 507       pht->SetPID(caloCluster->GetPid());
 508       Float_t pos[3]; caloCluster->GetPosition(pos) ;
 509       pht->SetX(pos[0]);
 510       pht->SetY(pos[1]);
 511       pht->SetZ(pos[2]);
 512       pht->SetPhi(atan2(pos[1],pos[0]));
 513       pht->SetTheta(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2]));
 514       pht->SetPt(caloCluster->E()*sin(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2])));
 515       pht->SetCharge(0);
 516       if(caloCluster->IsEMCAL()) pht->SetCaloType(AliJPhoton::kEMCALCalo);
 517       if(caloCluster->IsPHOS())  pht->SetCaloType(AliJPhoton::kPHOSCalo);
 518       pht->SetDistToBadChannel(caloCluster->GetDistanceToBadChannel());
 519       pht->SetDispersion(caloCluster->GetClusterDisp());
 520       pht->SetM20(caloCluster->GetM20());
 521       pht->SetM02(caloCluster->GetM02());
 522       pht->SetEmcCpvDist(caloCluster->GetEmcCpvDistance());
 523       pht->SetNCells(caloCluster->GetNCells());
 524       pht->SetCellsAbsId(caloCluster->GetCellsAbsId());
 525       Int_t imoduleID = GetSuperModuleNumber(caloCluster->IsEMCAL(), caloCluster->GetCellAbsId(0));
 526       pht->SetSuperModuleID(imoduleID);
 527   //    char buffer[20] = "photon";
 528  //     pht->PrintOut(buffer);
 529
 530       fPhotonList->SetNPhotons(++nPhotons);
 531     } // end if
 532   } //PHOS and EMCAL clusters
 533
 534 }
 535
 536 //______________________________________________________________________________
 537 void AliJCORRANTask::ReadAODCaloClusters(const AliAODEvent* aod)
 538 {
 539   // Read the AliAODCaloClusters and fill the list of AliJPhoton containers
 540   Int_t numberOfCaloClusters = aod->GetNCaloClusters() ;
 541   if(fDebug < 5) cout << "AOD::number of ESD caloclusters " << numberOfCaloClusters << endl;
 542   Short_t nPhotons = 0;
 543   for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 544        AliAODCaloCluster * caloCluster = aod->GetCaloCluster(icluster) ;
 545        if(!caloCluster) continue;
 546        if(caloCluster->GetNTracksMatched() > 0) continue;
 547       // we will not implement any PID cut here
 548       fPhotonList->AddAliJPhoton(nPhotons);
 549       AliJPhoton *pht = fPhotonList->GetAliJPhoton(nPhotons);
 550
 551       pht->SetE(caloCluster->E());
 552       pht->SetFlavor(kPhoton);
 553       pht->SetChi2(caloCluster->Chi2());
 554       pht->SetPID((Double_t*)caloCluster->PID());
 555       Float_t pos[3]; caloCluster->GetPosition(pos) ;
 556       pht->SetX(pos[0]);
 557       pht->SetY(pos[1]);
 558       pht->SetZ(pos[2]);
 559       pht->SetPhi(atan2(pos[1],pos[0]));
 560       pht->SetTheta(atan2(sqrt(pos[0]*pos[1]+pos[1]*pos[1]),pos[2]));
 561       pht->SetPt(caloCluster->E()*sin(atan2(sqrt(pos[0]*pos[0]+pos[1]*pos[1]),pos[2])));
 562       pht->SetCharge(0);
 563       if(caloCluster->IsEMCALCluster()) pht->SetCaloType(AliJPhoton::kEMCALCalo);
 564       if(caloCluster->IsPHOSCluster())  pht->SetCaloType(AliJPhoton::kPHOSCalo);
 565       pht->SetDistToBadChannel(caloCluster->GetDistToBadChannel());
 566       pht->SetDispersion(caloCluster->GetDispersion());
 567       pht->SetM20(caloCluster->GetM20());
 568       pht->SetM02(caloCluster->GetM02());
 569       pht->SetEmcCpvDist(caloCluster->GetEmcCpvDistance());
 570       pht->SetNCells(int(caloCluster->GetNCells()));
 571       pht->SetCellsAbsId(caloCluster->GetCellsAbsId());
 572       Int_t imoduleID = GetSuperModuleNumber(caloCluster->IsEMCALCluster(), caloCluster->GetCellAbsId(0));
 573       pht->SetSuperModuleID(imoduleID);
 574
 575       fPhotonList->SetNPhotons(++nPhotons);
 576
 577   } // clusters
 578 }
 579
 580 //______________________________________________________________________________
 581 void AliJCORRANTask::ReadESDHeader(const AliESDEvent *esd)
 582 {
 583     // Read the AliESDEvent and fill the list of AliJHeader containers
 584     Short_t nHeaders = 0;
 585     //create a header and fill it
 586     fHeaderList->AddAliJHeader(nHeaders);
 587     AliJHeader *hdr = fHeaderList->GetAliJHeader(nHeaders);
 588
 589     //cout << esd->GetRunNumber() <<"\t"<< esd->GetEventNumberInFile() << endl;
 590
 591     AliMultiplicity *fSPDMult =(AliMultiplicity *) esd->GetMultiplicity();
 592     hdr->SetSPDTrackletMult(fSPDMult->GetNumberOfTracklets());
 593
 594     hdr->SetEventID( esd->GetEventNumberInFile());
 595
 596     hdr->SetTriggerMaskAlice(esd->GetTriggerMask()); //ULong64_t
 597     hdr->SetTriggerMaskJCorran(ConvertTriggerMask(/*esd->GetTriggerMask()*/)); //UInt_t
 598
 599     const AliESDVertex * vtxESD = esd->GetVertex();
 600     hdr->SetZVertex(vtxESD->GetZv());
 601     hdr->SetZVertexErr(vtxESD->GetZRes());
 602
 603     fHeaderList->SetNHeaders(++nHeaders);
 604 }
 605
 606 //______________________________________________________________________________
 607 void AliJCORRANTask::ReadAODHeader(const AliAODEvent *aod)
 608 {
 609   //read AOD event header
 610     Short_t nHeaders = 0;
 611        //create a header and fill it
 612         fHeaderList->AddAliJHeader(nHeaders);
 613         AliJHeader *hdr = fHeaderList->GetAliJHeader(nHeaders);
 614
 615         //load aod event header
 616         AliAODHeader * aodh = aod->GetHeader();
 617
 618         hdr->SetCentrality(int(aodh->GetCentrality()));
 619
 620         hdr->SetTriggerMaskAlice(aodh->GetTriggerMask()); //ULong64_t
 621         hdr->SetEventType(aodh->GetEventType());
 622
 623         fHeaderList->SetNHeaders(++nHeaders);
 624 }
 625
 626 //______________________________________________________________________________
 627 void AliJCORRANTask::ReadFilter()
 628 {
 629    //read filter
 630    TList* olist = OutputTree()->GetUserInfo();
 631    AliAnalysisFilter* trackFilter = (AliAnalysisFilter*)olist->At(0);
 632
 633    if(!trackFilter) return;
 634
 635    TList* cutslist = trackFilter->GetCuts();
 636    //cout << "number of cuts " << cutslist->GetSize();
 637
 638    for(int ic = 0; ic< cutslist->GetSize(); ic++){
 639        AliESDtrackCuts* icuts = (AliESDtrackCuts*)cutslist->At(ic);
 640
 641        // read ESD track quality cuts
 642        Float_t trkcut[14];
 643        trkcut[0] = (Float_t)icuts->GetMinNClusterTPC();
 644        trkcut[1] = (Float_t)icuts->GetMinNClustersITS();
 645        trkcut[2] = (Float_t)icuts->GetMaxChi2PerClusterTPC();
 646        trkcut[3] = (Float_t)icuts->GetMaxChi2PerClusterITS();
 647        trkcut[4] = (Float_t)icuts->GetRequireTPCRefit();
 648        trkcut[5] = (Float_t)icuts->GetRequireITSRefit();
 649        trkcut[6] = (Float_t)icuts->GetAcceptKinkDaughters();
 650        icuts->GetMaxCovDiagonalElements(trkcut[7],trkcut[8],trkcut[9],trkcut[10],trkcut[11]);
 651        trkcut[12] = (Float_t)icuts->GetMaxNsigmaToVertex();
 652        trkcut[13] = (Float_t)icuts->GetRequireSigmaToVertex();
 653
 654        fTrackQACuts->Fill(trkcut);
 655
 656        // read track kinmatic cuts
 657        Float_t kinecut[14];
 658        icuts->GetPRange(kinecut[0], kinecut[1]);
 659        icuts->GetPtRange(kinecut[2], kinecut[3]);
 660        icuts->GetPxRange(kinecut[4], kinecut[5]);
 661        icuts->GetPyRange(kinecut[6], kinecut[7]);
 662        icuts->GetPzRange(kinecut[8], kinecut[9]);
 663        icuts->GetEtaRange(kinecut[10], kinecut[11]);
 664        icuts->GetRapRange(kinecut[12], kinecut[13]);
 665
 666        fTrackKineCuts->Fill(kinecut);
 667
 668
 669    } // end cuts loop
 670
 671    PostData(2,fQAList);
 672 }
 673 //______________________________________________________________________________
 674 Int_t AliJCORRANTask::GetSuperModuleNumber(bool isemcal, Int_t absId)
 675 {
 676   //get super module number
 677   if(isemcal){
 678     return fEMCALGeom->GetSuperModuleNumber(absId) ;
 679   } else {
 680     Int_t    relId[4];
 681     if ( absId >= 0) {
 682       fPHOSGeom->AbsToRelNumbering(absId,relId);
 683       return relId[0]-1;
 684     } else return -1;
 685   }//PHOS
 686
 687   return -1;
 688 }
 689
 690 //_____________________________________________________________________________
 691
 692 UInt_t AliJCORRANTask::ConvertTriggerMask(/*Long64_t alicetriggermask*/){
 693  //convert alice trigger mask to jcorran trigger mask
 694   UInt_t triggerMaskJC=0;
 695
 696   Bool_t isSelected = kTRUE;
 697   isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 698
 699   if(isSelected){ //tag collision candidates
 700     triggerMaskJC += (1<<kMinBiasTriggerBitJCorran);
 701   }
 702
 703   return triggerMaskJC;
 704 }
 705
 706
 707 //______________________________________________________________________________
 708 bool AliJCORRANTask::StoreDownscaledMinBiasEvent(){
 709
 710   bool isThisEventToBeStored = kFALSE;
 711   static Long_t evtNumber=0;
 712
 713   if( ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() ){
 714     //collision candidate
 715
 716     if(evtNumber% fDownscaling ==0){ isThisEventToBeStored = kTRUE;  } //store every 50th collision candidate event
 717     evtNumber++;
 718   }
 719   return isThisEventToBeStored;
 720 }
 721 //______________________________________________________________________________
 722 bool AliJCORRANTask::ContainsESDHighPtTrack(const AliESDEvent* esd){
 723   bool isThisEventToBeStored = kFALSE;
 724
 725   Int_t nt = esd->GetNumberOfTracks();
 726
 727   for(Int_t it = 0; it < nt; it++) {
 728     AliESDtrack *track = esd->GetTrack(it);
 729     //Does event contain high pt particle above 2 GeV
 730     //which fulfills loose track selection criteria
 731     if(track->Pt() > fLowerCutOnLPMom && fEsdTrackCuts->IsSelected(track)){
 732       isThisEventToBeStored = kTRUE;
 733       break;
 734     }
 735   }
 736
 737   return isThisEventToBeStored;
 738 }
 739
 740 //______________________________________________________________________________
 741 bool AliJCORRANTask::ContainsESDHighECaloClusters(const AliESDEvent* esd){
 742   bool isThisEventToBeStored = kFALSE;
 743
 744   Int_t numberOfCaloClusters = esd->GetNumberOfCaloClusters() ;
 745   // loop over all the Calo Clusters
 746   for(Int_t icluster = 0 ; icluster < numberOfCaloClusters ; icluster++) {
 747     AliESDCaloCluster *caloCluster = esd->GetCaloCluster(icluster) ;
 748     if(!caloCluster) continue;
 749     if(caloCluster->GetTrackMatched()==-1){
 750       //sotre calo clusters above 2 GeV
 751       if( caloCluster->E() > fLowerCutOnLeadingCaloClusterE){
 752         isThisEventToBeStored = kTRUE;
 753         break;
 754       }
 755     }
 756   }
 757   return isThisEventToBeStored;
 758 }
 759
 760
 761
 762
 763
 764
 765 /*
 766
 767
 768 void AliJCORRANTask::ReadMCTracks(AliMCEvent *fMC)
 769 {
 770   AliGenEventHeader* genHeader = fMC->GenEventHeader();
 771   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 772
 773   Double_t ptHard = 0;
 774   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
 775
 776     nTrials = pythiaGenHeader->Trials();
 777     ptHard  = pythiaGenHeader->GetPtHard();
 778     //cout << nTrials <<"\t"<< ptHard << endl;
 779
 780     AliStack *stack = fMC->Stack();
 781     Int_t np    = fMC->GetNumberOfTracks();
 782     Int_t nprim = stack->GetNtrack();
 783     if(np!=nprim) cout << "GetNumberOfTracks = "<< np <<"\t, stack = "<< nprim << endl;
 784     Short_t ntrack = 0;
 785     for(Int_t itrack = 0; itrack < nprim; itrack++){
 786       //create a new JMCTrack and fill the track info
 787       fMCTrackList->AddJMCTrack(ntrack);
 788       JMCTrack *ctrack = fMCTrackList->GetTrack(ntrack);
 789       TParticle *ipart = stack->Particle(itrack);
 790       Int_t icode= ipart->GetPdgCode();
 791       Float_t pt = ipart->Pt();
 792       Float_t eta = ipart->Eta();
 793       Float_t phi = ipart->Phi();
 794       Int_t   status  = ipart->GetStatusCode();
 795       Bool_t  isprim =  ipart->IsPrimary();
 796       // fill base class
 797       ctrack->SetE(ipart->Energy());
 798       ctrack->SetTheta(ipart->Theta());
 799       ctrack->SetPhi(phi);
 800       ctrack->SetPt(pt);
 801       // MC
 802       ctrack->SetPtHard(ptHard);
 803       ctrack->SetPdgCode(icode);
 804       ctrack->SetStatusCode(status);
 805       ctrack->SetIsPrimary(isprim);
 806       ctrack->SetProductionVertex(ipart->Vx(),ipart->Vx(),ipart->Vz());
 807
 808       //bool isInc = (status ==  1 && icode ==  22); //Inclusive
 809       bool ispi0 = (status == 11 && icode == 111); //kPizero
 810       bool isDgamma = (status == 6 || status == 7) && icode == 22; // Direct photon
 811       bool inPHOS  = (ispi0||isDgamma)&&fabs(eta)<0.12;
 812       bool inEMCAL = (ispi0||isDgamma)&&fabs(eta)<0.7;
 813       bool inTPC   = fabs(eta)<0.9;
 814       ctrack->SetMother(0,ipart->GetFirstMother());
 815       ctrack->SetMother(1,ipart->GetSecondMother());
 816       ctrack->SetDaughter(0,ipart->GetFirstDaughter());
 817       ctrack->SetDaughter(1,ipart->GetLastDaughter());
 818       ctrack->SetIsInPHOS(inPHOS);
 819       ctrack->SetIsInEMCAL(inEMCAL);
 820       ctrack->SetIsInTPC(inTPC);
 821
 822       fMCTrackList->SetNTracks(++ntrack);
 823     }// loop for al primary tracks
 824 }
 825 */