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