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