PWG4/macros/AddTaskPi0EMCALPbPb.C

   1
   2 AliAnalysisTaskParticleCorrelation *AddTaskPi0EMCALPbPb(
   3                                                     const TString data          = "AOD",
   4                                                     const TString calorimeter   = "EMCAL",
   5                                                     const Bool_t  printSettings = kFALSE,
   6                                                     const Bool_t  simulation    = kFALSE,
   7                                                     const Bool_t  outputAOD     = kFALSE,
   8                                                     const TString outputfile    = "",
   9                                                     const Int_t   year          = 2011,
  10                                                     const TString col           = "PbPb",
  11                                                     const TString trig          = "",
  12                                                     const TString clustersArray = ""
  13                                                     )
  14 {
  15   // Creates a PartCorr task, configures it and adds it to the analysis manager.
  16
  17   // Get the pointer to the existing analysis manager via the static access method.
  18
  19   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  20   if (!mgr) {
  21     ::Error("AddTask", "No analysis manager to connect to.");
  22     return NULL;
  23   }
  24
  25   // Check the analysis type using the event handlers connected to the analysis manager.
  26
  27   if (!mgr->GetInputEventHandler()) {
  28     ::Error("AddTask", "This task requires an input event handler");
  29     return NULL;
  30   }
  31   Bool_t kInputDataType = "AOD";
  32   if(!data.Contains("delta"))
  33     kInputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
  34
  35   Bool_t kUseKinematics = kFALSE;
  36   if(simulation) {
  37     kUseKinematics = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
  38     if (!kUseKinematics && data=="AOD" && kInputDataType != "ESD") kUseKinematics = kTRUE; //AOD primary should be available ...
  39   }
  40
  41   //cout<<"********* ACCESS KINE? "<<kUseKinematics<<endl;
  42
  43   // #### Configure analysis ####
  44
  45   AliAnaPartCorrMaker * maker = new AliAnaPartCorrMaker();
  46
  47   // General frame setting and configuration
  48   maker->SetReader   (ConfigureReader(data,kInputDataType,calorimeter,clustersArray,col,outputAOD,kUseKinematics,printSettings)   );
  49   maker->SetCaloUtils(ConfigureCaloUtils(year,col,clustersArray, simulation, printSettings));
  50
  51   // Analysis tasks setting and configuration
  52   Int_t n = 0;//Analysis number, order is important
  53   maker->AddAnalysis(ConfigurePhotonAnalysis(data,calorimeter,year,col,clustersArray,simulation,"",printSettings), n++); // Photon cluster selection
  54   maker->AddAnalysis(ConfigurePi0Analysis   (data,calorimeter,year,col,clustersArray,simulation,"",printSettings), n++); // Pi0 invariant mass analysis
  55
  56   maker->SetAnaDebug(-1)  ;
  57   maker->SwitchOnHistogramsMaker()  ;
  58   if(data.Contains("delta")) maker->SwitchOffAODsMaker() ;
  59   else                       maker->SwitchOnAODsMaker()  ;
  60
  61   if(printSettings) maker->Print("");
  62
  63   //printf("<< End Configuration of %d analysis for calorimeter %s >>\n",n, calorimeter.Data());
  64
  65   // Create task
  66
  67   AliAnalysisTaskParticleCorrelation * task = new AliAnalysisTaskParticleCorrelation (Form("PartCorr%s_Trig%s_Cl%s",calorimeter.Data(),trig.Data(),clustersArray.Data()));
  68   task->SetConfigFileName(""); //Don't configure the analysis via configuration file.
  69   //task->SetDebugLevel(-1);
  70   task->SetBranches("ESD:AliESDRun.,AliESDHeader"); //just a trick to get Constantin's analysis to work
  71   task->SetAnalysisMaker(maker);
  72   mgr->AddTask(task);
  73
  74   //Create containers
  75
  76   if(outputfile.Length()==0)outputfile = AliAnalysisManager::GetCommonFileName();
  77   TString name(Form("%s_Trig%s_Cl%s",calorimeter.Data(),trig.Data(),clustersArray.Data()));
  78   AliAnalysisDataContainer *cout_pc   = mgr->CreateContainer(name.Data(), TList::Class(),
  79                                                              AliAnalysisManager::kOutputContainer,
  80                                                              Form("%s:%s",outputfile.Data(),name.Data()));
  81
  82   AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer(Form("%sCuts",name.Data()), TList::Class(),
  83                                                              AliAnalysisManager::kParamContainer,
  84                                                              Form("%s:%sCuts",outputfile.Data(),name.Data()));
  85
  86   // Create ONLY the output containers for the data produced by the task.
  87   // Get and connect other common input/output containers via the manager as below
  88   //==============================================================================
  89   mgr->ConnectInput  (task, 0, mgr->GetCommonInputContainer());
  90   // AOD output slot will be used in a different way in future
  91   if(!data.Contains("delta")   && outputAOD) mgr->ConnectOutput (task, 0, mgr->GetCommonOutputContainer());
  92   mgr->ConnectOutput (task, 1, cout_pc);
  93   mgr->ConnectOutput (task, 2, cout_cuts);
  94
  95   return task;
  96 }
  97
  98 //____________________________________
  99 AliCaloTrackReader * ConfigureReader(TString kData,TString kInputDataType, TString kCalorimeter,
 100                                      TString kClusterArray, TString kCollisions,
 101                                      Bool_t kOutputAOD, Bool_t kUseKinematics, Bool_t kPrint)
 102 {
 103
 104   AliCaloTrackReader * reader = 0;
 105   if     (kData.Contains("AOD"))   reader = new AliCaloTrackAODReader();
 106   else if(kData=="ESD")            reader = new AliCaloTrackESDReader();
 107   else if(kData=="MC" &&
 108           kInputDataType == "ESD") reader = new AliCaloTrackMCReader();
 109
 110   reader->SetDebug(-1);//10 for lots of messages
 111
 112   //Delta AOD?
 113   //reader->SetDeltaAODFileName("");
 114   if(kOutputAOD) reader->SwitchOnWriteDeltaAOD()  ;
 115
 116   // MC settings
 117   if(kUseKinematics){
 118     if(kInputDataType == "ESD"){
 119       reader->SwitchOnStack();
 120       reader->SwitchOffAODMCParticles();
 121     }
 122     else if(kInputDataType == "AOD"){
 123       reader->SwitchOffStack();
 124       reader->SwitchOnAODMCParticles();
 125     }
 126   }
 127
 128   //------------------------
 129   // Detector input filling
 130   //------------------------
 131
 132   //Min cluster/track E
 133   reader->SetEMCALEMin(0.5);
 134   reader->SetEMCALEMax(1000);
 135   reader->SetPHOSEMin(0.3);
 136   reader->SetPHOSEMax(1000);
 137   reader->SetCTSPtMin(0.1);
 138   reader->SetCTSPtMax(1000);
 139
 140   reader->SwitchOffFiducialCut();
 141
 142   // Tracks
 143   reader->SwitchOffCTS();
 144   //gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/CreateTrackCutsPWG4.C");
 145   //AliESDtrackCuts * esdTrackCuts = CreateTrackCutsPWG4(10011004);
 146   //reader->SetTrackCuts(esdTrackCuts);
 147
 148   // Calorimeter
 149
 150   reader->SetEMCALClusterListName(kClusterArray);
 151   if(kClusterArray == "") {
 152     //printf("**************** Normal analysis **************** \n");
 153     reader->SwitchOnClusterRecalculation(); // Bad map removal
 154   }
 155   else {
 156     printf("**************** Input for analysis is Clusterizer %s **************** \n", kClusterArray.Data());
 157     reader->SwitchOffClusterRecalculation();
 158   }
 159
 160   if(kCalorimeter == "EMCAL") {
 161     reader->SwitchOnEMCALCells();
 162     reader->SwitchOnEMCAL();
 163   }
 164   if(kCalorimeter == "PHOS") {
 165     reader->SwitchOnPHOSCells();
 166     reader->SwitchOnPHOS();
 167   }
 168
 169   // for case data="deltaAOD", no need to fill the EMCAL/PHOS cluster lists
 170   if(kData.Contains("delta")){
 171     reader->SwitchOffEMCAL();
 172     reader->SwitchOffPHOS();
 173     reader->SwitchOffEMCALCells();
 174     reader->SwitchOffPHOSCells();
 175   }
 176
 177   //-----------------
 178   // Event selection
 179   //-----------------
 180
 181   //if(!kUseKinematics) reader->SetFiredTriggerClassName("CEMC7EGA-B-NOPF-CENTNOTRD"); // L1 Gamma
 182
 183   if     (kCollisions=="pp"  ) {
 184     reader->SwitchOffEventSelection();
 185     reader->SwitchOffV0ANDSelection() ;        // and besides v0 AND
 186     reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
 187     reader->SetZvertexCut(50.);                // Open cut
 188   }
 189   else if(kCollisions=="PbPb") {
 190     reader->SwitchOffEventSelection();         // remove pileup by default
 191     reader->SwitchOffV0ANDSelection() ;        // and besides v0 AND
 192     reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
 193     reader->SetZvertexCut(10.);                // Centrality defined in this range.
 194
 195     // Centrality
 196     reader->SetCentralityClass("V0M");
 197     reader->SetCentralityOpt(10);  // 10 centrality bins
 198     reader->SetCentralityBin(-1,-1); // Accept all events, if not select range
 199
 200     // Event plane (only used in AliAnaPi0 for the moment)
 201     reader->SetEventPlaneMethod("Q");
 202   }
 203
 204   if(kPrint) reader->Print("");
 205
 206   return reader;
 207
 208 }
 209
 210 //_______________________________________
 211 AliCalorimeterUtils* ConfigureCaloUtils(Int_t kYears, TString kCollisions, TString kClusterArray,
 212                                         Bool_t kSimulation, Bool_t kPrint)
 213 {
 214
 215   AliCalorimeterUtils *cu = new AliCalorimeterUtils;
 216   cu->SetDebug(-1);
 217
 218   if(kYears==2010) cu->SetEMCALGeometryName("EMCAL_FIRSTYEARV1");
 219   else             cu->SetEMCALGeometryName("EMCAL_COMPLETEV1");
 220
 221   // Remove clusters close to borders, at least max energy cell is 1 cell away
 222   cu->SetNumberOfCellsFromEMCALBorder(1);
 223   cu->SetNumberOfCellsFromPHOSBorder(2);
 224
 225   //if(kClusterArray == "")
 226   //  cu->SwitchOffRecalculateClusterTrackMatching(); // Done in clusterization
 227   //else
 228   //  cu->SwitchOnRecalculateClusterTrackMatching();
 229
 230   //EMCAL only settings
 231   AliEMCALRecoUtils * recou = cu->GetEMCALRecoUtils();
 232
 233
 234   if(kCollisions=="pp"  ) { // Do only for pp for the moment
 235     cu->SwitchOnCorrectClusterLinearity();
 236     if(!kSimulation) recou->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected);
 237     else             recou->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
 238   }
 239
 240   recou->SwitchOnRejectExoticCell();
 241   if(kClusterArray == "") recou->SwitchOnRejectExoticCluster();
 242   else                    recou->SwitchOffRejectExoticCluster();
 243
 244   if(kPrint) cu->Print("");
 245
 246   return cu;
 247
 248 }
 249
 250 //_____________________________________
 251 AliAnaPhoton* ConfigurePhotonAnalysis(TString kData, TString kCalorimeter, Int_t kYears, TString kCollisions, TString kClusterArray,
 252                                       Bool_t kSimulation, TString kTrig = "", Bool_t kPrint)
 253 {
 254
 255   AliAnaPhoton *anaphoton = new AliAnaPhoton();
 256   anaphoton->SetDebug(-1); //10 for lots of messages
 257
 258   // cluster selection cuts
 259
 260   anaphoton->SwitchOffFiducialCut();
 261
 262   anaphoton->SetCalorimeter(kCalorimeter);
 263
 264   if(kCalorimeter == "PHOS"){
 265     anaphoton->SetNCellCut(2);// At least 2 cells
 266     anaphoton->SetMinPt(0.3);
 267     anaphoton->SetMinDistanceToBadChannel(2, 4, 5);
 268   }
 269   else {//EMCAL
 270     anaphoton->SetNCellCut(1);// At least 2 cells
 271     anaphoton->SetMinPt(0.5); // avoid mip peak at E = 260 MeV
 272     anaphoton->SetMaxPt(1000);
 273     anaphoton->SetTimeCut(-1000,1000); // open cut, usual time window of [425-825] ns if time recalibration is off
 274     anaphoton->SetMinDistanceToBadChannel(1, 2, 3); // For filtered AODs, new releases.
 275   }
 276
 277   anaphoton->SwitchOnTrackMatchRejection() ;
 278
 279   //PID cuts (shower shape)
 280   anaphoton->SwitchOnCaloPID(); // do PID selection, unless specified in GetCaloPID, selection not based on bayesian
 281   AliCaloPID* caloPID = anaphoton->GetCaloPID();
 282   //Not used in bayesian
 283
 284   //EMCAL
 285   caloPID->SetEMCALLambda0CutMax(0.30);
 286   caloPID->SetEMCALLambda0CutMin(0.10);
 287
 288   caloPID->SetEMCALDEtaCut(0.025);
 289   caloPID->SetEMCALDPhiCut(0.05);
 290   // In case of official AODs when dX and dZ was not stored, open the cuts
 291   // and rely on having a match recorded. In case of reclusterization, try.
 292   if(kData=="AOD" && kClusterArray==""){
 293     caloPID->SetEMCALDEtaCut(2000);
 294     caloPID->SetEMCALDPhiCut(2000);
 295   }
 296
 297   //PHOS
 298   caloPID->SetPHOSDispersionCut(2.5);
 299   caloPID->SetPHOSRCut(2.);
 300   if(kData=="AOD") caloPID->SetPHOSRCut(2000.); // Open cut since dX, dZ not stored
 301
 302   if(kCalorimeter=="PHOS"){
 303     caloPID->SetHistoDEtaRangeAndNBins(-200, 200, 200); // dZ
 304     caloPID->SetHistoDPhiRangeAndNBins(-200, 200, 200); // dX
 305   }
 306
 307   //caloPID->SetTOFCut(10000000); // Not used, only to set PID bits
 308
 309   anaphoton->SwitchOffFillShowerShapeHistograms();  // Filled before photon shower shape selection
 310
 311   // Input / output delta AOD settings
 312
 313   if(!kData.Contains("delta")) {
 314     anaphoton->SetOutputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
 315     anaphoton->SetOutputAODClassName("AliAODPWG4ParticleCorrelation");
 316     //anaphoton->SetOutputAODClassName("AliAODPWG4Particle"); // use if no correlation done
 317   }
 318   else anaphoton->SetInputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
 319
 320   //Set Histograms name tag, bins and ranges
 321
 322   anaphoton->AddToHistogramsName("AnaPhoton_");
 323   SetHistoRangeAndNBins(anaphoton); // see method below
 324
 325   // Number of particle type MC histograms
 326   anaphoton->FillNOriginHistograms(8);
 327   anaphoton->FillNPrimaryHistograms(4);
 328
 329   if(kPrint) anaphoton->Print("");
 330
 331   return anaphoton;
 332
 333 }
 334
 335
 336
 337 //_______________________________
 338 AliAnaPi0* ConfigurePi0Analysis(TString kData, TString kCalorimeter, Int_t kYears, TString kCollisions, TString kClusterArray,
 339                                 Bool_t kSimulation, TString kTrig = "", Bool_t kPrint)
 340 {
 341
 342   AliAnaPi0 *anapi0 = new AliAnaPi0();
 343
 344   anapi0->SetDebug(-1);//10 for lots of messages
 345   if(kPrint) anapi0->Print("");
 346
 347   // Input delta AOD settings
 348   anapi0->SetInputAODName(Form("Photon%s_Trig%s_Cl%s",kCalorimeter.Data(), kTrig.Data(),kClusterArray.Data()));
 349
 350   // Calorimeter settings
 351   anapi0->SetCalorimeter(kCalorimeter);
 352   if(kCalorimeter=="PHOS") anapi0->SetNumberOfModules(3); //PHOS first year
 353   else {
 354     if(kYears==2010) anapi0->SetNumberOfModules(4); //EMCAL first year
 355     else             anapi0->SetNumberOfModules(10);
 356   }
 357
 358   //settings for pp collision mixing
 359   anapi0->SwitchOnOwnMix(); //Off when mixing done with general mixing frame
 360
 361   // Cuts
 362   if(kCalorimeter=="EMCAL") anapi0->SetPairTimeCut(70);
 363
 364   if     (kCollisions=="pp"  ) {
 365     printf("****** Configure Pi0 for pp analysis\n");
 366     anapi0->SetNCentrBin(1);
 367     anapi0->SetNZvertBin(10);
 368     anapi0->SetNRPBin(1);
 369     anapi0->SetNMaxEvMix(100);
 370   }
 371   else if(kCollisions=="PbPb") {
 372     printf("****** Configure Pi0 for PbPb analysis\n");
 373     anapi0->SetNCentrBin(5);
 374     anapi0->SetNZvertBin(3);
 375     anapi0->SetNRPBin(1);
 376     anapi0->SetNMaxEvMix(5);
 377   }
 378
 379   anapi0->SwitchOffMultipleCutAnalysis();
 380   anapi0->SwitchOffSMCombinations();
 381
 382   //Set Histograms name tag, bins and ranges
 383
 384   anapi0->AddToHistogramsName("AnaPi0_");
 385   SetHistoRangeAndNBins(anapi0, kCalorimeter, kYears, kCollisions, kSimulation); // see method below
 386
 387   return anapi0;
 388
 389 }
 390
 391 //________________________________________________________
 392 void SetHistoRangeAndNBins (AliAnaPartCorrBaseClass* ana, TString kCalorimeter,
 393                             Int_t kYears, TString kCollisions, Bool_t kSimulation)
 394 {
 395   // Set common bins for all analysis and MC histograms filling
 396
 397   if(kSimulation) ana->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet.
 398   else            ana->SwitchOffDataMC() ;
 399
 400   ana->SetHistoPtRangeAndNBins(0, 100, 250) ; // Energy and pt histograms
 401
 402   if(kCalorimeter=="EMCAL"){
 403     if(kYears==2010){
 404       ana->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 122*TMath::DegToRad(), 78) ;
 405       ana->SetHistoXRangeAndNBins(-230,90,120); // QA
 406       ana->SetHistoYRangeAndNBins(370,450,40);  // QA
 407     }
 408     else {
 409       ana->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 182*TMath::DegToRad(), 108) ;
 410       ana->SetHistoXRangeAndNBins(-600,90,200); // QA
 411       ana->SetHistoYRangeAndNBins(100,450,100); // QA
 412     }
 413
 414     ana->SetHistoEtaRangeAndNBins(-0.72, 0.72, 144) ;
 415   }
 416   else{
 417     ana->SetHistoPhiRangeAndNBins(260*TMath::DegToRad(), 320*TMath::DegToRad(), 60) ;
 418     ana->SetHistoEtaRangeAndNBins(-0.13, 0.13, 130) ;
 419
 420   }
 421
 422   ana->SetHistoShowerShapeRangeAndNBins(0, 3, 300);
 423
 424   // Invariant mass analysis
 425   ana->SetHistoMassRangeAndNBins(0., 1., 200) ;
 426   ana->SetHistoAsymmetryRangeAndNBins(0., 1. , 100) ;
 427
 428   // check if time calibration is on
 429   ana->SetHistoTimeRangeAndNBins(-1000.,1000,1000);
 430   ana->SetHistoDiffTimeRangeAndNBins(-200, 200, 800);
 431
 432   // QA, electron, charged
 433   ana->SetHistoPOverERangeAndNBins(0,10.,100);
 434   ana->SetHistodEdxRangeAndNBins(0.,200.,200);
 435
 436   // QA
 437   ana->SetHistoFinePtRangeAndNBins(0, 10, 200) ; // bining for fhAmpId
 438   ana->SetHistodRRangeAndNBins(0.,TMath::Pi(),150);
 439   ana->SetHistoRatioRangeAndNBins(0.,2.,100);
 440   ana->SetHistoVertexDistRangeAndNBins(0.,500.,500);
 441   ana->SetHistoNClusterCellRangeAndNBins(0,500,500);
 442   ana->SetHistoZRangeAndNBins(-400,400,200);
 443   ana->SetHistoRRangeAndNBins(400,450,25);
 444   ana->SetHistoV0SignalRangeAndNBins(0,5000,500);
 445   ana->SetHistoV0MultiplicityRangeAndNBins(0,5000,500);
 446   ana->SetHistoTrackMultiplicityRangeAndNBins(0,5000,500);
 447
 448 }
 449