PWG1/ITS/AliAnalysisTaskSPD.cxx

   1 #include "AliAnalysisTask.h"\r
   2 #include "AliAnalysisManager.h"\r
   3 #include "AliAnalysisDataContainer.h"\r
   4 #include "AliMultiplicity.h"\r
   5 #include "AliITSRecPoint.h"\r
   6 #include "AliITSDetTypeRec.h"\r
   7 #include "AliGeomManager.h"\r
   8 #include "AliESDVertex.h"\r
   9 #include "AliESDEvent.h"\r
  10 #include "AliMagF.h"\r
  11 #include <TSystem.h>\r
  12 #include <TTree.h>\r
  13 #include <TNtuple.h>\r
  14 #include <TH1I.h>\r
  15 #include <TH2F.h>\r
  16 #include <TGeoManager.h>\r
  17 #include <TChain.h>\r
  18 #include <TGeoGlobalMagField.h>\r
  19 #include "AliESDInputHandlerRP.h"\r
  20 #include "AliAnalysisTaskSPD.h"\r
  21 #include "AliITSsegmentationSPD.h"\r
  22 \r
  23 ClassImp(AliAnalysisTaskSPD)\r
  24   //______________________________________________________________________________\r
  25   const Int_t AliAnalysisTaskSPD::fgkDDLModuleMap[20][12] = {\r
  26     { 4, 5, 0, 1, 80, 81, 84, 85, 88, 89, 92, 93},\r
  27     {12,13, 8, 9, 96, 97,100,101,104,105,108,109},\r
  28     {20,21,16,17,112,113,116,117,120,121,124,125},\r
  29     {28,29,24,25,128,129,132,133,136,137,140,141},\r
  30     {36,37,32,33,144,145,148,149,152,153,156,157},\r
  31     {44,45,40,41,160,161,164,165,168,169,172,173},\r
  32     {52,53,48,49,176,177,180,181,184,185,188,189},\r
  33     {60,61,56,57,192,193,196,197,200,201,204,205},\r
  34     {68,69,64,65,208,209,212,213,216,217,220,221},\r
  35     {76,77,72,73,224,225,228,229,232,233,236,237},\r
  36     { 7, 6, 3, 2, 83, 82, 87, 86, 91, 90, 95, 94},\r
  37     {15,14,11,10, 99, 98,103,102,107,106,111,110},\r
  38     {23,22,19,18,115,114,119,118,123,122,127,126},\r
  39     {31,30,27,26,131,130,135,134,139,138,143,142},\r
  40     {39,38,35,34,147,146,151,150,155,154,159,158},\r
  41     {47,46,43,42,163,162,167,166,171,170,175,174},\r
  42     {55,54,51,50,179,178,183,182,187,186,191,190},\r
  43     {63,62,59,58,195,194,199,198,203,202,207,206},\r
  44     {71,70,67,66,211,210,215,214,219,218,223,222},\r
  45     {79,78,75,74,227,226,231,230,235,234,239,238}\r
  46   };\r
  47 \r
  48 //______________________________________________________________________________\r
  49 AliAnalysisTaskSPD::AliAnalysisTaskSPD() : \r
  50   AliAnalysisTaskSE(), \r
  51   fOutput(0),\r
  52   fSegSPD(0)\r
  53 {\r
  54   //\r
  55 }\r
  56 //______________________________________________________________________________\r
  57 AliAnalysisTaskSPD::AliAnalysisTaskSPD(const char *name) : \r
  58   AliAnalysisTaskSE(name), \r
  59   fOutput(0),\r
  60   fSegSPD(0)\r
  61 {\r
  62   DefineOutput(1, TList::Class());\r
  63 }\r
  64 \r
  65 \r
  66 //___________________________________________________________________________\r
  67 AliAnalysisTaskSPD::~AliAnalysisTaskSPD(){\r
  68   //\r
  69   if (fOutput) {\r
  70     delete fOutput;\r
  71     fOutput = 0;\r
  72   }\r
  73   if(fSegSPD) delete fSegSPD;\r
  74 }\r
  75 \r
  76 //________________________________________________________________________\r
  77 void AliAnalysisTaskSPD::Init()\r
  78 {\r
  79   // Initialization\r
  80 \r
  81   if(fDebug > 1) printf("AliAnalysisTaskSPD::Init() \n");\r
  82 }\r
  83 //___________________________________________________________________________\r
  84 \r
  85 void AliAnalysisTaskSPD::UserCreateOutputObjects() {\r
  86   //\r
  87   if(fDebug > 1) printf("AliAnalysisTaskSPD::UserCreateOutputObjects() \n");\r
  88 \r
  89   fSegSPD = new AliITSsegmentationSPD();\r
  90 \r
  91   OpenFile(1);\r
  92   fOutput = new TList();\r
  93 \r
  94   TH1I *nentries = new TH1I("events","events",1,0,1);\r
  95   fOutput->AddLast(nentries);\r
  96 \r
  97   TH2F *modulemap[240];\r
  98 \r
  99   for(Int_t imod =0; imod < 240; imod++){\r
 100     modulemap[imod] = new TH2F(Form("mod%i",imod),Form("cluster map for module %i",imod),800,-4,4,200,-1,1);\r
 101     modulemap[imod]->SetXTitle("Local Z (cm)");\r
 102     modulemap[imod]->SetYTitle("Local X (cm)");\r
 103     fOutput->AddLast(modulemap[imod]);\r
 104   }\r
 105 \r
 106   TH1F *hFOgood = new TH1F("hFOgood"," Events with FO and clusters ",1200,0,1200);\r
 107   hFOgood->SetXTitle("chipkey");\r
 108   fOutput->AddLast(hFOgood);\r
 109 \r
 110   TH1F *hFOnoisy = new TH1F("hFOnoisy"," Events with FO but no cluster",1200,0,1200);\r
 111   hFOnoisy->SetXTitle("chipkey");\r
 112   fOutput->AddLast(hFOnoisy);\r
 113   \r
 114   TH1F *hFiredChips = new TH1F("hFiredChips"," yield of fired chips",1200,0,1200);\r
 115   hFiredChips->SetXTitle("chipkey");\r
 116   fOutput->AddLast(hFiredChips);\r
 117 \r
 118   TH2F *hSPDphivsSPDeta= new TH2F("hSPDphivsSPDeta", "Tracklets - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());\r
 119   hSPDphivsSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");\r
 120   hSPDphivsSPDeta->GetYaxis()->SetTitle("#varphi [rad]");\r
 121   fOutput->AddLast(hSPDphivsSPDeta);\r
 122 \r
 123   TH1F *hSPDphiZpos= new TH1F("hSPDphiZpos", "Tracklets - #varphi (Z>0)",360,0.,2*TMath::Pi());\r
 124   hSPDphiZpos->SetXTitle("#varphi [rad]");\r
 125   fOutput->AddLast(hSPDphiZpos);\r
 126 \r
 127   TH1F *hSPDphiZneg= new TH1F("fHistSPDphivsSPDetaZneg", "Tracklets - #varphi (Z<0)",360,0.,2*TMath::Pi());\r
 128   hSPDphiZneg->SetXTitle("#varphi [rad]");\r
 129   fOutput->AddLast(hSPDphiZneg);\r
 130 \r
 131 }\r
 132 //______________________________________________________________________________\r
 133 void AliAnalysisTaskSPD::UserExec(Option_t */*option*/)\r
 134 {\r
 135   if(fDebug > 1) printf("AliAnalysisTaskSPD::Exec() \n");\r
 136   \r
 137 \r
 138   AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
 139   if(!hand) {\r
 140     printf("No AliESDInputHandlerRP \n");\r
 141     return;\r
 142   }  \r
 143 \r
 144   AliESDEvent *ESD = hand->GetEvent();\r
 145   if(!ESD) {\r
 146     printf("No AliESDEvent \n");\r
 147     return;\r
 148   }  \r
 149 \r
 150   TTree * treeRP = hand->GetTreeR("ITS");\r
 151   if(!treeRP) {\r
 152     printf("No RecPoints tree \n");\r
 153     return;\r
 154   }  \r
 155 \r
 156   // RecPoints info\r
 157 \r
 158   TClonesArray statITSrec("AliITSRecPoint");\r
 159   TClonesArray *ITSCluster = &statITSrec;\r
 160   TBranch* branch=treeRP->GetBranch("ITSRecPoints");\r
 161   if(!branch) return;\r
 162 \r
 163   branch->SetAddress(&ITSCluster);\r
 164 \r
 165   Bool_t isfiredchip[20][6][10];\r
 166   for(Int_t eq=0; eq<20; eq++){\r
 167     for(Int_t hs =0; hs<6; hs++){\r
 168       for(Int_t chip=0; chip<10; chip++){\r
 169         isfiredchip[eq][hs][chip]=kFALSE;\r
 170       }\r
 171     }\r
 172   }\r
 173 \r
 174 \r
 175   for(Int_t iMod=0;iMod<240;iMod++){\r
 176     branch->GetEvent(iMod);\r
 177     Int_t nrecp = statITSrec.GetEntries();\r
 178     for(Int_t irec=0;irec<nrecp;irec++) {\r
 179       AliITSRecPoint *recp = (AliITSRecPoint*)statITSrec.At(irec);\r
 180       Int_t lay=recp->GetLayer();\r
 181       if(lay>1) continue;\r
 182       Float_t local[3]={-1,-1};\r
 183       local[1]=recp->GetDetLocalX();\r
 184       local[0]=recp->GetDetLocalZ();\r
 185       //printf("local X %f   local Z %f   in module %i \n",local[0],local[1],iMod);\r
 186       ((TH2F*)fOutput->At(iMod+1))->Fill(local[0],local[1]);\r
 187       Int_t eq = GetOnlineEqIdFromOffline(iMod);\r
 188       Int_t hs = GetOnlineHSFromOffline(iMod);\r
 189       Int_t row, col;\r
 190       fSegSPD->LocalToDet(0.5,local[0],row,col);\r
 191       Int_t chip = GetOnlineChipFromOffline(iMod,col);\r
 192       isfiredchip[eq][hs][chip]=kTRUE;\r
 193     }\r
 194   }\r
 195 \r
 196   //  ESDs info\r
 197 \r
 198   // First looking at the FO bits\r
 199   const AliMultiplicity *mult = ESD->GetMultiplicity();\r
 200 \r
 201   for(Int_t eq=0; eq<20; eq++){\r
 202     for(Int_t hs =0; hs<6; hs++){\r
 203       for(Int_t chip=0; chip<10; chip++){\r
 204         Int_t key = GetOfflineChipKeyFromOnline(eq,hs,chip);\r
 205         if(mult->TestFastOrFiredChips(key) && isfiredchip[eq][hs][chip]) ((TH1F*)fOutput->At(241))->Fill(key);\r
 206         if(!isfiredchip[eq][hs][chip] && mult->TestFastOrFiredChips(key)) ((TH1F*)fOutput->At(242))->Fill(key);\r
 207         if(isfiredchip[eq][hs][chip]) ((TH1F*)fOutput->At(243))->Fill(key);\r
 208       }\r
 209     }\r
 210   }\r
 211 \r
 212   // Then looking at tracklets \r
 213   const AliESDVertex* vtxESD = ESD->GetVertex();\r
 214   if(vtxESD){\r
 215 \r
 216     Double_t esdvtx[3];\r
 217     vtxESD->GetXYZ(esdvtx);\r
 218 \r
 219     for(Int_t iTracklet =0; iTracklet < mult->GetNumberOfTracklets(); iTracklet++){\r
 220 \r
 221       Float_t phiTr= mult->GetPhi(iTracklet);\r
 222       Float_t etaTr =mult->GetEta(iTracklet);\r
 223 \r
 224       ((TH2F*)fOutput->At(244))->Fill(etaTr,phiTr);\r
 225 \r
 226       //layer 0\r
 227       //Float_t z = esdvtx[2] + 3.9 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- etaTr)));\r
 228       Int_t isZpos = (esdvtx[2]>0 && etaTr>0) ;\r
 229       Int_t isZneg = (esdvtx[2]<0 && etaTr<0) ; \r
 230 \r
 231       if(isZpos)  ((TH1F*)fOutput->At(245))->Fill(phiTr);\r
 232       if(isZneg)  ((TH1F*)fOutput->At(246))->Fill(phiTr);\r
 233     }\r
 234   }\r
 235 \r
 236   ((TH1I*)fOutput->At(0))->Fill(0); // monitoring plot\r
 237 \r
 238   PostData(1,fOutput);  \r
 239   return;\r
 240   \r
 241 }\r
 242 //______________________________________________________________________________\r
 243 void AliAnalysisTaskSPD::Terminate(Option_t */*option*/)\r
 244 {\r
 245   // Terminate analysis\r
 246   if(fDebug > 1) printf("AliAnalysisTaskSPD:: Terminate() \n");\r
 247   fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
 248   if (!fOutput) {     \r
 249     printf("ERROR: fOutput not available\n");\r
 250     //return;\r
 251   }\r
 252 }\r
 253 //________________________________________________________________________\r
 254 UInt_t AliAnalysisTaskSPD::GetOfflineModuleFromOnline(UInt_t eqId, UInt_t hs, UInt_t chip) {\r
 255   // online->offline (module)\r
 256   if (eqId<20 && hs<6 && chip<10) return fgkDDLModuleMap[eqId][hs*2+chip/5];\r
 257   else return 240;\r
 258 }\r
 259 //________________________________________________________________________\r
 260 UInt_t AliAnalysisTaskSPD::GetOfflineChipKeyFromOnline(UInt_t eqId, UInt_t hs, UInt_t chip) {\r
 261   // online->offline (chip key: 0-1199)\r
 262   if (eqId<20 && hs<6 && chip<10) {\r
 263     UInt_t module = GetOfflineModuleFromOnline(eqId,hs,chip);\r
 264     UInt_t chipInModule = ( chip>4 ? chip-5 : chip );\r
 265     if(eqId>9) chipInModule = 4 - chipInModule;  // side C only\r
 266     return (module*5 + chipInModule);\r
 267   } else return 1200;\r
 268 }\r
 269 //__________________________________________________________________________\r
 270 UInt_t AliAnalysisTaskSPD::GetOnlineEqIdFromOffline(UInt_t module) {\r
 271   // offline->online (eq)\r
 272   for (UInt_t eqId=0; eqId<20; eqId++) {\r
 273     for (UInt_t iModule=0; iModule<12; iModule++) {\r
 274       if (GetModuleNumber(eqId,iModule)==(Int_t)module) return eqId;\r
 275     }\r
 276   }\r
 277   return 20; // error\r
 278 }\r
 279 //__________________________________________________________________________\r
 280 UInt_t AliAnalysisTaskSPD::GetOnlineHSFromOffline(UInt_t module) {\r
 281   // offline->online (hs)\r
 282   for (UInt_t eqId=0; eqId<20; eqId++) {\r
 283     for (UInt_t iModule=0; iModule<12; iModule++) {\r
 284       if (GetModuleNumber(eqId,iModule)==(Int_t)module) return iModule/2;\r
 285     }\r
 286   }\r
 287   return 6; // error\r
 288 }\r
 289 //__________________________________________________________________________\r
 290 UInt_t AliAnalysisTaskSPD::GetOnlineChipFromOffline(UInt_t module, UInt_t colM) {\r
 291   // offline->online (chip)\r
 292   for (UInt_t eq=0; eq<20; eq++) {\r
 293     for (UInt_t iModule=0; iModule<12; iModule++) {\r
 294       if (GetModuleNumber(eq,iModule)==(Int_t)module) {\r
 295         if (module<80) {\r
 296           if (eq<10) { // side A\r
 297             return (159-colM)/32 + 5*(iModule%2);\r
 298           }\r
 299           else { // side C\r
 300             return colM/32 + 5*(iModule%2);\r
 301           }\r
 302         }\r
 303         else if (module<240) {\r
 304           if (eq<10) { // side A\r
 305             return colM/32 + 5*(iModule%2);\r
 306           }\r
 307           else { // side C\r
 308             return (159-colM)/32 + 5*(iModule%2);\r
 309           }\r
 310         }\r
 311       }\r
 312     }\r
 313   }\r
 314   return 10; // error\r
 315 }\r
 316 //__________________________________________________________________________\r
 317 Int_t AliAnalysisTaskSPD::GetModuleNumber(UInt_t iDDL, UInt_t iModule) {\r
 318   if (iDDL<20 && iModule<12) return fgkDDLModuleMap[iDDL][iModule];\r
 319   else return 240;\r
 320 }\r
 321 \r
 322 \r