update to master versions
[u/mrichter/AliRoot.git] / PWGCF / Correlations / JCORRAN / AliJCORRAN.cxx
CommitLineData
66be7134 1/**************************************************************************
2 * Copyright(c) 1998-2014, 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// Comment describing what this class does needed!
17
18// jcorran main class
19// used in local and grid execution
20// blah
21// blah
22// blah
23
24#include <TH1D.h>
25#include "AliJCORRAN.h"
26
27#include "AliJTrackCounter.h"
28#include <TClonesArray.h>
29
30#include "AliJCard.h"
31#include "AliJHistos.h"
32#include "AliJCorrelations.h"
33#include "AliJEventPool.h"
34#include "AliJDataManager.h"
35
36#include "AliJEventHeader.h"
37#include "AliJRunHeader.h"
38#include "AliJTrack.h"
39#include "AliJPhoton.h"
40#include "AliJMCTrack.h"
41#include "AliJConst.h"
42
43
44
45
46#include "AliJEfficiency.h"
47#include <iostream>
48
49ClassImp(AliJCORRAN)
50
51 AliJCORRAN::AliJCORRAN() :
52 TObject(),
53 fExecLocal(kTRUE),
54 fFirstEvent(kTRUE),
55 fjtrigg((particleType)-100),
56 fjassoc((particleType)-100),
57 fcard(0),
58 finputFile(0),
59 fInclusiveFile(""),
60 fevt(0),
61 fhistos(0),
62 fcorrelations(0),
63 fphotonPool(0),
64 fassocPool(0),
65 fphotonList(0),
66 fchargedHadronList(0),
67 fpizeroList(0),
68 ftriggList(0),
69 fassocList(0),
70 fpairList(0),
71 fpairCounterList(0),
72 finputList(0),
73 fdmg(0),
74 feventHeader(0),
75 frunHeader(0),
76 fnumberEvents(0),
77 fieout(0),
78 fEventCounter(0),
79 fcent(0),
80 fncBin(0),
81 fnPttBin(0),
82 fbTriggCorrel(0),
83 fbLPCorrel(0),
84 fbLPpairCorrel(0),
85 fTrackEtaRange(0),
86 flowerPtAssocBoarder(0),
87 fCentMultLow(0),
88 fCentMultHigh(0),
89 fEventBC(0),
90 fSQRTS(0),
91 fEfficiency(0),
92 fRunTable(0),
93 fIsolationR(0),
94 fHadronSelectionCut(0)
95{
96 // constructor
97}
98
99AliJCORRAN::AliJCORRAN(Bool_t execLocal) :
100 TObject(),
101 fExecLocal(execLocal),
102 fFirstEvent(kTRUE),
103 fjtrigg((particleType)-100),
104 fjassoc((particleType)-100),
105 fcard(0),
106 finputFile(0),
107 fInclusiveFile(""),
108 fevt(0),
109 fhistos(0),
110 fcorrelations(0),
111 fphotonPool(0),
112 fassocPool(0),
113 fphotonList(0),
114 fchargedHadronList(0),
115 fpizeroList(0),
116 ftriggList(0),
117 fassocList(0),
118 fpairList(0),
119 fpairCounterList(0),
120 finputList(0),
121 fdmg(0),
122 feventHeader(0),
123 frunHeader(0),
124 fnumberEvents(0),
125 fieout(0),
126 fEventCounter(0),
127 fcent(0),
128 fncBin(0),
129 fnPttBin(0),
130 fbTriggCorrel(0),
131 fbLPCorrel(0),
132 fbLPpairCorrel(0),
133 fTrackEtaRange(0),
134 flowerPtAssocBoarder(0),
135 fCentMultLow(0),
136 fCentMultHigh(0),
137 fEventBC(0),
138 fSQRTS(0),
139 fEfficiency(0),
140 fRunTable(0),
141 fIsolationR(0),
142 fHadronSelectionCut(0)
143{
144 // constructor
145}
146
147AliJCORRAN::~AliJCORRAN(){
148 // destructor
149}
150
151AliJCORRAN::AliJCORRAN(const AliJCORRAN& obj) :
152 TObject(),
153 fExecLocal(obj.fExecLocal),
154 fFirstEvent(obj.fFirstEvent),
155 fjtrigg(obj.fjtrigg),
156 fjassoc(obj.fjassoc),
157 fcard(obj.fcard),
158 finputFile(obj.finputFile),
159 fInclusiveFile(obj.fInclusiveFile),
160 fevt(obj.fevt),
161 fhistos(obj.fhistos),
162 fcorrelations(obj.fcorrelations),
163 fphotonPool(obj.fphotonPool),
164 fassocPool(obj.fassocPool),
165 fphotonList(obj.fphotonList),
166 fchargedHadronList(obj.fchargedHadronList),
167 fpizeroList(obj.fpizeroList),
168 ftriggList(obj.ftriggList),
169 fassocList(obj.fassocList),
170 fpairList(obj.fpairList),
171 fpairCounterList(obj.fpairCounterList),
172 finputList(obj.finputList),
173 fdmg(obj.fdmg),
174 feventHeader(obj.feventHeader),
175 frunHeader(obj.frunHeader),
176 fnumberEvents(obj.fnumberEvents),
177 fieout(obj.fieout),
178 fEventCounter(obj.fEventCounter),
179 fcent(obj.fcent),
180 fncBin(obj.fncBin),
181 fnPttBin(obj.fnPttBin),
182 fbTriggCorrel(obj.fbTriggCorrel),
183 fbLPCorrel(obj.fbLPCorrel),
184 fbLPpairCorrel(obj.fbLPpairCorrel),
185 fTrackEtaRange(obj.fTrackEtaRange),
186 flowerPtAssocBoarder(obj.flowerPtAssocBoarder),
187 fCentMultLow(obj.fCentMultLow),
188 fCentMultHigh(obj.fCentMultHigh),
189 fEventBC(obj.fEventBC),
190 fSQRTS(obj.fSQRTS),
191 fEfficiency(obj.fEfficiency),
192 fRunTable(obj.fRunTable),
193 fIsolationR(obj.fIsolationR),
194 fHadronSelectionCut(obj.fHadronSelectionCut)
195{
196 // copy constructor
197 JUNUSED(obj);
198}
199
200AliJCORRAN& AliJCORRAN::operator=(const AliJCORRAN& obj){
201 // copy constructor
202 JUNUSED(obj);
203 return *this;
204}
205
206
207void AliJCORRAN::Initialize() const{
208 // init
209
210}
211
212void AliJCORRAN::UserCreateOutputObjects(){
213 // local init
214
215
216 cout << "jcorran user create output objects ----------------" << endl;
217
218 fHadronSelectionCut =int ( fcard->Get("HadronSelectionCut"));
219 fIsolationR = fcard->Get("IsolationR");
220
221 fhistos = new AliJHistos( fcard );
222 fhistos->CreateEventTrackHistos();
223 fhistos->CreateAzimuthCorrHistos();
224 fhistos->CreateIAAMoons();
225 fhistos->CreateXEHistos();
226 fhistos->CreateXtHistos();
227 fhistos->CreatePairPtCosThetaStar();
228
229 fhistos->fHMG->Print();
230
231 fEventBC = (Int_t)(fcard->Get( "eventBC" ));
232 fSQRTS = 0.;
233
234
235 //TODO: inclusive fhistos shipping along
236
237 fcorrelations = new AliJCorrelations( fcard, fhistos);
238 cout<<endl<< " -----" <<endl;
239 if( fInclusiveFile.Length() ) {
240 fhistos->ReadInclusiveHistos(fInclusiveFile);
241 fcorrelations->SetSampligInclusive(); //kperp background and triangle. Default is flat
242 cout<<"Sampling kperp and triangle from " << fInclusiveFile <<endl;
243 } else cout << "Sampling kperp and triangle from flat" <<endl;
244 cout<< " -----" <<endl <<endl;
245
246 //==================================
247
248 //cout<<kParticleTypeStrName[kPhoton]<<" "<<kParticleTypeStrName[fjtrigg]<<endl;
249 // EventPool for Mixing
5ca97880 250 fphotonPool = new AliJEventPool( fcard, fhistos, fcorrelations, kJPhoton); // for pi0 mass
66be7134 251 fassocPool = new AliJEventPool( fcard, fhistos, fcorrelations, fjassoc);
252
5ca97880 253 fphotonList = new TClonesArray(kParticleProtoType[kJPhoton],1500);
66be7134 254 // TClonesArray *cellList = new TClonesArray("AliJCaloCell",1500);
5ca97880 255 fchargedHadronList = new TClonesArray(kParticleProtoType[kJHadron],1500);
256 fpizeroList = new TClonesArray(kParticleProtoType[kJPizero],1500);
66be7134 257 ftriggList = new TClonesArray(kParticleProtoType[fjtrigg],1500);
258 fassocList = new TClonesArray(kParticleProtoType[fjassoc],1500);
259 fpairList = new TClonesArray(kParticleProtoType[fjtrigg],1500);
260 fpairCounterList = new TClonesArray("AliJTrackCounter",1500);
261 finputList = NULL;
262 //TClonesArray *isolPizeroList = new TClonesArray("AliPhJPiZero",1500);
263
264 fdmg = new AliJDataManager(fcard, fhistos, fcorrelations, fExecLocal);
265 fdmg->SetExecLocal( fExecLocal );
266
267 //==== Read the Data files =====
268 if( fExecLocal ){
269 fdmg->ChainInputStream(finputFile);
270 // TODO: run header is not supposed to be here
271 // doesn't work fSQRTS = 2.*frunHeader->GetBeamEnergy();
272
273 // for grid running, numberEvents is filled by the encapsulating
274 // grid task, which has access to the input handlers and can
275 // extract event numbers out of it
276 fnumberEvents = fdmg->GetNEvents();
277 fieout = fnumberEvents/20;
278 frunHeader = fdmg->GetRunHeader();
279 cout<<"RunID = "<<frunHeader->GetRunNumber()<< " Looping over "<<fnumberEvents<<" events"<<endl;
280
281 } else {
282 fdmg->SetRunHeader( frunHeader );
283 frunHeader = fdmg->GetRunHeader();
284 }
285
286 //==== Efficiency ====
287 fEfficiency = new AliJEfficiency;
288 fEfficiency->SetMode( fcard->Get("EfficiencyMode") ); // 0:NoEff, 1:Period 2:RunNum 3:Auto
4431deef 289 if(fExecLocal) {
290 fEfficiency->SetDataPath("/mnt/flustre/alice/taxi_jcorran/2013/EFF/data"); // Efficiency root file location local or alien
291 } else {
292 fEfficiency->SetDataPath("alien:///alice/cern.ch/user/d/djkim/legotrain/efficieny/data"); // Efficiency root file location local or alien
293 }
66be7134 294
295 //-------------------------------------------------------------------
296 fEventCounter=0;
297
298 fcent = -1;
299 fncBin = fcard->GetNoOfBins( kCentrType );
300 fnPttBin = fcard->GetNoOfBins( kTriggType );
301
302 fbTriggCorrel = fcard->Get("CorrelationType")==0;
303 fbLPCorrel = fcard->Get("CorrelationType")==1;
304 fbLPpairCorrel = fcard->Get("CorrelationType")==2;
305 fTrackEtaRange = fcard->Get("EtaRange");
306 flowerPtAssocBoarder = fcard->GetBinBorder(kAssocType, 0);
307
308 //==== Limit the no of track for each fcentrality ======
309 fCentMultLow = new TF1("fCentMultLow","[0]*sqrt([1]-x)+[2]", 0, 75);
310 fCentMultHigh = new TF1("fCentMultHigh","[0]*sqrt([1]-x)+[2]", 0, 90);
311
312 // fast parameter load
313
4431deef 314
66be7134 315 fhistos->fHMG->WriteConfig();
316 fFirstEvent = kTRUE;
317 fevt = -1;
318
319 cout << "end of jcorran user create output objects ----------------" << endl;
320
321}
322
323void AliJCORRAN::UserExec(){
324 // event loop
325 fevt++;
326
327 // TODO Add train mode if( !fExecLocal && fFirstEvent ){
328 if( 0 && !fExecLocal && fFirstEvent ){
329 fdmg->ChainInputStream("");
330 fieout = fnumberEvents/20;
331 if (fieout<1) fieout=1;
332 cout<<"RunID = "<<frunHeader->GetRunNumber()<< " Looping over "<<fnumberEvents<<" events"<<endl;
333
334 //==== Efficiency ====
335 // TODO run can be different in a job?
336 fhistos->fhEventPerRun->Fill(fnumberEvents>0?log(fnumberEvents):1);
337 fFirstEvent = kFALSE;
338 }
339
340 // TODO
341 if( fExecLocal ) {
342 if(fevt % fieout == 0) cout << fevt << "\t" << int(float(fevt)/fnumberEvents*100) << "%" << endl ;
343 }
344
345 if( fFirstEvent ) {
346 //==== RunTable
347 fRunTable = & AliJRunTable::GetSpecialInstance();
348 fRunTable->SetRunNumber( frunHeader->GetRunNumber() );
349 fSQRTS = fRunTable->GetBeamEnergy(fRunTable->GetPeriod());
350 cout << "sqrts = "<< fSQRTS << ",runnumber = "<< frunHeader->GetRunNumber() << endl;
351 fEfficiency->SetRunNumber( frunHeader->GetRunNumber() );
66be7134 352 fEfficiency->Load();
353 fFirstEvent = kFALSE;
354 }
355
356 if(fRunTable->IsHeavyIon()){
357 fCentMultLow->SetParameters( fcard->Get("CentMultCutLow",0), fcard->Get("CentMultCutLow",1), fcard->Get("CentMultCutLow",2) );
358 fCentMultHigh->SetParameters(fcard->Get("CentMultCutHigh",0), fcard->Get("CentMultCutHigh",1), fcard->Get("CentMultCutHigh",2) );
359 //fCentMultLow->Print();
360 //fCentMultHigh->Print();
361 }
362
363 fdmg->LoadEvent(fevt); // to be here to load EventHeader
364 fhistos->fhEvents->Fill( 0 );
365
366 if(!fdmg->IsGoodEvent()) return; // Vertex cut applied in IsGoodEvent and histo saved there too
367
368 feventHeader = fdmg->GetEventHeader();
369 double zVert = feventHeader->GetZVertex();
370 //----------------------------------------------------------
371
372 UInt_t triggermaskJCorran = feventHeader->GetTriggerMaskJCorran();
373 //cout << fevt <<"\t"<< zVert << "\t"<< triggermaskJCorran << endl;
374
375 if( fdmg->IsSelectedTrigger((int) triggermaskJCorran))
376 fhistos->fhEvents->Fill( 5 );
377
378 // select only some BC%4
379 if( feventHeader->GetBunchCrossNumber() % 4 != fEventBC &&
380 fEventBC > -1 )
381 return;
382
383 if( fdmg->IsSelectedTrigger((int) triggermaskJCorran))
384 fhistos->fhEvents->Fill( 6 );
385
386 if(fRunTable->IsHeavyIon() || fRunTable->IsPA()){
387 fcent = feventHeader->GetCentrality();
388 }else {
389 fcent = 1; //ntracks;
390 }
391
392 //cout<<"evt="<<fevt <<" c="<< fcent <<endl;
393 int cBin = fcard->GetBin(kCentrType, fcent);
394 if(cBin<0) return;
395
396 if( fdmg->IsSelectedTrigger((int) triggermaskJCorran))
397 fhistos->fhEvents->Fill( 7 );
398
399 int zBin = fcard->GetBin(kZVertType, zVert); //should be alway >0; checked in fdmg->IsGoodEvent()
400
401 // do not fill MB in case of MB mixing
402 if( fdmg->IsSelectedTrigger((int) triggermaskJCorran))
403 fhistos->fhZVert[cBin]->Fill(zVert);
404
405 //temporary fevent skip
406 //Trigger to be selected from the JCorran trigger mask is specified in the fCard
407 // if(! fdmg->IsSelectedTrigger((int) triggermaskJCorran))
408 // continue;
409
410 fEventCounter++;
411
412 //==== QA Event
413 fhistos->fhV0AMult[cBin]->Fill( feventHeader->GetV0AMult() );
414
415 //------------------------------------------------------------------
416 // Triggers and associated
417 //----------------------ooooo---------------------------------------
418
5ca97880 419 if(fjtrigg==kJPizero || fjassoc==kJPizero || fjtrigg==kJPhoton || fjassoc==kJPhoton){
66be7134 420 } // pizero || photon
5ca97880 421 if(fjtrigg==kJHadron || fjassoc==kJHadron){
66be7134 422 fchargedHadronList->Clear();
5ca97880 423 fdmg->RegisterList(fchargedHadronList, NULL, cBin, zBin, kJHadron);
66be7134 424 // apply efficiencies
425
426 for( int i = 0; i < fchargedHadronList->GetEntries(); i++ ){
427
428 AliJBaseTrack *triggTr = (AliJBaseTrack*)fchargedHadronList->At(i);
429 double ptt = triggTr->Pt();
430
431 double effCorr = 1./fEfficiency->GetCorrection(ptt, fHadronSelectionCut, fcent); // here you generate warning if ptt>30
432 fhistos->fhTrackingEfficiency[cBin]->Fill( ptt, 1./effCorr );
433 triggTr->SetTrackEff( 1./effCorr );
434 }
435 }
436
437 //---- assign input list ----
5ca97880 438 if(fjtrigg==kJPizero) finputList = fpizeroList;
439 else if(fjtrigg==kJHadron) finputList = fchargedHadronList;
440 else if(fjtrigg==kJPhoton) finputList = fphotonList;
66be7134 441 int noAllTriggTracks = finputList->GetEntries();
442 int noAllChargedTracks = fchargedHadronList->GetEntries();
443 fhistos->fhChargedMult[cBin]->Fill(noAllChargedTracks);
444 fhistos->fhChargedMultCent->Fill(fcent, noAllChargedTracks>0 ? log(noAllChargedTracks) : 0);
445
446 //---------------------------------------------
447 //-- Check fcentrality outliers ---
448 //-- and count enevents only here ---
449 //---------------------------------------------
450 if(fRunTable->IsHeavyIon() && fCentMultLow->GetParameter(0) >0 ){
451 double logMult = noAllChargedTracks>0 ? log(noAllChargedTracks) : 0 ;
452 if( fcent<fCentMultLow->GetParameter(1) && logMult <fCentMultLow->Eval(fcent) ) return;
453 if( fCentMultHigh->Eval(fcent) < logMult) return;
454 }
455 fhistos->fhChargedMultCut[cBin]->Fill(noAllChargedTracks);
456 fhistos->fhCentr->Fill(fcent);
457 fhistos->fhiCentr->Fill(cBin);
458
459 // ------------------------------
460 // --- calculate e-b-e vn
461 // ------------------------------
462 if(fRunTable->IsHeavyIon()){
463 double vdelta[2] = {0};
464 int vdeltaNorm = 0;
465 for(int it1=0; it1<noAllTriggTracks-1; it1++){
466 AliJBaseTrack *ftk1 = (AliJBaseTrack*)finputList->At(it1);
467 if(ftk1->Pt()<flowerPtAssocBoarder) continue;
468 for(int it2=it1+1; it2<noAllTriggTracks; it2++){
469 AliJBaseTrack *ftk2 = (AliJBaseTrack*)finputList->At(it2);
470 if(ftk2->Pt()<flowerPtAssocBoarder) continue;
471 if(fabs(ftk1->Eta() - ftk2->Eta())<1.0) continue;
472 double fdphi = ftk1->DeltaPhi(*ftk2);
473 fhistos->fhVN[cBin]->Fill(fdphi);
474 vdelta[0] += cos(2*fdphi);
475 vdelta[1] += cos(3*fdphi);
476 //cout<< ftk1->Pt() <<" "<< ftk2->Pt() <<" "<< fabs(ftk1->Eta() - ftk2->Eta()) <<" "<< 2*fdphi <<" "<< 3*fdphi <<endl;
477 vdeltaNorm++;
478 }
479 }
480 vdelta[0] = vdeltaNorm>0 ? vdelta[0]/vdeltaNorm : 0;
481 vdelta[1] = vdeltaNorm>0 ? vdelta[1]/vdeltaNorm : 0;
482 fhistos->fhVdelta2[cBin]->Fill(vdelta[0]*100);
483 fhistos->fhVdelta3[cBin]->Fill(vdelta[1]*100);
484 if(vdelta[0]>0) fhistos->fpV2->Fill( fcent, sqrt(vdelta[0])*100 );
485 if(vdelta[1]>0) {
486 fhistos->fpV3->Fill( fcent, sqrt(vdelta[1])*100 );
487 fcard->SetEventV3kv(vdelta[1]);
488 }
489 fhistos->fpVdeltaNorm->Fill( fcent, vdeltaNorm );
490 }
491
492 //----------------------------------------------------
493 //----- Generate trigg list and find LP --
494 //----- Fiducial cut should be used in AliJCorrelations--
495 //----------------------------------------------------
496 AliJTrackCounter *lpTrackCounter = new AliJTrackCounter(), *lpPairCounter = new AliJTrackCounter();
497 AliJBaseTrack *lPTr = NULL;
498 int noTriggs=0;
499 ftriggList->Clear();
500 for(int itrack=0; itrack<noAllTriggTracks; itrack++){
501 AliJBaseTrack *triggTr = (AliJBaseTrack*)finputList->At(itrack);
502 triggTr->SetTriggBin( fcard->GetBin(kTriggType, triggTr->Pt()) );
503
504 double ptt = triggTr->Pt();
505 double etat = triggTr->Eta();
506 //TODO iCut == 0;
507
508 double effCorr = 1.0/triggTr->GetTrackEff();
509
510 if( ptt>flowerPtAssocBoarder ){
511 //FK//double effCorr = 1./fcard->TrackEfficiency(ptt, fcent); // here you generate warning if ptt>30
512 //double effCorr = 1./fcard->TrackEfficiency(ptt, etat, cBin); // here you generate warning if ptt>30
513 fhistos->fhChargedPt[cBin]->Fill(ptt, effCorr );
514 fhistos->fhChargedPtNoCorr[cBin]->Fill( ptt );
515 fhistos->fhChargedEta->Fill(triggTr->Eta(), effCorr);
516 //fhistos->fhChargedPtJacek[cBin]->Fill(ptt, effCorr );
517 fhistos->fhChargedPtJacek[cBin]->Fill(ptt, ptt>0 ? 1./ptt*effCorr : 0); //One CANNOT do 1/ptt here!! First unfold.
518 if( -0.8<etat && etat<-0.2) fhistos->fhChargedPtJacekEta[cBin][0]->Fill(ptt, ptt>0 ? 1./ptt*effCorr : 0);
519 if( -0.2<etat && etat< 0.3) fhistos->fhChargedPtJacekEta[cBin][1]->Fill(ptt, ptt>0 ? 1./ptt*effCorr : 0);
520 if( 0.3<etat && etat< 0.8) fhistos->fhChargedPtJacekEta[cBin][2]->Fill(ptt, ptt>0 ? 1./ptt*effCorr : 0);
521 fhistos->fhChargedPtFiete->Fill(ptt, effCorr );
522 }
523
524 if( !triggTr->IsInTriggerBin() ) continue;
525 //FK//triggTr->SetTrackEff( fcard->TrackEfficiency(triggTr->Pt(), fcent) );
526 int iptt = triggTr->GetTriggBin();
527 fhistos->fhIphiTrigg[cBin][iptt]->Fill( triggTr->Phi(), effCorr);
528 fhistos->fhIetaTrigg[cBin][iptt]->Fill( triggTr->Eta(), effCorr);
529
530 if( ptt > lpTrackCounter->Pt() ) {
531 lpTrackCounter->Store(noTriggs, ptt, iptt);
532 lPTr = triggTr;
533 }
534 //cout <<"1 ptt="<< triggTr->Pt() <<" bin="<< triggTr->GetTriggBin() << " st=" << triggTr->GetStatus() << " trigg eff=" << triggTr->GetTrackEff() << endl;
535 new ((*ftriggList)[noTriggs++]) AliJBaseTrack(*triggTr);
536 fhistos->fhTriggMult[cBin][iptt]->Fill(noAllTriggTracks);
537 }
538
539 //----------------------------------------------------
540 //----- Find sum of two leading particles ------------
541 //----------------------------------------------------
542 if(fbLPpairCorrel){
543 int noPairs=0;
544 fpairList->Clear();
545 fpairCounterList->Clear();
546 for(int ii=0;ii<noTriggs-1;ii++){
547 AliJBaseTrack *triggTr1 = (AliJBaseTrack*)ftriggList->At(ii);
548 for(int jj=ii+1;jj<noTriggs;jj++){
549 AliJBaseTrack *triggTr2 = (AliJBaseTrack*)ftriggList->At(jj);
550 TLorentzVector lVPair = triggTr1->GetLorentzVector()+triggTr2->GetLorentzVector();
551 double fdphi = DeltaPhi(triggTr1->Phi(), triggTr2->Phi());
552 new ((*fpairList)[noPairs]) AliJBaseTrack(lVPair);
553 new ((*fpairCounterList)[noPairs++])
554 AliJTrackCounter( triggTr1->GetID(), triggTr2->GetID(), ii, jj, fdphi );
555 }
556 }
557
558 // ----- Find the Leading Pair --------------------------
559 AliJBaseTrack *pairTr = NULL;
560 for(int ii=0;ii<noPairs;ii++){
561 pairTr = (AliJBaseTrack*)fpairList->At(ii);
562 AliJTrackCounter *pairCounter = (AliJTrackCounter*)fpairCounterList->At(ii);
563 //cout<<pairTr->Pt()<<endl; pairCounter->Print();
564 if(pairTr->Pt() > lpPairCounter->Pt() && fabs(pairCounter->GetPairDPhi())<1.0) {
565 int iptt = fcard->GetBin(kTriggType, pairTr->Pt());
566 lpPairCounter = pairCounter;
567 lpPairCounter->Store(ii, pairTr->Pt(), iptt);
568 }
569 }
570 if(lpPairCounter->Exists()){
5ca97880 571 pairTr->SetParticleType(kJHadron);
66be7134 572 //double effCorr = 1./fcard->TrackEfficiency(lpTrackCounter->GetLPpt());
573 fhistos->fhLPpairPt->Fill(pairTr->Pt());
574 }
575 }
576
577 //--------------------------------------------------
578 //--- Generate assoc list and pool ---
579 //--------------------------------------------------
580 fassocList->Clear();
581 int noAssocs=0;
582 double sumPtAroundLP = 0;
5ca97880 583 if(fjassoc==kJPizero) finputList = fpizeroList;
584 else if(fjassoc==kJHadron) finputList = fchargedHadronList;
585 else if(fjassoc==kJPhoton) finputList = fphotonList;
66be7134 586
587 int noAllAssocTracks = finputList->GetEntries();
588
589
590 for(int itrack=0;itrack<noAllAssocTracks; itrack++){
591
592 AliJBaseTrack *assocTr = (AliJBaseTrack*)finputList->At(itrack);
593 assocTr->SetAssocBin( fcard->GetBin(kAssocType, assocTr->Pt()) );
594
595 if(assocTr->IsInAssocBin()){
596
597 int ipta = assocTr->GetAssocBin();
598 double effCorr = 1.0/assocTr->GetTrackEff();
599 fhistos->fhIphiAssoc[cBin][ipta]->Fill( assocTr->Phi(), effCorr);
600 fhistos->fhIetaAssoc[cBin][ipta]->Fill( assocTr->Eta(), effCorr);
601 new ((*fassocList)[noAssocs++]) AliJBaseTrack(*assocTr);
602 }
603 // determine an activity in the cone around trigger
604 if(lpTrackCounter->Exists() && (lPTr->GetID()!=assocTr->GetID()) && (0.5 < assocTr->Pt())){
605 double dPhi = DeltaPhi( assocTr->Phi(), lPTr->Phi() );
606 double dEta = assocTr->Eta() - lPTr->Eta();
607 if(fIsolationR > sqrt(dPhi*dPhi+dEta*dEta)) sumPtAroundLP += assocTr->Pt();//FK//
608 }
609 }
610
611 FillXtHistos(finputList, lpTrackCounter);
612
613 //-----------------------------------------------
614 // Set the isolation flag
615 //-----------------------------------------------
616 if( lpTrackCounter->Exists() ){
617 //double effCorr = 1./fcard->TrackEfficiency(lpTrackCounter->Pt(), fcent);
618 //TODO ??? AGAIN AGAIN???
619 double effCorr = 1./fEfficiency->GetCorrection(lpTrackCounter->Pt(), fHadronSelectionCut, fcent );
620 fhistos->fhLPpt->Fill(lpTrackCounter->Pt(), effCorr);
621 fhistos->fhLPeta->Fill(lPTr->Eta(), effCorr);
622 fhistos->fhBkgActivity[lpTrackCounter->GetPtBin()]->Fill(sumPtAroundLP/lpTrackCounter->Pt());
623
624 if( sumPtAroundLP/lpTrackCounter->Pt()< fcard->GetCutOnBkgActivity() && ///relative activity
625 ( fabs(lPTr->Eta()) < (fTrackEtaRange - fIsolationR) ) ){ //fiducial cut
626
627 fhistos->fhIsolatedLPpt->Fill(lpTrackCounter->Pt(), effCorr );
628
629 AliJBaseTrack* lpTrigger = (AliJBaseTrack*) ftriggList->At(lpTrackCounter->GetIndex());
630 lpTrigger->SetIsIsolated(1);
631 }
632 }
633
634 fhistos->fhAssocMult->Fill(noAssocs);
635 //cout<"Triggs = "<<<noTriggs<<" Assoc = "<<noAssocs<<" all = "<<noAllTracks<<endl;
636
637 //============================================================
638 //there is obviously a problem when you do "fixed" correlation.
639 //In this case the noTriggs==0 condition is never fullfilled
640 //============================================================
641 //if(noTriggs==0 || ((fjtrigg==kPizero) && (fjassoc==kPizero)) )
642 //if(leadingPt<1.5) //should be fixed
643 //if(noTriggs==0 && noAssocs>0 ){}
644 if(noAssocs>0 ) fassocPool->AcceptList(fassocList, fcent, zVert, noAllChargedTracks, fevt);
645
646 //------------------------------------------------------------------
647 // Do the Correlation
648 //----------------------ooooo---------------------------------------
649 int noTriggTracs=-1;
650 noTriggTracs = fbTriggCorrel ? noTriggs : 1;
651 if(fbLPCorrel && !lpTrackCounter->Exists()) return;
652 if(fbLPpairCorrel && !lpPairCounter->Exists()) return;
653 AliJBaseTrack *triggTr = NULL;
654
655 for(int ii=0;ii<noTriggTracs;ii++){ // trigger loop
656 if (fbTriggCorrel) triggTr = (AliJBaseTrack*)ftriggList->At(ii);
657 if (fbLPCorrel) triggTr = (AliJBaseTrack*)ftriggList->At(lpTrackCounter->GetIndex());
658 if (fbLPpairCorrel) triggTr = (AliJBaseTrack*)fpairList->At(lpPairCounter->GetIndex());
659 //for(int ii=0;ii<noIsolPizero;ii++){ // trigger loop }
660 //AliJBaseTrack *triggTr = (AliJBaseTrack*)isolPizeroList->At(ii);
661 double ptt = triggTr->Pt();
662 int iptt = triggTr->GetTriggBin();
663 if(iptt<0) {
664 cout<<"Not registered trigger ! I better stop here." <<endl;
665 exit(-1);
666 }
667 double effCorr = 1.0/triggTr->GetTrackEff();
668 fhistos->fhTriggPtBin[cBin][zBin][iptt]->Fill(ptt, effCorr);//inclusive
669
670 if(triggTr->GetIsIsolated()>0) fhistos->fhTriggPtBinIsolTrigg[kReal][cBin][iptt]->Fill(ptt, effCorr);
671
672 for(int jj=0;jj<noAssocs;jj++){ // assoc loop
673 AliJBaseTrack *assocTr = (AliJBaseTrack*)fassocList->At(jj);
674 //assocTr->PrintOut("assoc track");
675 if(fbLPpairCorrel &&
676 (assocTr->GetID()==lpPairCounter->GetPairTrackID(0) ||
677 assocTr->GetID()==lpPairCounter->GetPairTrackID(1)) ) continue;
678 //-------------------------------------------------------------
679 fcorrelations->FillAzimuthHistos(kReal, cBin, zBin, triggTr, assocTr);
680 //-------------------------------------------------------------
681 } // end assoc loop
682 } // end trigg loop
683 // == mix trigg with assoc
684 if (fbLPpairCorrel)
685 fassocPool->Mix(fpairList, kAzimuthFill, fcent, zVert, noAllChargedTracks, fevt);
686 else
687 fassocPool->Mix(ftriggList, kAzimuthFill, fcent, zVert, noAllChargedTracks, fevt);
688
689 //--------------------------------------------------------------
690 // End of the Correlation
691 //--------------------------------------------------------------
692
693
694
695}
696
697void AliJCORRAN::Terminate() {
698 // termination
699
700 /* TODO
701 for (int hic = 0;hic < fcard->GetNoOfBins(kCentrType);hic++){
702 ScaleNotEquidistantHisto( fhistos->fhChargedPt[hic], 1);
703 ScaleNotEquidistantHisto( fhistos->fhChargedPtNoCorr[hic], 1);
704 ScaleNotEquidistantHisto( fhistos->fhChargedPtJacek[hic], 1);
705 }
706 ScaleNotEquidistantHisto( fhistos->fhLPpt, 1);
707 ScaleNotEquidistantHisto( fhistos->fhLPpairPt, 1);
708 ScaleNotEquidistantHisto( fhistos->fhIsolatedLPpt, 1);
709 ScaleNotEquidistantHisto( fhistos->fhChargedPtFiete, 1);
710 */
711
712 // cout<<"MB's="<<noMB<<" "<<" ERT's="<<noERT<<endl;
713 fcorrelations->PrintOut();
714 fassocPool->PrintOut();
715 fEfficiency->Write();
716 fhistos->fHMG->Print();
717 //fhistos->fHMG->Write();
718}
719
720particleType AliJCORRAN::GetParticleType(char *inchar){
721 // part type
722 for(int i=0;i<kNumberOfParticleTypes;i++) {
723 //cout<<"<"<<inchar<<"> <"<<particleTypeStr[i]<<"> "<<strcmp(inchar,particleTypeStr[i])<<" "<<(particleType)i<<endl;
724 if(strcmp(inchar,kParticleTypeStrName[i])==0) return (particleType)i;
725 }
726 std::cout<<"Particle type <"<<inchar<<"> not recognized"<<std::endl;
727 exit(1);
728}
729
730double AliJCORRAN::DeltaPhi(double phi1, double phi2) {
731 // dphi
732 double res = atan2(sin(phi1-phi2), cos(phi1-phi2));
733 //return res>-kJPi/3.0 ? res : kJTwoPi+res ;
734 return res>-kJPi*9./20. ? res : kJTwoPi+res ;
735}
736
737void AliJCORRAN::ScaleNotEquidistantHisto(TH1D *hid, const double sc=1){
738 // histo scaler
739 for(int i=1;i<= hid->GetNbinsX();i++){
740 hid->SetBinContent(i,hid->GetBinContent(i)*sc/hid->GetBinWidth(i));
741 hid->SetBinError(i,hid->GetBinError(i)*sc/hid->GetBinWidth(i));
742 }
743}
744
745// Fill xT histograms (Esko)
746void AliJCORRAN::FillXtHistos(TClonesArray *inputList, AliJTrackCounter *lpTrackCounter){
747 // Here should be a comment describing this method
748 JUNUSED(lpTrackCounter);
749
750 enum xTtype { kInclusive=0, kIsolated=1, kIsolatedLP=2} ;
751 double lowerPtCut = 0.2;
752 double isolationR = 0.4;
753 double fCutOnBkgActivity = 0.10;
754 int cBin = fcard->GetBin(kCentrType, fcent);
755 double sqrts = fSQRTS;
756 int noAllTriggTracks = inputList->GetEntries();
757 //cout << "Entering Esko's analysis loop" << endl;
758 //cout << "Sqrts = " << sqrts << " pT = " << ptt << " xT = " << xtt << endl;
759 for(int itrack=0; itrack<noAllTriggTracks; itrack++){
760 AliJBaseTrack *triggTr = (AliJBaseTrack*)inputList->At(itrack);
761 double sumPtAroundTrigger = 0;
762 double ptt = triggTr->Pt();
763 double xtt = 2.0 * ptt / (1.0 * sqrts);
764 double effCorr = 1.0/triggTr->GetTrackEff();
765 fhistos->fhPtForXt[kInclusive][cBin]->Fill(ptt, ptt>0 ? 1./ptt*effCorr : 0);
766 fhistos->fhXt[kInclusive][cBin]->Fill(xtt, effCorr );
767 fhistos->fhXtWeighted[kInclusive][cBin]->Fill(xtt, effCorr*1.0/xtt );
768
769 for(int jj=0;jj<inputList->GetEntriesFast();jj++){ // assoc loop
770 AliJBaseTrack *assocTr = (AliJBaseTrack*)inputList->At(jj);
771 if(triggTr->GetID()==assocTr->GetID()) continue;
772 double pta = assocTr->Pt();
773 // determine the activity in the cone around trigger (Esko)
774 if( lowerPtCut < pta ){
775 double dPhi = DeltaPhi( assocTr->Phi(), triggTr->Phi() );
776 double dEta = assocTr->Eta() - triggTr->Eta();
777 if(isolationR > sqrt(dPhi*dPhi+dEta*dEta)) sumPtAroundTrigger += assocTr->Pt();//FK//
778 }
779 //cout << " Assoc number " << assocCounter++ << endl;
780 }
781 // If pT sum is below the limit, fill to the isolated histos
782 if( sumPtAroundTrigger/ptt < fCutOnBkgActivity && ///relative activity
783 ( fabs(triggTr->Eta()) < (fTrackEtaRange - isolationR) ) ){ //fiducial cut
784 // pT and xT
785 // kInclusive kIsolated, kLPIsolated
786 fhistos->fhPtForXt[kIsolated][cBin]->Fill(ptt,effCorr*1.0/ptt);
787 fhistos->fhXt[kIsolated][cBin]->Fill(xtt, effCorr );
788 fhistos->fhXtWeighted[kIsolated][cBin]->Fill(xtt,effCorr*1.0/xtt);
789 }
790 //cout<<"pT sum around trigger = " << sumPtAroundTrigger << endl;
791 }
792
793} // end FillXtHistos