76257f278c983cc71f180533246d33d6108ab769
[u/mrichter/AliRoot.git] / 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 "AliESDtrackCuts.h"
60 #include "AliPHOSGeoUtils.h"
61 #include "AliEMCALGeoUtils.h"
62 #include "AliESDtrackCuts.h"
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 //______________________________________________________________________________
79 AliJCORRANTask::AliJCORRANTask() :   
80   AliAnalysisTaskSE("PWG4JCORRAN"),
81   fInputFormat(0),
82   fEsdTrackCuts(0), 
83   fDownscaling(1),
84   fLowerCutOnLPMom(0),
85   fLowerCutOnLeadingCaloClusterE(0), 
86   fLowerCutOnCaloClusterE(0.2),
87   fIsRealOrMC(0),
88   fAODName("jcorran.root"),
89   fTrackList(0x0),
90   fMCTrackList(0x0),
91   fPhotonList(0x0),
92   fHeaderList(0x0),
93   fAliRunHeader(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) : 
109   AliAnalysisTaskSE(name), 
110   fInputFormat(inputformat),  
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"),
118   fTrackList(0x0),
119   fMCTrackList(0x0),
120   fPhotonList(0x0),
121   fHeaderList(0x0),
122   fAliRunHeader(0x0),
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
134  
135
136 }
137
138 //____________________________________________________________________________
139 AliJCORRANTask::AliJCORRANTask(const AliJCORRANTask& ap) :
140   AliAnalysisTaskSE(ap.GetName()), 
141   fInputFormat(ap.fInputFormat),
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),
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   fEMCALGeom(ap.fEMCALGeom)
156
157   // cpy ctor
158 }
159
160 //_____________________________________________________________________________
161 AliJCORRANTask& AliJCORRANTask::operator= (const AliJCORRANTask& ap)
162 {
163 // assignment operator
164
165   this->~AliJCORRANTask();
166   new(this) AliJCORRANTask(ap);
167   return *this;
168 }
169
170 //______________________________________________________________________________
171 AliJCORRANTask::~AliJCORRANTask()
172 {
173   // destructor 
174   
175   delete fTrackList;
176   if(fMCTrackList) delete fMCTrackList;
177   delete fPhotonList;
178   delete fHeaderList;
179   delete fAliRunHeader;
180   if(fPHOSGeom)  delete fPHOSGeom;
181   if(fEMCALGeom) delete fEMCALGeom;
182
183 }
184
185 //________________________________________________________________________
186 void AliJCORRANTask::UserCreateOutputObjects()
187 {  
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  
199   // create the jcorran output deltaAOD
200   //if(fDebug > 5) cout << "AliJCORRANTask UserCreateOutputObjects----------------------"<<endl;
201   
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   }
216   
217 }
218
219 //______________________________________________________________________________
220 void 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())); 
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   }
236
237   fTrackList->Reset();
238   fMCTrackList->Reset();
239   fPhotonList->Reset();
240   fHeaderList->Reset();
241  
242
243   static int runId=-1;
244  
245   if(fInputFormat=="ESD"){
246    //  if(fDebug > 5) cout <<"--------- Reading ESD --------"<< endl; 
247      AliESDEvent* esd = (AliESDEvent*)InputEvent();
248
249      AliMCEvent* mcEvent = NULL; 
250      if(fIsRealOrMC)  mcEvent = MCEvent();
251
252
253      //=========== FILL AND STORE RUN HEADER   (ONLY ONCE PER RUN) =============
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
277        //Store Run header
278        (OutputTree()->GetUserInfo())->Add(fAliRunHeader);  //FK// 
279     }
280
281     if(storeEvent){ 
282       //-------------- reset all the arrays -------------
283          //store event only when it is downscaled min bias
284       // or contais high pt hadron
285       // or contains high energy cluster in EMCAL or PHOS
286       ReadESDTracks(esd);
287       ReadESDCaloClusters(esd);
288       ReadESDHeader(esd);
289       if(fIsRealOrMC) ReadMCTracks(mcEvent);
290     }
291   }else if( fInputFormat == "AODout" || fInputFormat == "AODin") {
292   
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 -------------
302          
303       ReadAODTracks(aod);
304       ReadAODCaloClusters(aod);
305       ReadAODHeader(aod);
306     }
307     
308   }else{
309     cout << "Error: Not correct InputDataFormat especified " << endl;
310     return;
311   }
312 }
313
314 //______________________________________________________________________________
315 void AliJCORRANTask::Init()
316 {
317   // Intialisation of parameters
318   AliInfo("Doing initialization") ; 
319
320 }
321
322 //______________________________________________________________________________
323 void AliJCORRANTask::Terminate(Option_t *)
324 {
325   // Processing when the event loop is ended
326   OutputTree()->Print(); 
327   if(fInputFormat == "AODout") ReadFilter(); // change it to save this info also from AODin !!!! 
328
329   ((AliJRunHeader *) (OutputTree()->GetUserInfo())->First())->PrintOut();  
330
331   cout<<"PWG4JCORRAN Analysis DONE !!"<<endl; 
332
333
334 }
335
336 //______________________________________________________________________________
337 void AliJCORRANTask::ReadESDTracks(const AliESDEvent * esd)
338 {
339
340     // Read the AliESDtrack and fill the list of AliJTrack containers
341     Int_t nt = esd->GetNumberOfTracks();
342    // if(fDebug < 5) cout << "ESD::NumberOfTracks = " << nt << endl;
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);
352         if(! fEsdTrackCuts->IsSelected(track)) continue; //apply quality selection criteria
353
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();
367         Float_t tpcChi2PerCluster = 0.;
368         if(nClust>0.) tpcChi2PerCluster = track->GetTPCchi2()/Float_t(nClust);
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
391       //  Int_t itsLabel = track->GetITSLabel(); //FK//
392       //  Int_t tpcLabel = track->GetTPCLabel(); //FK//   
393
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);
403         ctrack->SetFlavor(kNone);//kHadron);
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
424       //  ctrack->SetITSLabel(itsLabel);//FK//
425       //  ctrack->SetTPCLabel(tpcLabel);//FK//
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(track->GetType() != AliAODTrack::kPrimary) continue; // only primaries 
445         if(! fEsdTrackCuts->IsSelected(track)) continue; //apply loose selection criteria
446         //create a new AliJTrack and fill the track info
447         fTrackList->AddAliJTrack(ntrk);
448         AliJTrack *ctrack = fTrackList->GetAliJTrack(ntrk);
449
450         ctrack->SetPtot(track->P());
451         ctrack->SetPt(track->Pt());
452         ctrack->SetTheta(track->Theta());
453         ctrack->SetPhi(track->Phi());
454         ctrack->SetEta(track->Eta());
455         ctrack->SetPID((Double_t*)track->PID());
456         ctrack->SetFlavor(kNone); //kHadron);
457         ctrack->SetCharge(track->Charge());
458         ctrack->SetChi2perNDF(track->Chi2perNDF());
459         ctrack->SetChi2Trig(track->GetChi2MatchTrigger());
460         ctrack->SetRecFlags(track->GetFlags());
461       
462        
463        // ctrack->SetITSLabel(track->GetLabel());//FK//?
464       //  ctrack->SetTPCLabel(track->GetLabel());//FK//?
465         fTrackList->SetNTracks(++ntrk);
466
467      } // end tracks loop
468 }
469
470 //______________________________________________________________________________
471 void 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){
482       if(caloCluster->E()<fLowerCutOnCaloClusterE) continue;                  //FK//
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);
509       
510       fPhotonList->SetNPhotons(++nPhotons);
511     } // end if 
512   } //PHOS and EMCAL clusters
513
514 }
515
516 //______________________________________________________________________________
517 void 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;
527       // we will not implement any PID cut here      
528       fPhotonList->AddAliJPhoton(nPhotons);
529       AliJPhoton *pht = fPhotonList->GetAliJPhoton(nPhotons);
530       
531       pht->SetE(caloCluster->E());
532       pht->SetFlavor(kPhoton);
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 //______________________________________________________________________________
561 void 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);
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 //______________________________________________________________________________
585 void AliJCORRANTask::ReadAODHeader(const AliAODEvent *aod)
586 {
587   //read AOD event header
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());
601         
602   fHeaderList->SetNHeaders(++nHeaders);
603 }
604
605 //______________________________________________________________________________
606 Int_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
624 UInt_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 //______________________________________________________________________________
640 bool 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 }
653
654 //______________________________________________________________________________
655 bool AliJCORRANTask::ContainsESDHighPtTrack(){
656
657   bool isThisEventToBeStored = kFALSE; //initialize return value
658
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       }
692     }
693   }
694
695   return isThisEventToBeStored;
696 }
697
698 //______________________________________________________________________________
699 bool AliJCORRANTask::ContainsESDHighECaloClusters(){
700   bool isThisEventToBeStored = kFALSE; //initialize return value
701
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
736       if( caloCluster->E() > fLowerCutOnLeadingCaloClusterE){
737         isThisEventToBeStored = kTRUE;
738         break; 
739       }
740     }
741   }
742
743   return isThisEventToBeStored;
744 }
745
746
747 //______________________________________________________________________________
748
749 void 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
756
757   //nTrials = pythiaGenHeader->Trials();
758   //ptHard  = pythiaGenHeader->GetPtHard();
759
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;
766
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     }
773
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);
807       
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
828 bool 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;
836 }
837
838 //______________________________________________________________________________
839
840
841
842
843
844