PWG2/SPECTRA/AliAnalysisTaskCentral.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 //  ------------------------------------------
  17 //  analysis task for azimuthal isotropic
  18 //  expansion in highly central collisions
  19 //  author: Cristian Andrei
  20 //          acristian@niham.nipne.ro
  21 //  ------------------------------------------
  22
  23 #include "TChain.h"
  24 #include "TTree.h"
  25 #include "TH1F.h"
  26 #include "TList.h"
  27 #include "TObjArray.h"
  28 #include "TString.h"
  29 #include "TFile.h"
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliMCEvent.h"
  33 #include "AliMCEventHandler.h"
  34 #include "AliESDEvent.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliESDtrackCuts.h"
  37 #include "AliCFContainer.h"
  38 #include "AliAnalysisDataContainer.h"
  39 #include "AliAnalysisDataSlot.h"
  40
  41 #include "AliAnalysisCentralCutMC.h"
  42 #include "AliAnalysisCentralCutESD.h"
  43 #include "AliAnalysisCentralCutEvtMC.h"
  44 #include "AliAnalysisCentralCutEvtESD.h"
  45 #include "AliAnalysisCentralExtrapolate.h"
  46 #include "AliAnalysisTaskCentral.h"
  47
  48
  49 ClassImp(AliAnalysisTaskCentral)
  50
  51
  52 //________________________________________________________________________
  53 AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)
  54   : AliAnalysisTask(name, "")
  55   ,fESD(0), fMC(0)
  56   ,fNoEvt(0)
  57   ,fCFContainerPi(0)
  58   ,fCFContainerK(0)
  59   ,fCFContainerP(0)
  60   ,fSim(kFALSE)
  61   ,fOutList(NULL)
  62
  63 {
  64   // Constructor
  65         printf("AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)\n");
  66
  67   // Define input and output slots here
  68         DefineInput(0, TChain::Class());
  69
  70         DefineOutput(0, TList::Class());
  71
  72
  73         for(Int_t i=0; i<10; i++){
  74                 fCutsList[i] = 0;
  75     }
  76
  77         InitCuts(); //initialize the analysis specific cuts
  78     if (!fCutsList) {
  79                 Printf("ERROR: fCutsList not available");
  80                 return;
  81     }
  82
  83         printf("AliAnalysisTaskCentral::Constructor\n");
  84 }
  85
  86
  87 //________________________________________________________________________
  88 AliAnalysisTaskCentral::~AliAnalysisTaskCentral()
  89 {
  90 // Destructor
  91 // Delete the created objects
  92
  93         if(fESD) delete fESD;
  94         if(fMC) delete fMC;
  95
  96         if(fCutsList){
  97         delete [] fCutsList;
  98         }
  99
 100         if(fNoEvt) delete fNoEvt;
 101
 102         if(fCFContainerPi) delete fCFContainerPi;
 103         if(fCFContainerPi) delete fCFContainerK;
 104         if(fCFContainerPi) delete fCFContainerP;
 105
 106         if(fOutList) delete fOutList;
 107
 108         printf("AliAnalysisTaskCentral::Destructor\n");
 109 }
 110
 111 //________________________________________________________________________
 112 void AliAnalysisTaskCentral::ConnectInputData(Option_t *) {
 113 // get the event from the input chain
 114
 115     TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 116         printf("tree(%p)\n", (void*)tree);
 117
 118     if (!tree) {
 119         Printf("ERROR: Could not read chain from input slot 0");
 120     } else {
 121         AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 122         if (!esdH) {
 123             Printf("ERROR: Could not get ESDInputHandler");
 124         }
 125         else fESD = esdH->GetEvent();
 126
 127         AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 128         if (!mcH) {
 129             Printf("ERROR: Could not get MCInputHandler");
 130         }
 131         else fMC = mcH->MCEvent();
 132     }
 133 }
 134
 135 //________________________________________________________________________
 136 void AliAnalysisTaskCentral::CreateOutputObjects(){
 137 //Creates the output list
 138
 139 //      fNoEvt = new TH1F(Form("%s_NoEvt",GetName()), "Number of processed events", 1, 0.0, 1.0);
 140
 141         fNoEvt = new TH1F("TaskCentral_NoEvt", "Number of processed events", 1, 0.0, 1.0);
 142 //set here the:
 143     //min and max variable values (to do: move these to a separate cuts class)
 144     const Double_t ptmin  = 0.0 ; //bins in pt
 145     const Double_t ptmax  = 5.0 ; //bins in y
 146     const Double_t etamin  = -0.9 ; //bins in pt
 147     const Double_t etamax  = 0.9 ; //bins in y
 148     //apropiate number of bins for histos and CF Container
 149     const Int_t nbinpt  = 25 ; //bins in pt
 150     const Int_t nbineta  = 25 ; //bins in y
 151
 152 //Correction Framework CONTAINER DEFINITION
 153     //the sensitive variables, their indices
 154         UInt_t ipt = 0;
 155         UInt_t ieta  = 1;
 156     //Setting up the container grid...
 157     const UInt_t nstep = 2 ; //number of selection steps MC & ESD
 158     const Int_t nvar = 2 ; //number of variables on the grid:pt,eta
 159
 160     //arrays for the number of bins in each dimension
 161     //nbin is defined above =
 162     Int_t iBin[nvar];
 163     iBin[0]=nbinpt;
 164     iBin[1]=nbineta;
 165
 166     //arrays for lower bounds :
 167     Double_t *binLim1=new Double_t[nbinpt+1];
 168     Double_t *binLim2=new Double_t[nbineta+1];
 169
 170     //values for bin lower bounds
 171     for(Int_t i=0; i<=nbinpt; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/nbinpt*(Double_t)i ;
 172     for(Int_t i=0; i<=nbineta; i++) binLim2[i]=(Double_t)etamin  + (etamax-etamin)  /nbineta*(Double_t)i ;
 173
 174         //container creation
 175         fCFContainerPi = new AliCFContainer("TaskCentral_CFCont_Pi","container for tracks",nstep,nvar,iBin);
 176         fCFContainerK = new AliCFContainer("TaskCentral_CFCont_K","container for tracks",nstep,nvar,iBin);
 177         fCFContainerP = new AliCFContainer("TaskCentral_CFCont_P","container for tracks",nstep,nvar,iBin);
 178
 179     //setting the bin limits
 180     fCFContainerPi -> SetBinLimits(ipt,binLim1);
 181     fCFContainerPi -> SetBinLimits(ieta,binLim2);
 182
 183         fCFContainerK -> SetBinLimits(ipt,binLim1);
 184         fCFContainerK -> SetBinLimits(ieta,binLim2);
 185
 186         fCFContainerP -> SetBinLimits(ipt,binLim1);
 187         fCFContainerP -> SetBinLimits(ieta,binLim2);
 188
 189 //outout list creation
 190         fOutList = new TList();
 191         fOutList->Add(fNoEvt);
 192         fOutList->Add(fCFContainerPi);
 193         fOutList->Add(fCFContainerK);
 194         fOutList->Add(fCFContainerP);
 195 }
 196
 197
 198 //________________________________________________________________
 199 void AliAnalysisTaskCentral::InitCuts(){
 200
 201 //Create and set cuts
 202
 203 //------------------EVENT LEVEL CUTS--------------------
 204 //-----------------------MC-----------------------------
 205         AliAnalysisCentralCutEvtMC *evMCCuts = new AliAnalysisCentralCutEvtMC();
 206         evMCCuts->SetMultiplicityRange(1,100);
 207         evMCCuts->SetDirectivityRange(0.0, 1.0);
 208     evMCCuts->SetDirUnitRange(0.0, 1.0);
 209
 210 //-----------------------ESD----------------------------
 211         AliAnalysisCentralCutEvtESD *evESDCuts = new AliAnalysisCentralCutEvtESD();
 212         evESDCuts->SetMultiplicityRange(1,100);
 213         evESDCuts->SetDirectivityRange(0.0, 1.0);
 214     evESDCuts->SetSPDMultiplicityRange(1,20000);
 215     evESDCuts->SetSPDDirectivityRange(0.0, 1.0);
 216
 217         TObjArray* mcEventCuts = new TObjArray();
 218         mcEventCuts->AddLast(evMCCuts);
 219
 220         TObjArray* esdEventCuts = new TObjArray();
 221         esdEventCuts->AddLast(evESDCuts);
 222
 223
 224 //------------------PARTICLE LEVEL CUTS-----------------
 225 //-------------------General MC Cuts--------------------
 226         AliAnalysisCentralCutMC *mcCutsGen = new AliAnalysisCentralCutMC();
 227         mcCutsGen->SetOnlyPrimaries(kTRUE);
 228         mcCutsGen->SetPtRange(0.2,4.0);
 229         mcCutsGen->SetEtaRange(-0.2,0.2);
 230
 231 //-------------------Specific MC Cuts-------------------
 232         AliAnalysisCentralCutMC *mcCutsPi = new AliAnalysisCentralCutMC();
 233         mcCutsPi->SetPDGCode(211); //211 pion;  321 kaon; 2212 proton
 234
 235         AliAnalysisCentralCutMC *mcCutsK = new AliAnalysisCentralCutMC();
 236         mcCutsK->SetPDGCode(321); //211 pion;  321 kaon; 2212 proton
 237
 238         AliAnalysisCentralCutMC *mcCutsP = new AliAnalysisCentralCutMC();
 239         mcCutsP->SetPDGCode(2212); //211 pion;  321 kaon; 2212 proton
 240
 241 //--------each task has its own mcList of cuts----------
 242         TObjArray *mcListGen = new TObjArray(); //general MC cuts
 243         mcListGen->AddLast(mcCutsGen);
 244
 245         TObjArray *mcListPi = new TObjArray(); //task pt pions
 246         mcListPi->AddLast(mcCutsPi);
 247
 248         TObjArray *mcListK = new TObjArray(); //task pt kaons
 249         mcListK->AddLast(mcCutsK);
 250
 251         TObjArray *mcListP = new TObjArray(); //task pt protons
 252         mcListP->AddLast(mcCutsP);
 253
 254
 255 //-------------------General ESD Cuts-------------------
 256         AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "Loose");
 257         esdCutsGen->SetMinNClustersTPC(50);
 258         esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
 259         esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
 260         esdCutsGen->SetRequireTPCRefit(kTRUE);
 261         esdCutsGen->SetAcceptKinkDaughters(kFALSE);
 262         esdCutsGen->SetMaxNsigmaToVertex(2.0);
 263         esdCutsGen->SetRequireSigmaToVertex(kTRUE);
 264         esdCutsGen->SetPtRange(0.2,4.0);
 265         esdCutsGen->SetEtaRange(-0.2,0.2);
 266
 267 //-------------------Specific ESD Cuts------------------
 268         AliAnalysisCentralCutESD *esdCutsPi = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
 269         esdCutsPi->SetPIDtype("Custom");
 270         esdCutsPi->SetPriorFunctions(kFALSE);
 271         esdCutsPi->SetPartType(kPiPlus);
 272
 273         AliAnalysisCentralCutESD *esdCutsK = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
 274         esdCutsK->SetPIDtype("Custom");
 275         esdCutsK->SetPriorFunctions(kFALSE);
 276         esdCutsK->SetPartType(kKPlus);
 277
 278         AliAnalysisCentralCutESD *esdCutsP = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
 279         esdCutsP->SetPIDtype("Custom");
 280         esdCutsP->SetPriorFunctions(kFALSE);
 281         esdCutsP->SetPartType(kProton);
 282
 283 //--------each task has its own esdList of cuts---------
 284         TObjArray* esdListGen = new TObjArray(); //general ESD track cuts
 285         esdListGen->AddLast(esdCutsGen);
 286
 287         TObjArray* esdListPi = new TObjArray();
 288         esdListPi->AddLast(esdCutsPi);
 289
 290         TObjArray* esdListK = new TObjArray();
 291         esdListK->AddLast(esdCutsK);
 292
 293         TObjArray* esdListP = new TObjArray();
 294         esdListP->AddLast(esdCutsP);
 295
 296 //------set the cuts to the RIGHT! fCutsList slots-------
 297 // // event level cuts
 298         SetCuts(0, mcEventCuts);
 299         SetCuts(1, esdEventCuts);
 300
 301  //particle level cuts
 302         SetCuts(2, mcListGen);
 303         SetCuts(3, mcListPi);
 304         SetCuts(4, mcListK);
 305         SetCuts(5, mcListP);
 306
 307         SetCuts(6, esdListGen);
 308         SetCuts(7, esdListPi);
 309         SetCuts(8, esdListK);
 310         SetCuts(9, esdListP);
 311
 312 }
 313
 314
 315 //_______________________________________________________________________
 316 void AliAnalysisTaskCentral::SendEvent(TObject *obj) const{
 317
 318 // Some cuts (ie MC IsPrimary) need the MC Event info
 319
 320     if (!fCutsList) {
 321                 printf("No particle cut list found!\n\n");
 322                 return;
 323     }
 324     else {
 325                 for(Int_t isel=0;isel< 10; isel++){
 326                         if(!fCutsList[isel]) continue;
 327                         TObjArrayIter iter(fCutsList[isel]);
 328                         AliAnalysisCuts *cut = 0;
 329
 330                         while ((cut = (AliAnalysisCuts*)iter.Next())) {
 331
 332                                 TString cutName=cut->GetName();
 333                                 if(!cutName){
 334                                         printf("No cutname!\n");
 335                                         return;
 336                                 }
 337
 338                                 Bool_t checkCut=cutName.Contains("AliAnalysisCentralCutMC");
 339
 340                                 if(checkCut){
 341                                         AliAnalysisCentralCutMC *newcut = dynamic_cast<AliAnalysisCentralCutMC *>(cut);
 342                                         newcut->ReceiveEvt(obj);
 343                                 }
 344                         }
 345                 }
 346     }
 347
 348 }
 349
 350
 351 //________________________________________________________________________
 352 Bool_t AliAnalysisTaskCentral::CheckCuts(Int_t no, TObject *obj) const{
 353
 354 // For each cut run IsSelected();
 355 //      printf("AliAnalysisTaskCentral::CheckCuts IN\n");
 356
 357     if(no > 10){
 358                 printf("\nAliAnalysisTaskCentral::CheckCuts -> Cut number is not ok! \n");
 359                 return kFALSE;
 360     }
 361
 362     if(!fCutsList[no]){
 363                 printf("AliAnalysisTaskCentral::CheckCuts -> cuts list problem! \n");
 364                 return kFALSE;
 365     }
 366
 367     TObjArrayIter iter(fCutsList[no]);
 368     AliAnalysisCuts *cut = 0;
 369
 370     while((cut = (AliAnalysisCuts*)iter.Next())){
 371                 if(!cut->IsSelected(obj)){
 372 //                      printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
 373                         return kFALSE;
 374                 }
 375     }
 376
 377 //      printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
 378     return kTRUE;
 379 }
 380
 381
 382 //________________________________________________________________________
 383 void AliAnalysisTaskCentral::Exec(Option_t *) {
 384
 385 // Main loop
 386 // Called for each event
 387
 388         Double_t pt, eta;
 389         Double_t w = 1.0;
 390         const Int_t nvar=2; //number of variables on the grid:pt,vtx
 391         Double_t value[nvar];//used to fill the CF container
 392
 393         if(fSim){  // if running on simulations -> look at MC Truth
 394
 395                 if (!fMC) {
 396                         Printf("ERROR: fMC not available");
 397                         return;
 398                 }
 399
 400                 if(CheckCuts(0, fMC)){ //check event level cuts
 401
 402                         SendEvent(fMC);
 403
 404                         // MC loop
 405                         for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) {
 406
 407                                 AliMCParticle *particle  = dynamic_cast<AliMCParticle*>(fMC->GetTrack(ipart));
 408
 409                                 if(!particle){
 410                                                 printf("\nMCParticle pointer is null!!!\n");
 411                                                 continue;
 412                                 }
 413
 414                                 if (!CheckCuts(2, particle)) continue; //check the MC general particle cuts
 415
 416                                 pt = particle->Pt();
 417                                 eta = particle->Eta();
 418
 419                                 if(pt>0) w = 1.0/pt; //invariant distribution
 420
 421                                 value[0]=pt;
 422                                 value[1]=eta;
 423
 424                                 if(CheckCuts(3, particle)){  //fill the right container for each particle
 425                                         fCFContainerPi->Fill(value,0,w);
 426                                 }
 427                                 else if(CheckCuts(4, particle)){
 428                                         fCFContainerK->Fill(value,0,w);
 429                                 }
 430                                 else if(CheckCuts(5, particle)){
 431                                         fCFContainerP->Fill(value,0,w);
 432                                 }
 433                         } //end MC particle loop
 434                 }
 435         }
 436         else{  // if we DONT run in simulated data we fill the MC step of the CFCont with 0
 437                 value[0]=0;
 438                 value[1]=0;
 439
 440                 fCFContainerPi->Fill(value,0);
 441                 fCFContainerK->Fill(value,0);
 442                 fCFContainerP->Fill(value,0);
 443         }
 444
 445
 446     if (!fESD) {
 447         Printf("ERROR: fESD not available");
 448         return;
 449     }
 450
 451     if(CheckCuts(1, fESD)){
 452
 453                 Printf("There are %d ESD tracks in this event\n", fESD->GetNumberOfTracks());
 454
 455                 // ESD loop
 456                 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 457
 458                         AliESDtrack* track = fESD->GetTrack(iTracks);
 459
 460                         if (!track) {
 461                                         Printf("ERROR: Could not receive track %d", iTracks);
 462                                 continue;
 463                         }
 464
 465                         if(!CheckCuts(6, track)) continue; //check general ESD track cuts
 466
 467                         pt = track->Pt();
 468                         eta = track->Eta();
 469
 470                         if(pt>0) w = 1.0/pt; // invariant
 471
 472                         value[0]=pt;
 473                         value[1]=eta;
 474
 475                         if(CheckCuts(7, track)){
 476                                 fCFContainerPi->Fill(value,1,w);
 477                         }
 478                         else if(CheckCuts(8, track)){
 479                                 fCFContainerK->Fill(value,1,w);
 480                         }
 481                         else if(CheckCuts(9, track)){
 482                                 fCFContainerP->Fill(value,1,w);
 483                         }
 484
 485                 } //end ESD track loop
 486
 487                 fNoEvt->Fill(0); //get the number of analyzed events
 488         }
 489
 490   // Post output data.
 491         PostData(0, fOutList);
 492
 493 }
 494
 495 //________________________________________________________________________
 496 void AliAnalysisTaskCentral::Terminate(Option_t *) {
 497
 498 // Called once at the end of the query
 499         printf("\n\n****************************************\n");
 500     printf("\tAliAnalysisCentralExtrapolate Terminate... \n");
 501
 502         TList *outList = dynamic_cast<TList*>(GetOutputData(0));
 503         if(!outList){
 504                 printf("Unable to get output list! \n");
 505                 return;
 506         }
 507
 508         AliAnalysisCentralExtrapolate *extPi = new AliAnalysisCentralExtrapolate("extrapolationpi");
 509         extPi->SetInputList(outList);  //send the outlist to the extrapolation class
 510         extPi->SetParticle("kPiPlus"); //set the particle type
 511         extPi->ApplyEff();             //correct the pt distribution !!HAS TO RUN BEFORE extrapolation!!
 512         extPi->BoltzmannFit();         //fit and extrapolate using Boltzmann-Gibbs Blast wave model
 513         extPi->TsallisFit();           //fit and extrapolate using Tsallis Blast wave model
 514         TList *extOutListPi = extPi->GetOutputList();
 515
 516         AliAnalysisCentralExtrapolate *extK = new AliAnalysisCentralExtrapolate("extrapolationK");
 517         extK->SetInputList(outList);
 518         extK->SetParticle("kKPlus");
 519         extK->ApplyEff();
 520         extK->BoltzmannFit();
 521         extK->TsallisFit();
 522         TList *extOutListK = extK->GetOutputList();
 523
 524         AliAnalysisCentralExtrapolate *extP = new AliAnalysisCentralExtrapolate("extrapolationP");
 525         extP->SetInputList(outList);
 526         extP->SetParticle("kProton");
 527         extP->ApplyEff();
 528         extP->BoltzmannFit();
 529         extP->TsallisFit();
 530         TList *extOutListP = extP->GetOutputList();
 531
 532
 533     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 534         if(!mgr){
 535                 printf("Unable to get AnalysisManager! \n");
 536         return;
 537         }
 538
 539         TString anaType;
 540         mgr->GetAnalysisTypeString(anaType);
 541
 542         if(anaType.Contains("local")){
 543                 fOutList->Add(extOutListPi);
 544                 fOutList->Add(extOutListK);
 545                 fOutList->Add(extOutListP);
 546         }
 547         else{
 548                 AliAnalysisDataContainer *cont = GetOutputSlot(0)->GetContainer();
 549                 if(!cont){
 550                         printf("Unable to get DataContainer! \n");
 551                         return;
 552                 }
 553
 554                 printf("file name = %s\n", cont->GetFileName());
 555                 TFile file(cont->GetFileName(),"update");
 556
 557                 file.cd("PWG2Central");
 558
 559                 gFile->WriteObject(extOutListPi,"pion_list","SingleKey");
 560                 gFile->WriteObject(extOutListK,"kaon_list","SingleKey");
 561                 gFile->WriteObject(extOutListP,"proton_list","SingleKey");
 562                 file.Close();
 563         }
 564
 565
 566 }