update to master versions
[u/mrichter/AliRoot.git] / PWGCF / Correlations / JCORRAN / AliJCORRAN.cxx
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
49 ClassImp(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
99 AliJCORRAN::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
147 AliJCORRAN::~AliJCORRAN(){
148     // destructor
149 }
150
151 AliJCORRAN::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
200 AliJCORRAN& AliJCORRAN::operator=(const AliJCORRAN& obj){
201     // copy constructor
202     JUNUSED(obj);
203     return *this;
204 }
205
206
207 void AliJCORRAN::Initialize() const{
208     // init
209
210 }
211
212 void 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
250     fphotonPool  = new AliJEventPool( fcard, fhistos, fcorrelations, kJPhoton);  // for pi0 mass
251     fassocPool   = new AliJEventPool( fcard, fhistos, fcorrelations, fjassoc);
252
253     fphotonList = new TClonesArray(kParticleProtoType[kJPhoton],1500);      
254     //     TClonesArray *cellList = new TClonesArray("AliJCaloCell",1500);      
255     fchargedHadronList  = new TClonesArray(kParticleProtoType[kJHadron],1500);      
256     fpizeroList = new TClonesArray(kParticleProtoType[kJPizero],1500);      
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
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         }
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
314
315         fhistos->fHMG->WriteConfig();
316         fFirstEvent = kTRUE;
317         fevt = -1;
318
319         cout << "end of jcorran user create output objects ----------------" << endl;
320
321 }
322
323 void 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() );
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
419         if(fjtrigg==kJPizero || fjassoc==kJPizero || fjtrigg==kJPhoton || fjassoc==kJPhoton){
420         } // pizero || photon
421         if(fjtrigg==kJHadron || fjassoc==kJHadron){
422                 fchargedHadronList->Clear();
423                 fdmg->RegisterList(fchargedHadronList, NULL, cBin, zBin, kJHadron);
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 ---- 
438         if(fjtrigg==kJPizero)      finputList = fpizeroList;  
439         else if(fjtrigg==kJHadron) finputList = fchargedHadronList;
440         else if(fjtrigg==kJPhoton) finputList = fphotonList;
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()){
571                         pairTr->SetParticleType(kJHadron);
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;
583         if(fjassoc==kJPizero) finputList = fpizeroList;  
584         else if(fjassoc==kJHadron) finputList = fchargedHadronList;
585         else if(fjassoc==kJPhoton) finputList = fphotonList;
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
697 void 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
720 particleType  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
730 double 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
737 void 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)
746 void 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