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