MUON/AliMUONClusterInput.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 /* $Id$ */
  17
  18 #include <TClonesArray.h>
  19 #include <TMinuit.h>
  20
  21 #include "AliRun.h"
  22 #include "AliMUON.h"
  23 #include "AliMUONChamber.h"
  24 #include "AliMUONConstants.h"
  25 #include "AliMUONClusterInput.h"
  26 #include "AliMUONMathieson.h"
  27 #include "AliMUONRawCluster.h"
  28 #include "AliMUONDigit.h"
  29 #include "AliLog.h"
  30
  31 ClassImp(AliMUONClusterInput)
  32
  33 AliMUONClusterInput* AliMUONClusterInput::fgClusterInput = 0;
  34 TMinuit* AliMUONClusterInput::fgMinuit = 0;
  35 AliMUONMathieson* AliMUONClusterInput::fgMathieson = 0;
  36
  37 AliMUONClusterInput::AliMUONClusterInput()
  38   : TObject(),
  39     fCluster(0),
  40     fChargeCorrel(1.),
  41     fSegmentationType(1),
  42     fDetElemId(0)
  43
  44 {
  45   fDigits[0]=0;
  46   fDigits[1]=0;
  47   fSegmentation[0]=0;
  48   fSegmentation[1]=0;
  49   fSegmentation2[0]=0;
  50   fSegmentation2[1]=0;
  51 }
  52
  53 AliMUONClusterInput* AliMUONClusterInput::Instance()
  54 {
  55 // return pointer to the singleton instance
  56     if (fgClusterInput == 0) {
  57         fgClusterInput = new AliMUONClusterInput();
  58         fgMinuit = new TMinuit(8);
  59     }
  60
  61     return fgClusterInput;
  62 }
  63
  64 AliMUONClusterInput::~AliMUONClusterInput()
  65 {
  66 // Destructor
  67     delete fgMinuit;
  68     delete fgMathieson;
  69 }
  70
  71 AliMUONClusterInput::AliMUONClusterInput(const AliMUONClusterInput& clusterInput):TObject(clusterInput)
  72 {
  73 // Protected copy constructor
  74
  75   AliFatal("Not implemented.");
  76 }
  77
  78 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig1, TClonesArray* dig2)
  79 {
  80 // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
  81     fChamber=chamber;
  82     fDigits[0]=dig1;
  83     fDigits[1]=dig2;
  84     fNDigits[0]=dig1->GetEntriesFast();
  85     fNDigits[1]=dig2->GetEntriesFast();
  86
  87     AliMUON *pMUON;
  88     AliMUONChamber* iChamber;
  89
  90     pMUON = (AliMUON*) gAlice->GetModule("MUON");
  91     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
  92       AliFatal("Wrong segmentation type");
  93
  94     iChamber =  &(pMUON->Chamber(chamber));
  95     fgMathieson = new AliMUONMathieson();
  96
  97     fSegmentation[0]=iChamber->SegmentationModel(1);
  98     fSegmentation[1]=iChamber->SegmentationModel(2);
  99
 100     fNseg = 2;
 101     if (chamber < AliMUONConstants::NTrackingCh()) {
 102       if (chamber > 1 ) {
 103         fgMathieson->SetPitch(AliMUONConstants::Pitch());
 104         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
 105         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
 106         fChargeCorrel = AliMUONConstants::ChargeCorrel();
 107       } else {
 108         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
 109         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
 110         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
 111         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
 112       }
 113     }
 114 }
 115
 116 void AliMUONClusterInput::SetDigits(Int_t chamber, TClonesArray* dig)
 117 {
 118 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
 119     fDigits[0]=dig;
 120     AliMUON *pMUON;
 121     AliMUONChamber* iChamber;
 122
 123     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 124     if ((fSegmentationType = pMUON->WhichSegmentation()) != 1)
 125       AliFatal("Wrong segmentation type");
 126
 127     iChamber =  &(pMUON->Chamber(chamber));
 128     fSegmentation[0]=iChamber->SegmentationModel(1);
 129
 130     fNseg=1;
 131 }
 132 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig1, TClonesArray* dig2)
 133 {
 134   // Set pointer to digits with corresponding segmentations and responses (two cathode planes)
 135     fChamber = chamber;
 136     fDetElemId = idDE;
 137     fDigits[0]  = dig1;
 138     fDigits[1]  = dig2;
 139     fNDigits[0] = dig1->GetEntriesFast();
 140     fNDigits[1] = dig2->GetEntriesFast();
 141
 142     AliMUON *pMUON;
 143     AliMUONChamber* iChamber;
 144
 145     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 146     iChamber =  &(pMUON->Chamber(chamber));
 147
 148     fgMathieson = new AliMUONMathieson();
 149     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
 150       AliFatal("Wrong segmentation type");
 151
 152     fSegmentation2[0]=iChamber->SegmentationModel2(1);
 153     fSegmentation2[1]=iChamber->SegmentationModel2(2);
 154
 155     fNseg = 2;
 156     if (chamber < AliMUONConstants::NTrackingCh()) {
 157       if (chamber > 1 ) {
 158         fgMathieson->SetPitch(AliMUONConstants::Pitch());
 159         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3());
 160         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3());
 161         fChargeCorrel = AliMUONConstants::ChargeCorrel();
 162       } else {
 163         fgMathieson->SetPitch(AliMUONConstants::PitchSt1());
 164         fgMathieson->SetSqrtKx3AndDeriveKx2Kx4(AliMUONConstants::SqrtKx3St1());
 165         fgMathieson->SetSqrtKy3AndDeriveKy2Ky4(AliMUONConstants::SqrtKy3St1());
 166         fChargeCorrel = AliMUONConstants::ChargeCorrelSt1();
 167       }
 168     }
 169 }
 170
 171 void AliMUONClusterInput::SetDigits(Int_t chamber, Int_t idDE, TClonesArray* dig)
 172 {
 173 // Set pointer to digits with corresponding segmentations and responses (one cathode plane)
 174
 175     fChamber = chamber;
 176     fDetElemId = idDE;
 177     fDigits[0] = dig;
 178
 179     AliMUON *pMUON;
 180     AliMUONChamber* iChamber;
 181
 182     pMUON = (AliMUON*) gAlice->GetModule("MUON");
 183     iChamber =  &(pMUON->Chamber(chamber));
 184     if ((fSegmentationType = pMUON->WhichSegmentation()) != 2)
 185       AliFatal("Wrong segmentation type");
 186
 187     fSegmentation2[0]=iChamber->SegmentationModel2(1);
 188
 189     fNseg=1;
 190 }
 191
 192 void  AliMUONClusterInput::SetCluster(AliMUONRawCluster* cluster)
 193 {
 194 // Set the current cluster
 195     //PH printf("\n %p \n", cluster);
 196     fCluster=cluster;
 197     Float_t qtot;
 198     Int_t   i, cath, ix, iy;
 199     AliMUONDigit* digit;
 200     fNmul[0]=cluster->GetMultiplicity(0);
 201     fNmul[1]=cluster->GetMultiplicity(1);
 202     //PH printf("\n %p %p ", fDigits[0], fDigits[1]);
 203
 204     for (cath=0; cath<2; cath++) {
 205         qtot=0;
 206         for (i=0; i<fNmul[cath]; i++) {
 207             // pointer to digit
 208             digit =(AliMUONDigit*)
 209                 (fDigits[cath]->UncheckedAt(cluster->GetIndex(i,cath)));
 210             // pad coordinates
 211             ix = digit->PadX();
 212             iy = digit->PadY();
 213             // pad charge
 214             fCharge[i][cath] = digit->Signal();
 215             // pad centre coordinates
 216 //          fSegmentation[cath]->GetPadCxy(ix, iy, x, y);
 217             // globals kUsed in fitting functions
 218             fix[i][cath]=ix;
 219             fiy[i][cath]=iy;
 220             // total charge per cluster
 221             qtot+=fCharge[i][cath];
 222             // Current z
 223             Float_t xc, yc;
 224             if (fSegmentationType == 1)
 225               fSegmentation[cath]->GetPadC(ix,iy,xc,yc,fZ);
 226             else
 227               fSegmentation2[cath]->GetPadC(fDetElemId,ix,iy,xc,yc,fZ);
 228
 229         } // loop over cluster digits
 230         fQtot[cath]=qtot;
 231         fChargeTot[cath]=Int_t(qtot);
 232     }  // loop over cathodes
 233 }
 234
 235
 236
 237 Float_t AliMUONClusterInput::DiscrChargeS1(Int_t i,Double_t *par)
 238 {
 239 // Compute the charge on first cathod only.
 240 return DiscrChargeCombiS1(i,par,0);
 241 }
 242
 243 Float_t AliMUONClusterInput::DiscrChargeCombiS1(Int_t i,Double_t *par, Int_t cath)
 244 {
 245 // par[0]    x-position of cluster
 246 // par[1]    y-position of cluster
 247
 248     Float_t q1;
 249     if (fSegmentationType == 1) {
 250
 251       fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
 252       //  First Cluster
 253       fSegmentation[cath]->SetHit(par[0],par[1],fZ);
 254       q1 = fgMathieson->IntXY(fSegmentation[cath]);
 255
 256     } else {
 257
 258       fSegmentation2[cath]->SetPad(fDetElemId, fix[i][cath], fiy[i][cath]);
 259       //  First Cluster
 260       fSegmentation2[cath]->SetHit(fDetElemId, par[0],par[1],fZ);
 261       q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 262    }
 263
 264    Float_t value = fQtot[cath]*q1;
 265    return value;
 266 }
 267
 268
 269 Float_t AliMUONClusterInput::DiscrChargeS2(Int_t i,Double_t *par)
 270 {
 271 // par[0]    x-position of first  cluster
 272 // par[1]    y-position of first  cluster
 273 // par[2]    x-position of second cluster
 274 // par[3]    y-position of second cluster
 275 // par[4]    charge fraction of first  cluster
 276 // 1-par[4]  charge fraction of second cluster
 277
 278   Float_t q1, q2;
 279
 280   if (fSegmentationType == 1) {
 281
 282     fSegmentation[0]->SetPad(fix[i][0], fiy[i][0]);
 283     //  First Cluster
 284     fSegmentation[0]->SetHit(par[0],par[1],fZ);
 285     q1 = fgMathieson->IntXY(fSegmentation[0]);
 286
 287     //  Second Cluster
 288     fSegmentation[0]->SetHit(par[2],par[3],fZ);
 289     q2 = fgMathieson->IntXY(fSegmentation[0]);
 290
 291   } else {
 292
 293     fSegmentation2[0]->SetPad(fDetElemId, fix[i][0], fiy[i][0]);
 294     //  First Cluster
 295     fSegmentation2[0]->SetHit(fDetElemId, par[0],par[1],fZ);
 296     q1 = fgMathieson->IntXY(fSegmentation[0]);
 297
 298     //  Second Cluster
 299     fSegmentation2[0]->SetHit(fDetElemId,par[2],par[3],fZ);
 300     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[0]);
 301   }
 302
 303   Float_t value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
 304   return value;
 305 }
 306
 307 Float_t AliMUONClusterInput::DiscrChargeCombiS2(Int_t i,Double_t *par, Int_t cath)
 308 {
 309 // par[0]    x-position of first  cluster
 310 // par[1]    y-position of first  cluster
 311 // par[2]    x-position of second cluster
 312 // par[3]    y-position of second cluster
 313 // par[4]    charge fraction of first  cluster - first cathode
 314 // 1-par[4]  charge fraction of second cluster
 315 // par[5]    charge fraction of first  cluster - second cathode
 316
 317   Float_t q1, q2;
 318
 319   if (fSegmentationType == 1) {
 320
 321     fSegmentation[cath]->SetPad(fix[i][cath], fiy[i][cath]);
 322     //  First Cluster
 323     fSegmentation[cath]->SetHit(par[0],par[1],fZ);
 324     q1 = fgMathieson->IntXY(fSegmentation[cath]);
 325
 326     //  Second Cluster
 327     fSegmentation[cath]->SetHit(par[2],par[3],fZ);
 328     q2 = fgMathieson->IntXY(fSegmentation[cath]);
 329
 330   } else {
 331     fSegmentation2[cath]->SetPad(fDetElemId,fix[i][cath], fiy[i][cath]);
 332     //  First Cluster
 333     fSegmentation2[cath]->SetHit(fDetElemId,par[0],par[1],fZ);
 334     q1 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 335
 336     //  Second Cluster
 337     fSegmentation2[cath]->SetHit(fDetElemId,par[2],par[3],fZ);
 338     q2 = fgMathieson->IntXY(fDetElemId, fSegmentation2[cath]);
 339   }
 340   Float_t value;
 341   if (cath==0) {
 342     value = fQtot[0]*(par[4]*q1+(1.-par[4])*q2);
 343   } else {
 344     value = fQtot[1]*(par[5]*q1+(1.-par[5])*q2);
 345   }
 346   return value;
 347 }
 348
 349 AliMUONClusterInput& AliMUONClusterInput
 350 ::operator = (const AliMUONClusterInput& rhs)
 351 {
 352 // Protected assignement operator
 353
 354   if (this == &rhs) return *this;
 355
 356   AliFatal("Not implemented.");
 357
 358   return *this;
 359 }